BACKGROUND
breast cancer  is the most frequent cancer in women worldwide  <cit> , and its development and progression are known to rely strongly on the stimulation by female sexual hormones, especially estrogens.
estrogen receptors α  and β  are transcription factors that mediate the actions of estrogens in target cells  <cit> .
ligand binding to erα or erβ induces receptor dimerization, either as homodimers  or heterodimers   <cit> , and promotes its translocation to the nucleus and binding of target chromatin sites via estrogen response elements  and other regulatory elements on dna  <cit> .
although encoded by different genes, erα and erβ share the same general modular protein structure of the nuclear hormone receptor superfamily, and have almost 100% amino acid sequence homology in their dna-binding domain.
they also show 59% amino acid homology in their ligand-binding domains  <cit> .
the two ers are thus quite similar in sequence and structure, but erβ has considerably different and, in most cases, opposite biological effects compared to erα in bc cells, both in vivo and in vitro.
erα regulates the transcription of hundreds of genes  <cit> , enhancing bc cell growth, proliferation and survival in response to estrogens  <cit> .
the specific role of erβ and its impact in bc are unclear.
this er subtype is expressed in 70% of human breast tumors in combination with erα, even if some human breast tumors express only erβ  <cit> .
several reports have suggested that erβ has anti-proliferative action in bc cells, by increasing the expression of anti-proliferative genes and/or decreasing the expression of proliferative and anti-apoptotic genes  <cit>  and the erα/erβ ratio determines the cell-specific response to estrogen.
in bc this ratio is higher than in normal tissues, due to up-regulation of erα and down-regulation of erβ  <cit> .
loss of erβ mrna levels in cancer can occur as a result of promoter methylation  <cit> .
these observations have suggested a positive prognostic value of this receptor subtype  <cit> .
however, several studies have reported also a negative prognostic value for erβ expression  <cit> , making the overall contribution of this receptor isoform to bc biology unclear.
the exact mechanism of the antagonism between erβ and erα is only partially known.
as the two receptors share only 30% homology in their transactivation domain af <dig>  <cit> , it is likely that they show different patterns of interaction with coregulatory proteins.
indeed, we have previously reported that the protein interactomes of both erβ and erα show significant differences of the protein complexes engaged by the two er subtypes  <cit> .
a particularly interesting subset of interacting proteins, with only partially overlapping interaction patterns between the two receptors  <cit> , comprises factors involved in rna maturation and splicing   <cit> .
alternative splicing is a mechanism by which cells can increase the variability of their proteomes by changing the composition of transcribed genes through differential choice of exons to be included in the final mrna molecule  <cit> .
almost 90% of human genes show alternative splicing during development, cell differentiation and disease  <cit> .
recent studies have shown the existence of cancer-specific splicing events by which transformed cells switch from the adult isoform of the gene to a more embryonic one, contributing to the cancer phenotype  <cit> .
alternative splicing events have been monitored in bc and in numerous tumor types  <cit> , and erα itself has been reported to induce alternative splicing of a specific set of genes  <cit> .
here, we investigated the ability of erβ to regulate mrna maturation and splicing in hormone-responsive bc cells.
to this purpose, we performed high-throughput rna sequencing  analysis of human mcf- <dig> cell lines stably transfected with erβ, and compared them with the wild type line, expressing only erα, upon 17β-estradiol  stimulation.
RESULTS
high-throughput sequencing in erβ + and erβ- human bc cell lines
we have previously established and characterized subclones of the human bc cell line mcf- <dig> expressing human erβ fused to a tandem affinity purification  tag, and have shown that the addition of the tap-tag at either the n- or the c-terminus of the protein  does not alter significantly the receptor function, nor its ability to activate transcription or to antagonize erα-dependent transcription  <cit> .
in order to study the early events of hormone-induced pre-mrna maturation, we used these stable cell lines to perform deep-sequencing analysis of estrogen-induced transcriptional events shortly after stimulation with 17β-estradiol , to focus mostly on primary transcriptional events  <cit> .
for comparison, we also performed the experiment in wild-type mcf- <dig> cells, which do not express endogenous erβ.
almost  <dig> million reads/replicate were aligned against the reference human genome for erβ- and erβ + bc cell lines.
the number of reads for genes and isoforms were normalized to “fragment per kilobases of exon per million of mapped reads” .
in order to analyze genes and isoforms, we set  <dig>  fpkm  as the expression level threshold.
in this way, we identified  <dig>  ,  <dig>   and  <dig>   genes as expressed.
the criteria for considering genes and isoforms as significantly regulated by estradiol were: fpkm value ≥ <dig>  in at least one analyzed condition, q-value  ≤ <dig>  and |fold-change|  ≥ <dig>  .
as shown in figure  <dig>   <dig> ,  <dig>   and  <dig>   genes were detected as significantly regulated by e <dig> in these bc cell lines.
expression of erβ in mcf <dig> cells significantly affected the estrogen-dependent gene expression profile: the regulation of around  <dig> genes  was lost in both cell lines expressing erβ, while a large number of genes which were not regulated in wt cells became significantly regulated in ct-erβ  and nt-erβ  clones .
the genes regulated consistently in both erβ + lines are reported in additional file 2: tables s1-d. interestingly, expression of erβ in this bc cell line had a stronger effect on inhibited genes than on the activated ones: regulation of 40% of genes inhibited by estradiol in wt cells was lost in both erβ + cell lines, versus 14% for estrogen-activated genes.
gene ontology analysis revealed that among the most enriched functions in the group of genes whose regulation by estradiol was lost in both erβ + cells there were, as expected, dna replication, recombination and repair, as well as cell cycle and cell morphology .figure  <dig> venn diagram of expressed genes in erβ + and wt breast cancer cell lines.
the venn diagram shows the number of genes regulated by estradiol in the three cell lines, as indicated: wt ; ct-erβ ; nt-erβ .
the pie charts in the lower panels specify the direction of regulation by estradiol  of the genes whose regulation is present in wt cells but is lost in both erβ + cell lines , or of the genes whose regulation by estradiol is not present in the wt cells but appears in both erβ + lines .
estrogen-dependent splicing events in erβ + and erβ- human bc cell lines
events of exon skipping, mutually exclusive exons, alternative start, stop splice site and intron retention were annotated using the multivariate analysis of transcript splicing  software  <cit>  .
exon skipping appeared to be the predominant splice event in all cell lines analyzed.
mats reveled in detail  <dig>  ,  <dig>   and  <dig>  exon skipping events induced by estradiol, associated with  <dig>  ,  <dig>   and  <dig>  genes .
five hundred seventy-five events were common to all the cell lines analyzed while  <dig> showed opposite exon inclusion level in erβ + lines compared to erβ- wt cells.
we also observed high levels of retained intron and mutually exclusive exons events, confirming a complex and significant effect of ers on the regulation of rna splicing in these cell lines .figure  <dig> annotation of splice events in erα + and erα + erβ + bc cell lines.
the bar plot shows the number of all alternative splicing events occurring in the cell lines analyzed.
inclusion and exclusion behavior for each event are shown .
genes whose regulation has opposite direction in the erβ + lines compared to the wt mcf- <dig>  the heat map on the right side shows the gene expression fold changes induced by estradiol.
the matrix on the left side shows in black those genes for which a splicing event was detected in at least one of the cell lines.
the nomenclature for the cell lines is the following: wt ; ct-erβ ; nt-erβ .
to focus on the differences in splicing patterns between erβ + and erβ- cell lines, we first looked at the genes which were regulated by estradiol in opposite direction in erβ + cells versus wt cells .
among  <dig> regulated genes,  <dig> also underwent estradiol-induced alternative splicing in at least one of the cell lines, confirming that pre-mrna maturation was regulated concurrently with transcription in a significant fraction of erβ-regulated genes .
these erβ-regulated genes undergoing alternative splicing included transcriptional regulators , enzymes involved in nucleoside/nucleotide metabolism , actin remodeling and cellular transport processes , and protein translation .
estrogen-induced differential splicing has been reported also in genes independently on transcriptional regulation  <cit> .
in order to find the splicing events differentially regulated by estradiol in erβ + compared to wt cells, we scanned the whole list of expressed genes for splicing patterns occurring differentially in erβ + vs erβ- cells.
in addition, we focused on those splicing events whose occurrence significantly altered the ratio between different isoforms of the same gene.
therefore, for each isoform we calculated the percentage of gene expression associated with that particular isoform  and selected for further analysis only those isoforms in which estradiol induced a change in fpkm ratio of at least 10% either in wt mcf- <dig> or in both erβ + cell lines.
other criteria of inclusion were: occurrence of at least one splicing event as detected by mats; at least one isoform of the gene with fpkm value ≥ <dig>  in at least one analyzed condition, q-value ≤ <dig>  either in wt or in both erβ + cell lines; regulation in both erβ + lines in opposite directions compared to the wt.
in this way, we identified the quantitatively most relevant splicing events differentially regulated by e <dig> in erβ + versus erβ- cells, including  <dig> genes whose isoform composition changed significantly after e <dig> stimulation in an erβ-dependent fashion .
among these, we found genes involved in apoptosis , lipid metabolism , nutrient transport , transmembrane receptor signaling , notch signaling , as well as some non-coding rnas .
an example of a gene whose splicing pattern was affected by estradiol in an erβ-dependent fashion is reported in figure 4a, showing the gene sgk <dig>  which encodes a serum and glucocorticoid-induced serine/threonine protein kinase involved in ion transport affecting many cellular processes such as cell growth, proliferation, survival, apoptosis and migration  <cit> .
the gene was induced by estradiol in both wt and erβ + cells; however, expression of erβ in the absence of hormonal stimulation induced a promoter usage switch and retention of the first intron causing an alternative translation start site and therefore the expression of an isoform with a different n-terminal sequence , compared to the major isoform expressed .
it has been reported that alternative isoforms at the n-terminus can affect sgk <dig> localization and protein stability: the erβ-specific form with intron-retention misses a very crucial sequence involved in targeting to the endoplasmic reticulum as well as in proteasomal degradation  <cit> .figure  <dig> selected isoform switches affected by the expression of erβ.
the erβ-dependent differential splicing events that affect most prominently the balance of different isoforms for each gene were identified by using the following parameters:  isoform ratio  changing of at least 10% after estradiol stimulation, either in the wt cells or in both the erβ + cells in at least one of the isoforms of the gene;  at least one isoform of the gene significantly regulated by estradiol  either in the wt cells or in both the erβ + cells;  isoform ratio changing in opposite direction in the wt cells compared to the erβ + cells;  at least one splicing event identified by mats analysis.
thirty-five genes satisfied all the selection requirements; of these, only  <dig> to  <dig> isoforms were selected for presentation, according to the expression levels and the regulation by estradiol.
left panel: heat map of regulation of the selected genes by estradiol in: ct-erβ ; nt-erβ ; wt .
right panel: heat map with the change in isoform ratio  between e2-treated and non-treated cells in the indicated cell lines; for each gene, at least two different isoforms are presented, to show estrogen-induced switch from one isoform to the other.
this data suggests that expression of erβ causes switches in estradiol-induced splicing patterns, potentially affecting expression or function or er targets.
estrogen-dependent alternative promoter usage in erβ + vs erβ- bc cells
as the usage of alternative promoters is a major determinant of protein diversity, even more than alternative splicing  <cit> , we next focused on utilization of multiple promoters.
we grouped the primary transcripts of a gene based on the promoter used, and subsequently tested changes in primary transcript abundance by measuring the square root of the jensen-shannon divergence that occurred within and between the analyzed groups.
finally, we investigated the potential promoter switch regulation for the genes comprising more than one differentially expressed transcript initiating from distinct genomic loci.
there were  <dig> ,  <dig>  and  <dig>  distinct promoter-switching genes  in erβ + and erβ- bc cell lines, respectively .
of the  <dig> promoter-switching genes recorded in wt mcf- <dig> cells,  <dig> did not show promoter switch in both erβ + cell lines, while  <dig> new promoter-switching events, not present in wt cells, were detected in both erβ + lines.
these  <dig> erβ-specific promoter-switching genes are involved in important cellular functions known to be controlled by e <dig> in bc cells, such as transcription , dna metabolism and repair , pre-mrna maturation and splicing , translation , protein ubiquitination and proteasome pathway , cytoskeleton and cytokinesis , membrane metabolism, remodeling and intracellular transport , cell adhesion and polarity , signal transduction .
as an example, we focused on the psd <dig> gene, predicted to be a nucleotide exchange factor for adp ribosylation factor   <dig>  a member of the ras family involved in vesicular trafficking, remodeling of membrane lipids, and signal transduction  <cit> .
as show in figure 4b, for this gene we found  <dig> distinct primary transcripts  whose usage ratios changed with e <dig> stimulus  compared to the control untreated cells .
in erβ + cells, estradiol induced a switch from promoter tss <dig>  to the downstream promoter tss <dig>  resulting in a shorter transcript , which is predicted to undergo nonsense-mediated decay .
estrogen-dependent splice ratios in erβ + vs erβ- cells
to obtain a comprehensive view of the estrogen-induced differences in the splice ratios between wt, ct-erβ and nt-erβ cells, we employed cuffdiff v <dig> .
<dig>  <cit> , which calculates the changes in splice isoforms abundance, by quantifying the square root of the jensen-shannon divergence, considering each primary transcripts able to produce multiple isoforms.
we determined  <dig> ,  <dig>  and  <dig>  differentially spliced genes  with a fdr value < <dig>  .
of the  <dig> spliced genes detected in wt cells, sixty-nine lacked splicing events in both erβ + cell lines, suggesting inhibition of erα-dependent splicing by erβ.
moreover,  <dig> erβ-specific dsgs, in which no estrogen-induced splicing was recorded in wt cells, were found in both erβ + cell lines.
these include genes involved in mitosis and cytokinesis , cell cycle control , transcription , protein folding , ubiquitination and sumoylation , and signal transduction .
correlation between er binding and er-dependent splicing
in order to identify direct splicing targets of erα and erβ, we next investigated the presence of erα and erβ binding sites in genomic locations close to the identified dsgs .
to verify the specificity of er binding sites, we created a heat-map based on previously published chip-seq data  <cit> .
the binding densities of erα and erβ were clustered according to the seqminer platform  <cit> .
in the clustering shown in , each line represents a genomic location of a binding site with its surrounding ± <dig>  kb region.
in the left panel, erα binding sites were used as a reference to collect a chip-seq tag densities window in erβ + and erβ- cell lines, while in the right panel erβ binding sites were used as a reference.
the heat map, representing the clustered density matrix, confirmed that erβ presence modified a significant number of erα binding sites.
in order to investigate the correlation between erα and erβ dna binding and splicing events, we compared our rna-seq data with the chip-seq data we had previously obtained in these same cell lines  <cit> .
we considered the binding sites included within regions spanning 10 kb upstream or downstream of all dsgs in each cell line.
the circos plot  <cit>  in figure  <dig> shows the dsgs for which estradiol challenge induced at least one splicing event in both erβ + cell lines with a nearby binding site for erα and/or erβ.
the outer ring  reports the erα binding sites, while the inner ring  shows the erβ binding sites.
based on erα and erβ binding, we distinguished three different dsgs groups.
group  <dig> dsgs were associated with both erα and erβ binding sites .
group  <dig> includes dsgs associated with erα binding sites exclusively  and group  <dig> dsgs associated with erβ binding sites exclusively .
the vast majority of genes exhibited both binding sites , confirming the competing role of erα and erβ.
interestingly, among the putative direct erβ splicing targets we found many genes involved in transcription: transcription factors , histone methyltransferases  and acetyltransferase , other transcriptional regulators .
other interesting functions associated with these genes are: apoptosis , autophagy , splicing , protein ubiquitination  cytoskeleton and cytokinesis , primary cilium biogenesis  <cit>  , intracellular trafficking/membrane trafficking , cell adhesion .figure  <dig> correlation between erβ-specific splicing and er binding.
circos plot of differential spliced genes  common to both erβ + cell lines which contain at least one er binding event within a window of  <dig> kb around the gene.
the outer ring shows chromosome ideograms with the relative genes located in their respective chromosomal locations, with the following color code: the genes that have both erβ and erα binding are in black; the genes which have only erβ binding are in red; the genes which have only erα binding are in blue.
the two internal rings represent erα and erβ binding events in blue and red, respectively.
taken together, these analyses suggest the strong effect of ers binding on alternative splicing and confirming the key role of erβ in bc cells.
discussion
in this study we investigated for the first time the effects of erβ expression on estrogen-dependent pre-mrna maturation and splicing.
we used two different subclones of erα-positive and estrogen-dependent mcf- <dig> cells stably expressing erβ, and we found that expression of erβ in these cells significantly affected the estradiol-dependent early transcriptional program and splicing pattern.
in particular, introduction of erβ caused loss of regulation of 25% of the estrogen-regulated genes.
moreover, a comparison between the erβ + and the wt erβ- cells showed that expression of erβ caused loss of erα-induced promoter switching in 75% of the genes, and of erα-induced splicing in 72% of the genes.
besides affecting erα-dependent transcription and splicing, expression of erβ caused the appearance of erβ-specific estrogen-responsive transcription  , promoter-switching  and differential splicing  events, counting only the events that were present in both erβ + lines.
the biological functions of the erβ-specific genes included dna replication and repair, cell cycle, apoptosis and autophagy, dna transcription, lipid metabolism, membrane metabolism, intracellular trafficking, mrna maturation and translation, protein ubiquitination and sumoylation, cell signaling, confirming that a wide variety of cellular processes are affected by erβ.
based on our data, there are at least three different mechanisms by which erβ can affect erα-dependent transcription:  competition with erα for binding to target gene promoters, in the form of competitive binding or heterodimerization, which can alter the recruitment of coregulators: the comparison between chip-seq and rna-seq data showed that the majority of the primary targets identified in this experiment were directly targeted by both erα and erβ , as expected also from data in the literature  <cit> ;  gain of new binding sites that are not bound by erα alone: a subset of binding sites appeared in the erβ + cells but were not present in wt cells, again consistently with a previous report comparing erα and erβ genomic binding patterns in mcf- <dig> cells  <cit> ;  secondary effects: expression of erβ induced transcription and splicing of transcriptional regulators  and of splicing factors, which can in turn affect estrogen-induced transcription and pre-mrna maturation.
for instance, erβ induced promoter switching in the ppm1g gene, encoding a protein phosphatase responsible for dephosphorylation of pre-mrna splicing factors  <cit> ; the erβ-specific alternate protein isoform from this gene  has an additional  <dig> amino acids in the n-terminal catalytic domain compared to the main isoform , which are likely to modulate the enzymatic function.
even if the sole detection of splicing events does not give information on the biological consequences of erβ effects on rna splicing, it is tempting to speculate on the possible implications of erβ-dependent splicing events on bc cells.
changes in splicing patterns can affect biological processes by many different mechanisms, including gain-of-function or functional switches, altered cellular localization, dominant negative effect, changes in protein/mrna stability.
for instance, expression of erβ induced a promoter switch in the usp <dig> gene, resulting in a shorter transcript  compared to the major isoform expressed in wt cells .
this erβ-specific isoform is missing the n-terminal  <dig> amino acids, a region of the protein of critical importance for interaction with substrates.
this suggests the possibility of a switch in substrate affinity induced by erβ.
this is particularly interesting as usp <dig> is a deubiquitinating enzyme, responsible for removing ubiquitin chains from both the tumor suppressor p <dig> and its negative regulator mdm <dig>  <cit> , which instead is an ubiquitin ligase inducing degradation of p <dig>  as usp <dig> binds both p <dig> and mdm <dig> with the same n-terminal domain  <cit> , the overall effect of its enzymatic activity is highly dependent on the relative affinity for the two targets.
the erβ-induced switch may alter this equilibrium, thus modulating such a relevant aspect of cancer biology as p <dig> stability.
in the case of ifnγ receptor  <dig>  erα induced expression of the full length transcript , while erβ favored a switch toward truncated forms  lacking the transmembrane domain, and therefore predicted to be secreted as soluble forms in the extracellular environment, with the potential of acting as dominant negative modulators of interferon signaling.
in another example , the erβ-induced splicing switch did not alter the open reading frame but caused the expression of a longer mrna  compared to the isoform that was favored by erα , possibly affecting the rate of translation or stability of the mrna.
another way of finding hints to the biological consequences of erβ-specific splicing events described here is the reconstruction of pathways whose genes are differentially spliced following erβ introduction in the cells.
for instance, we found many genes involved in the notch signaling pathway differentially spliced by erβ: the above-mentioned psen <dig> is the catalytic subunit of the γ-secretase complex, responsible for intramembrane proteolysis of transmembrane receptors including notch, resulting in the release of notch intracellular domain , which migrates to the nucleus and regulates transcription of target genes  <cit> .
also, ncor <dig> is bound to the unliganded cbf- <dig> transcription factor, a primary effector of notch signaling which acts as a repressor in unstimulated cells, but is converted to an activator  after binding the nicd  <cit> .
furthermore, poglut <dig> has recently been shown to be an endoplasmic reticulum o-glycosyl-transferase responsible for glycosylation of notch and required for its function  <cit> , and slc35c <dig> is an endoplasmic reticulum transporter responsible for accumulation of gdp-fucose, which is used for notch fucosylation, required for full activation  <cit> .
finally, sgk <dig> has recently been shown to be a negative regulator of notch signaling by inhibiting γ-secretase activity and promoting notch degradation  <cit> , and magi <dig> has been shown to recruit notch ligand dll <dig> to cadherin-based adherens junctions, stabilizing it on the cell surface  <cit> .
it is worth noting that the predominant notch receptor expressed in these cell lines was notch <dig>  which was induced by e <dig> in wt cells , while its basal expression was lowered and its up-regulation smoothened in the erβ + cells.
a possible modulation of the notch pathway by erβ is especially interesting as notch is a known regulator of breast development and maintenance of breast stem cells  <cit> ; alterations in the notch pathway have been involved in breast carcinogenesis, and in particular the notch pathway has been implicated in the development of triple negative bc , a particularly aggressive form of bc which does not express erα, progesterone receptor  or her <dig>  and which has shown resistance to all known therapies  <cit> .
targeting notch signaling has been proposed in tnbc.
as these cancers are erα-negative, hormonal treatment is not currently used for these patients; however, erβ could be expressed in up to 50% of tnbcs  <cit> , and its expression in tnbc has been associated with better prognosis  <cit> .
therefore, erβ may represent a potential new therapeutic target in tnbc.
the interrelation between erβ and notch in the development and prognosis of triple-negative bc should be investigated further in future research.
CONCLUSIONS
in conclusion, whole-genome analysis of early transcription evens and mrna processing associated with erβ confirmed a relevant role for this receptor in modulating erα-dependent transcription and splicing, but also identified novel, erβ-specific transcription and splicing events, confirming a wide range of actions of erβ in the biology of bc.
the data reported here confirm the complexity of estrogen action in bc cells and provide a comprehensive description of the effects of erα and by erβ on early transcription and splicing in hormone-responsive bc cells.
more importantly, they provide a starting point to identify the events of erβ-dependent splicing which are most significant for cancer biology.

 breast cancer is the most frequent cancer in women worldwide, and its development and progression are known to rely strongly on the stimulation by female sexual hormones, especially estrogens. estrogen receptors α and β are transcription factors that mediate the actions of estrogens in target cells. ligand binding to erα or erβ induces receptor dimerization, either as homodimers or heterodimers, and promotes its translocation to the nucleus and binding of target chromatin sites via estrogen response elements and other regulatory elements on dna. although encoded by different genes, erα and erβ share the same general modular protein structure of the nuclear hormone receptor superfamily, and have almost 100% amino acid sequence homology in their dna-binding domain. they also show 59% amino acid homology in their ligand-binding domains. the two ers are thus quite similar in sequence and structure, but erβ has considerably different and, in most cases, opposite biological effects compared to erα in bc cells, both in vivo and in vitro. erα regulates the transcription of hundreds of genes, enhancing bc cell growth, proliferation and survival in response to estrogens. the specific role of erβ and its impact in bc are unclear. this er subtype is expressed in 70% of human breast tumors in combination with erα, even if some human breast tumors express only erβ. several reports have suggested that erβ has anti-proliferative action in bc cells, by increasing the expression of anti-proliferative genes and/or decreasing the expression of proliferative and anti-apoptotic genes and the erα/erβ ratio determines the cell-specific response to estrogen. in bc this ratio is higher than in normal tissues, due to up-regulation of erα and down-regulation of erβ. loss of erβ mrna levels in cancer can occur as a result of promoter methylation. these observations have suggested a positive prognostic value of this receptor subtype. however, several studies have reported also a negative prognostic value for erβ expression, making the overall contribution of this receptor isoform to bc biology unclear. the exact mechanism of the antagonism between erβ and erα is only partially known. as the two receptors share only 30% homology in their transactivation domain af, it is likely that they show different patterns of interaction with coregulatory proteins. indeed, we have previously reported that the protein interactomes of both erβ and erα show significant differences of the protein complexes engaged by the two er subtypes. a particularly interesting subset of interacting proteins, with only partially overlapping interaction patterns between the two receptors, comprises factors involved in rna maturation and splicing. alternative splicing is a mechanism by which cells can increase the variability of their proteomes by changing the composition of transcribed genes through differential choice of exons to be included in the final mrna molecule. almost 90% of human genes show alternative splicing during development, cell differentiation and disease. recent studies have shown the existence of cancer-specific splicing events by which transformed cells switch from the adult isoform of the gene to a more embryonic one, contributing to the cancer phenotype. alternative splicing events have been monitored in bc and in numerous tumor types, and erα itself has been reported to induce alternative splicing of a specific set of genes. here, we investigated the ability of erβ to regulate mrna maturation and splicing in hormone-responsive bc cells. to this purpose, we performed high-throughput rna sequencing analysis of human mcf- cell lines stably transfected with erβ, and compared them with the wild type line, expressing only erα, upon 17β-estradiol stimulation.  high-throughput sequencing in erβ + and erβ- human bc cell lines we have previously established and characterized subclones of the human bc cell line mcf- expressing human erβ fused to a tandem affinity purification tag, and have shown that the addition of the tap-tag at either the n- or the c-terminus of the protein does not alter significantly the receptor function, nor its ability to activate transcription or to antagonize erα-dependent transcription. in order to study the early events of hormone-induced pre-mrna maturation, we used these stable cell lines to perform deep-sequencing analysis of estrogen-induced transcriptional events shortly after stimulation with 17β-estradiol, to focus mostly on primary transcriptional events. for comparison, we also performed the experiment in wild-type mcf- cells, which do not express endogenous erβ. almost million reads/replicate were aligned against the reference human genome for erβ- and erβ + bc cell lines. the number of reads for genes and isoforms were normalized to “fragment per kilobases of exon per million of mapped reads”. in order to analyze genes and isoforms, we set fpkm as the expression level threshold. in this way, we identified, and genes as expressed. the criteria for considering genes and isoforms as significantly regulated by estradiol were: fpkm value ≥ in at least one analyzed condition, q-value ≤ and |fold-change| ≥. as shown in figure , and genes were detected as significantly regulated by e in these bc cell lines. expression of erβ in mcf cells significantly affected the estrogen-dependent gene expression profile: the regulation of around genes was lost in both cell lines expressing erβ, while a large number of genes which were not regulated in wt cells became significantly regulated in ct-erβ and nt-erβ clones. the genes regulated consistently in both erβ + lines are reported in additional file 2: tables s1-d. interestingly, expression of erβ in this bc cell line had a stronger effect on inhibited genes than on the activated ones: regulation of 40% of genes inhibited by estradiol in wt cells was lost in both erβ + cell lines, versus 14% for estrogen-activated genes. gene ontology analysis revealed that among the most enriched functions in the group of genes whose regulation by estradiol was lost in both erβ + cells there were, as expected, dna replication, recombination and repair, as well as cell cycle and cell morphology.figure venn diagram of expressed genes in erβ + and wt breast cancer cell lines. the venn diagram shows the number of genes regulated by estradiol in the three cell lines, as indicated: wt; ct-erβ; nt-erβ. the pie charts in the lower panels specify the direction of regulation by estradiol of the genes whose regulation is present in wt cells but is lost in both erβ + cell lines, or of the genes whose regulation by estradiol is not present in the wt cells but appears in both erβ + lines. estrogen-dependent splicing events in erβ + and erβ- human bc cell lines events of exon skipping, mutually exclusive exons, alternative start, stop splice site and intron retention were annotated using the multivariate analysis of transcript splicing software. exon skipping appeared to be the predominant splice event in all cell lines analyzed. mats reveled in detail, and exon skipping events induced by estradiol, associated with, and genes. five hundred seventy-five events were common to all the cell lines analyzed while showed opposite exon inclusion level in erβ + lines compared to erβ- wt cells. we also observed high levels of retained intron and mutually exclusive exons events, confirming a complex and significant effect of ers on the regulation of rna splicing in these cell lines.figure annotation of splice events in erα + and erα + erβ + bc cell lines. the bar plot shows the number of all alternative splicing events occurring in the cell lines analyzed. inclusion and exclusion behavior for each event are shown. genes whose regulation has opposite direction in the erβ + lines compared to the wt mcf- the heat map on the right side shows the gene expression fold changes induced by estradiol. the matrix on the left side shows in black those genes for which a splicing event was detected in at least one of the cell lines. the nomenclature for the cell lines is the following: wt; ct-erβ; nt-erβ. to focus on the differences in splicing patterns between erβ + and erβ- cell lines, we first looked at the genes which were regulated by estradiol in opposite direction in erβ + cells versus wt cells. among regulated genes, also underwent estradiol-induced alternative splicing in at least one of the cell lines, confirming that pre-mrna maturation was regulated concurrently with transcription in a significant fraction of erβ-regulated genes. these erβ-regulated genes undergoing alternative splicing included transcriptional regulators, enzymes involved in nucleoside/nucleotide metabolism, actin remodeling and cellular transport processes, and protein translation. estrogen-induced differential splicing has been reported also in genes independently on transcriptional regulation. in order to find the splicing events differentially regulated by estradiol in erβ + compared to wt cells, we scanned the whole list of expressed genes for splicing patterns occurring differentially in erβ + vs erβ- cells. in addition, we focused on those splicing events whose occurrence significantly altered the ratio between different isoforms of the same gene. therefore, for each isoform we calculated the percentage of gene expression associated with that particular isoform and selected for further analysis only those isoforms in which estradiol induced a change in fpkm ratio of at least 10% either in wt mcf- or in both erβ + cell lines. other criteria of inclusion were: occurrence of at least one splicing event as detected by mats; at least one isoform of the gene with fpkm value ≥ in at least one analyzed condition, q-value ≤ either in wt or in both erβ + cell lines; regulation in both erβ + lines in opposite directions compared to the wt. in this way, we identified the quantitatively most relevant splicing events differentially regulated by e in erβ + versus erβ- cells, including genes whose isoform composition changed significantly after e stimulation in an erβ-dependent fashion. among these, we found genes involved in apoptosis, lipid metabolism, nutrient transport, transmembrane receptor signaling, notch signaling, as well as some non-coding rnas. an example of a gene whose splicing pattern was affected by estradiol in an erβ-dependent fashion is reported in figure 4a, showing the gene sgk which encodes a serum and glucocorticoid-induced serine/threonine protein kinase involved in ion transport affecting many cellular processes such as cell growth, proliferation, survival, apoptosis and migration. the gene was induced by estradiol in both wt and erβ + cells; however, expression of erβ in the absence of hormonal stimulation induced a promoter usage switch and retention of the first intron causing an alternative translation start site and therefore the expression of an isoform with a different n-terminal sequence, compared to the major isoform expressed. it has been reported that alternative isoforms at the n-terminus can affect sgk localization and protein stability: the erβ-specific form with intron-retention misses a very crucial sequence involved in targeting to the endoplasmic reticulum as well as in proteasomal degradation.figure selected isoform switches affected by the expression of erβ. the erβ-dependent differential splicing events that affect most prominently the balance of different isoforms for each gene were identified by using the following parameters: isoform ratio changing of at least 10% after estradiol stimulation, either in the wt cells or in both the erβ + cells in at least one of the isoforms of the gene; at least one isoform of the gene significantly regulated by estradiol either in the wt cells or in both the erβ + cells; isoform ratio changing in opposite direction in the wt cells compared to the erβ + cells; at least one splicing event identified by mats analysis. thirty-five genes satisfied all the selection requirements; of these, only to isoforms were selected for presentation, according to the expression levels and the regulation by estradiol. left panel: heat map of regulation of the selected genes by estradiol in: ct-erβ; nt-erβ; wt. right panel: heat map with the change in isoform ratio between e2-treated and non-treated cells in the indicated cell lines; for each gene, at least two different isoforms are presented, to show estrogen-induced switch from one isoform to the other. this data suggests that expression of erβ causes switches in estradiol-induced splicing patterns, potentially affecting expression or function or er targets. estrogen-dependent alternative promoter usage in erβ + vs erβ- bc cells as the usage of alternative promoters is a major determinant of protein diversity, even more than alternative splicing, we next focused on utilization of multiple promoters. we grouped the primary transcripts of a gene based on the promoter used, and subsequently tested changes in primary transcript abundance by measuring the square root of the jensen-shannon divergence that occurred within and between the analyzed groups. finally, we investigated the potential promoter switch regulation for the genes comprising more than one differentially expressed transcript initiating from distinct genomic loci. there were, and distinct promoter-switching genes in erβ + and erβ- bc cell lines, respectively. of the promoter-switching genes recorded in wt mcf- cells, did not show promoter switch in both erβ + cell lines, while new promoter-switching events, not present in wt cells, were detected in both erβ + lines. these erβ-specific promoter-switching genes are involved in important cellular functions known to be controlled by e in bc cells, such as transcription, dna metabolism and repair, pre-mrna maturation and splicing, translation, protein ubiquitination and proteasome pathway, cytoskeleton and cytokinesis, membrane metabolism, remodeling and intracellular transport, cell adhesion and polarity, signal transduction. as an example, we focused on the psd gene, predicted to be a nucleotide exchange factor for adp ribosylation factor a member of the ras family involved in vesicular trafficking, remodeling of membrane lipids, and signal transduction. as show in figure 4b, for this gene we found distinct primary transcripts whose usage ratios changed with e stimulus compared to the control untreated cells. in erβ + cells, estradiol induced a switch from promoter tss to the downstream promoter tss resulting in a shorter transcript, which is predicted to undergo nonsense-mediated decay. estrogen-dependent splice ratios in erβ + vs erβ- cells to obtain a comprehensive view of the estrogen-induced differences in the splice ratios between wt, ct-erβ and nt-erβ cells, we employed cuffdiff v., which calculates the changes in splice isoforms abundance, by quantifying the square root of the jensen-shannon divergence, considering each primary transcripts able to produce multiple isoforms. we determined, and differentially spliced genes with a fdr value <. of the spliced genes detected in wt cells, sixty-nine lacked splicing events in both erβ + cell lines, suggesting inhibition of erα-dependent splicing by erβ. moreover, erβ-specific dsgs, in which no estrogen-induced splicing was recorded in wt cells, were found in both erβ + cell lines. these include genes involved in mitosis and cytokinesis, cell cycle control, transcription, protein folding, ubiquitination and sumoylation, and signal transduction. correlation between er binding and er-dependent splicing in order to identify direct splicing targets of erα and erβ, we next investigated the presence of erα and erβ binding sites in genomic locations close to the identified dsgs. to verify the specificity of er binding sites, we created a heat-map based on previously published chip-seq data. the binding densities of erα and erβ were clustered according to the seqminer platform. in the clustering shown in, each line represents a genomic location of a binding site with its surrounding ±  kb region. in the left panel, erα binding sites were used as a reference to collect a chip-seq tag densities window in erβ + and erβ- cell lines, while in the right panel erβ binding sites were used as a reference. the heat map, representing the clustered density matrix, confirmed that erβ presence modified a significant number of erα binding sites. in order to investigate the correlation between erα and erβ dna binding and splicing events, we compared our rna-seq data with the chip-seq data we had previously obtained in these same cell lines. we considered the binding sites included within regions spanning 10 kb upstream or downstream of all dsgs in each cell line. the circos plot in figure  shows the dsgs for which estradiol challenge induced at least one splicing event in both erβ + cell lines with a nearby binding site for erα and/or erβ. the outer ring reports the erα binding sites, while the inner ring shows the erβ binding sites. based on erα and erβ binding, we distinguished three different dsgs groups. group dsgs were associated with both erα and erβ binding sites. group includes dsgs associated with erα binding sites exclusively and group dsgs associated with erβ binding sites exclusively. the vast majority of genes exhibited both binding sites, confirming the competing role of erα and erβ. interestingly, among the putative direct erβ splicing targets we found many genes involved in transcription: transcription factors, histone methyltransferases and acetyltransferase, other transcriptional regulators. other interesting functions associated with these genes are: apoptosis, autophagy, splicing, protein ubiquitination cytoskeleton and cytokinesis, primary cilium biogenesis, intracellular trafficking/membrane trafficking, cell adhesion.figure correlation between erβ-specific splicing and er binding. circos plot of differential spliced genes common to both erβ + cell lines which contain at least one er binding event within a window of kb around the gene. the outer ring shows chromosome ideograms with the relative genes located in their respective chromosomal locations, with the following color code: the genes that have both erβ and erα binding are in black; the genes which have only erβ binding are in red; the genes which have only erα binding are in blue. the two internal rings represent erα and erβ binding events in blue and red, respectively. taken together, these analyses suggest the strong effect of ers binding on alternative splicing and confirming the key role of erβ in bc cells. discussion in this study we investigated for the first time the effects of erβ expression on estrogen-dependent pre-mrna maturation and splicing. we used two different subclones of erα-positive and estrogen-dependent mcf- cells stably expressing erβ, and we found that expression of erβ in these cells significantly affected the estradiol-dependent early transcriptional program and splicing pattern. in particular, introduction of erβ caused loss of regulation of 25% of the estrogen-regulated genes. moreover, a comparison between the erβ + and the wt erβ- cells showed that expression of erβ caused loss of erα-induced promoter switching in 75% of the genes, and of erα-induced splicing in 72% of the genes. besides affecting erα-dependent transcription and splicing, expression of erβ caused the appearance of erβ-specific estrogen-responsive transcription, promoter-switching and differential splicing events, counting only the events that were present in both erβ + lines. the biological functions of the erβ-specific genes included dna replication and repair, cell cycle, apoptosis and autophagy, dna transcription, lipid metabolism, membrane metabolism, intracellular trafficking, mrna maturation and translation, protein ubiquitination and sumoylation, cell signaling, confirming that a wide variety of cellular processes are affected by erβ. based on our data, there are at least three different mechanisms by which erβ can affect erα-dependent transcription: competition with erα for binding to target gene promoters, in the form of competitive binding or heterodimerization, which can alter the recruitment of coregulators: the comparison between chip-seq and rna-seq data showed that the majority of the primary targets identified in this experiment were directly targeted by both erα and erβ, as expected also from data in the literature; gain of new binding sites that are not bound by erα alone: a subset of binding sites appeared in the erβ + cells but were not present in wt cells, again consistently with a previous report comparing erα and erβ genomic binding patterns in mcf- cells; secondary effects: expression of erβ induced transcription and splicing of transcriptional regulators and of splicing factors, which can in turn affect estrogen-induced transcription and pre-mrna maturation. for instance, erβ induced promoter switching in the ppm1g gene, encoding a protein phosphatase responsible for dephosphorylation of pre-mrna splicing factors; the erβ-specific alternate protein isoform from this gene has an additional amino acids in the n-terminal catalytic domain compared to the main isoform, which are likely to modulate the enzymatic function. even if the sole detection of splicing events does not give information on the biological consequences of erβ effects on rna splicing, it is tempting to speculate on the possible implications of erβ-dependent splicing events on bc cells. changes in splicing patterns can affect biological processes by many different mechanisms, including gain-of-function or functional switches, altered cellular localization, dominant negative effect, changes in protein/mrna stability. for instance, expression of erβ induced a promoter switch in the usp gene, resulting in a shorter transcript compared to the major isoform expressed in wt cells. this erβ-specific isoform is missing the n-terminal amino acids, a region of the protein of critical importance for interaction with substrates. this suggests the possibility of a switch in substrate affinity induced by erβ. this is particularly interesting as usp is a deubiquitinating enzyme, responsible for removing ubiquitin chains from both the tumor suppressor p and its negative regulator mdm, which instead is an ubiquitin ligase inducing degradation of p as usp binds both p and mdm with the same n-terminal domain, the overall effect of its enzymatic activity is highly dependent on the relative affinity for the two targets. the erβ-induced switch may alter this equilibrium, thus modulating such a relevant aspect of cancer biology as p stability. in the case of ifnγ receptor erα induced expression of the full length transcript, while erβ favored a switch toward truncated forms lacking the transmembrane domain, and therefore predicted to be secreted as soluble forms in the extracellular environment, with the potential of acting as dominant negative modulators of interferon signaling. in another example, the erβ-induced splicing switch did not alter the open reading frame but caused the expression of a longer mrna compared to the isoform that was favored by erα, possibly affecting the rate of translation or stability of the mrna. another way of finding hints to the biological consequences of erβ-specific splicing events described here is the reconstruction of pathways whose genes are differentially spliced following erβ introduction in the cells. for instance, we found many genes involved in the notch signaling pathway differentially spliced by erβ: the above-mentioned psen is the catalytic subunit of the γ-secretase complex, responsible for intramembrane proteolysis of transmembrane receptors including notch, resulting in the release of notch intracellular domain, which migrates to the nucleus and regulates transcription of target genes. also, ncor is bound to the unliganded cbf- transcription factor, a primary effector of notch signaling which acts as a repressor in unstimulated cells, but is converted to an activator after binding the nicd. furthermore, poglut has recently been shown to be an endoplasmic reticulum o-glycosyl-transferase responsible for glycosylation of notch and required for its function, and slc35c is an endoplasmic reticulum transporter responsible for accumulation of gdp-fucose, which is used for notch fucosylation, required for full activation. finally, sgk has recently been shown to be a negative regulator of notch signaling by inhibiting γ-secretase activity and promoting notch degradation, and magi has been shown to recruit notch ligand dll to cadherin-based adherens junctions, stabilizing it on the cell surface. it is worth noting that the predominant notch receptor expressed in these cell lines was notch which was induced by e in wt cells, while its basal expression was lowered and its up-regulation smoothened in the erβ + cells. a possible modulation of the notch pathway by erβ is especially interesting as notch is a known regulator of breast development and maintenance of breast stem cells; alterations in the notch pathway have been involved in breast carcinogenesis, and in particular the notch pathway has been implicated in the development of triple negative bc, a particularly aggressive form of bc which does not express erα, progesterone receptor or her and which has shown resistance to all known therapies. targeting notch signaling has been proposed in tnbc. as these cancers are erα-negative, hormonal treatment is not currently used for these patients; however, erβ could be expressed in up to 50% of tnbcs, and its expression in tnbc has been associated with better prognosis. therefore, erβ may represent a potential new therapeutic target in tnbc. the interrelation between erβ and notch in the development and prognosis of triple-negative bc should be investigated further in future research.  in conclusion, whole-genome analysis of early transcription evens and mrna processing associated with erβ confirmed a relevant role for this receptor in modulating erα-dependent transcription and splicing, but also identified novel, erβ-specific transcription and splicing events, confirming a wide range of actions of erβ in the biology of bc. the data reported here confirm the complexity of estrogen action in bc cells and provide a comprehensive description of the effects of erα and by erβ on early transcription and splicing in hormone-responsive bc cells. more importantly, they provide a starting point to identify the events of erβ-dependent splicing which are most significant for cancer biology. 
taken together, these results demonstrate that erβ significantly affects estrogen-induced early transcription and mrna splicing in hormone-responsive bc cells, providing novel information on the biological role of erβ in these tumors. guided by these findings, here we performed rna sequencing to investigate in depth the differences in the early transcriptional events and rna splicing patterns induced by estradiol in cells expressing erα alone or erα and erβ. estrogens play an important role in breast cancer development and progression; when the two isoforms of the estrogen receptor are co-expressed each of them mediate specific effects of these hormones in bc cells. we have previously described the erβ and erα interactomes from bc cells, identifying specific and distinct patterns of protein interactions for the two receptors. we identified several splicing events induced by erα alone and by erα + erβ, demonstrating for the first time that erβ significantly affects estrogen-induced splicing in bc cells, as revealed by modification of a subset of erα-dependent splicing by erβ, as well as by the presence of splicing isoforms only in erβ + cells. interestingly, we identified putative direct targets of erβ-mediated alternative splicing by correlating the genomic locations of erβ and erα binding sites with estradiol-induced differential splicing in the corresponding genes. as clinical evidence regarding a prognostic role for this receptor subtype in hormone-responsive bc is still limited and conflicting, more knowledge is required on the biological functions of erβ in cancer cells.
 estrogens play an important role in breast cancer development and progression; when the two isoforms of the estrogen receptor are co-expressed each of them mediate specific effects of these hormones in bc cells. erβ has been suggested to exert an antagonist role toward the oncogenic activities of erα, and for this reason it is considered an oncosuppressor. as clinical evidence regarding a prognostic role for this receptor subtype in hormone-responsive bc is still limited and conflicting, more knowledge is required on the biological functions of erβ in cancer cells. we have previously described the erβ and erα interactomes from bc cells, identifying specific and distinct patterns of protein interactions for the two receptors. in particular, we identified factors involved in mrna splicing and maturation as important components of both erα and erβ pathways. guided by these findings, here we performed rna sequencing to investigate in depth the differences in the early transcriptional events and rna splicing patterns induced by estradiol in cells expressing erα alone or erα and erβ.  exon skipping was the most abundant splicing event in the post-transcriptional regulation by estradiol. we identified several splicing events induced by erα alone and by erα + erβ, demonstrating for the first time that erβ significantly affects estrogen-induced splicing in bc cells, as revealed by modification of a subset of erα-dependent splicing by erβ, as well as by the presence of splicing isoforms only in erβ + cells. in particular, we observed that erβ + bc cell lines exhibited around 2-fold more splicing events than the erβ- cells. interestingly, we identified putative direct targets of erβ-mediated alternative splicing by correlating the genomic locations of erβ and erα binding sites with estradiol-induced differential splicing in the corresponding genes.  taken together, these results demonstrate that erβ significantly affects estrogen-induced early transcription and mrna splicing in hormone-responsive bc cells, providing novel information on the biological role of erβ in these tumors. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords breast cancerestrogen receptor betaalternative splicingalternative promotersrnaseqissue-copyright-statement© the author 2015 