Differential effects of follicle-stimulating hormone glycoforms on the transcriptome profile of cultured rat granulosa cells as disclosed by RNA-seq

It has been documented that variations in glycosylation on glycoprotein hormones, confer distinctly different biological features to the corresponding glycoforms when multiple in vitro biochemical readings are analyzed. We here applied next generation RNA sequencing to explore changes in the transcriptome of rat granulosa cells exposed for 0, 6, and 12 h to 100 ng/ml of four highly purified follicle-stimulating hormone (FSH) glycoforms, each exhibiting different glycosylation patterns: human pituitary FSH18/21 and equine FSH (eqFSH) (hypo-glycosylated), and human FSH24 and chinese-hamster ovary cell-derived human recombinant FSH (recFSH) (fully-glycosylated). Total RNA from triplicate incubations was prepared from FSH glycoform-exposed cultured granulosa cells obtained from DES-pretreated immature female rats, and RNA libraries were sequenced in a HighSeq 2500 sequencer (2 × 125 bp paired-end format, 10–15 × 106 reads/sample). The computational workflow focused on investigating differences among the four FSH glycoforms at three levels: gene expression, enriched biological processes, and perturbed pathways. Among the top 200 differentially expressed genes, only 4 (0.6%) were shared by all 4 glycoforms at 6 h, whereas 118 genes (40%) were shared at 12 h. Follicle-stimulating hormone glycocoforms stimulated different patterns of exclusive and associated up regulated biological processes in a glycoform and time-dependent fashion with more shared biological processes after 12 h of exposure and fewer treatment-specific ones, except for recFSH, which exhibited stronger responses with more specifically associated processes at this time. Similar results were found for down-regulated processes, with a greater number of processes at 6 h or 12 h, depending on the particular glycoform. In general, there were fewer downregulated than upregulated processes at both 6 h and 12 h, with FSH18/21 exhibiting the largest number of down-regulated associated processes at 6 h while eqFSH exhibited the greatest number at 12 h. Signaling cascades, largely linked to cAMP-PKA, MAPK, and PI3/AKT pathways were detected as differentially activated by the glycoforms, with each glycoform exhibiting its own molecular signature. These data extend previous observations demonstrating glycosylation-dependent differential regulation of gene expression and intracellular signaling pathways triggered by FSH in granulosa cells. The results also suggest the importance of individual FSH glycoform glycosylation for the conformation of the ligand-receptor complex and induced signalling pathways.


INTRODUCTION
Follicle-stimulating hormone is produced by the anterior pituitary gland (AP) as different glycoforms, defined by the presence or abscence of glycans on the hormone-specific β-subunit.
As other glycoprotein hormones, this gonadotropin is composed of two subunits, the  and  subunits, associated with each other through non-covalent interactions; the -subunit, is common to all glycoprotein hormones [luteinizing hormone (LH), choriogonadotropin (CG) and thyroid-stimulating hormone (TSH)], whereas the -subunit, confers specificity for binding and action to each gonadotropin at its cognate receptor [1].In human FSH (hFSH), four Asn residues, two in FSH (N52 and N78) and two in FSH (N7 and N24) are targets for Nlinked glycosylation [2] (Fig. 1).In glycoprotein hormones, the oligosaccharide in position Asn 52 is involved in the activation of the receptor/signal transducer (G protein) system and biological response [3][4][5][6], an effect mediated through the stabilization of the conformation of the FSH dimer bound to the FSH receptor (FSHR) [7][8][9].Meanwhile, glycans attached to FSH play a major role in defining the circulatory half-life and in vivo bioactivity of the gonadotropin [10], albeit it has been shown that they also impact on FSH-mediated intracellular signaling [11][12][13][14].
Human FSH macroheterogeneity occurs at either one or both FSH N-glycosylation sites (Fig. 1), determining differences in apparent molecular weights by gel electrophoresis and immunoblotting [15][16][17].In fact, Western blots of FSH recovered after gel electrophoresis, yielded two FSH bands, a 24 kDa band exhibiting both FSH-subunit Asn 7 and Asn 24 N-linked glycans (designated as fully-o tetra-glycosylated 24kDa-FSH) and a 21kDa band in which the Asn 24 glycan is absent (hypo-or tri-glycosylated 21kDa-FSH) [15,18], resulting in FSH heterodimers designated as FSH 24 and FSH 21 , respectively.Purified Hypo-glycosylated preparations include an additional hypo-glycosylated variant lacking the Asn 7 glycan, with its corresponding heterodimer designated as FSH 18 (18kDa-FSH (Fig. 1c), hence the designation FSH 18/21 .Although an additional hypo-glycosylated 15kDa-FSH variant also has been detected in human pituitaries, this form assembles poorly with FSH and consequently very low levels of secretion of the corresponding heterodimer (FSH 15 glycoform) are detected [19].Thus, three hFSH glycoforms (FSH 24 , FSH 21 , and FSH 18 ) appear to be the physiologically relevant FSH variants in humans.The eFSHα subunit has 4 additional residues at the N-terminus, accounting for the difference in numbering.The glycan at position β7 in eqFSH (black dotted square) is absent in the bulk (90%) of the molecules contained in highly purified preparations [27].Glycans in recFSH were taken from Mastrangeli et al [28].Human FSH glycoform models created with the GLYCAM web tool are shown on the right with the subunits rendered as cartoons using PyMol.
The same color scheme is employed for subunits and glycans, which are rendered as spheres.Partially obscured Asn 24 glycan is colored dark blue to distinquish it from Asn 7 glycan (light blue).
Similar to the influence of microheterogeneity (i.e.variations in the structure of the carbohydrates and complexity of the oligosaccharides on gonadotropins) on the ability of the gonadotropins to activate the FSHR and trigger intracellular signaling [20][21][22], FSH macroheterogeneity also contributes to its bioactivity [14,23,24].In fact, in vitro and in vivo studies have shown that pituitary and recombinant hFSH glycoforms display differential FSHR binding kinetics and bioactivity [12,14,25,26].Further, in a recent study employing cAMP accumulation, -arrestin-mediated ERK1/2 activation, and intracellular calcium (iCa2+) accumulation as read-outs in FSH-stimulated HEK-293 cells, we found that in addition to determining the intensity of the biological response at the target cell, the presence or abscence of glycans attached to FSH conferred some degree of biased agonism to the different FSH glycoforms [24].
Follicle-stimulating hormone stimulation triggers activation of a complex array of diverssignaling cascades mediated not only by the canonical Gs/cAMP/PKA pathway but also by other G proteins and receptor interacting proteins [29,30].Activation of this signaling network and particular signaling modules most probably occurs through stabilization of different FSHR conformations by FSH glycoforms possessing distinct glycosylation patters [26].Given the differential effects of FSH glycoforms on particular FSH-stimulated read outs [24], we here applied next generation sequencing (NGS) to explore the impact of differential glycosylation on the transcriptome of cultured rat granulosa cells exposed during different times to a fixed dose of four highly purified FSH glycoform preparations, each exhibiting different glycosylation patterns: tetra-and tri-glycosylated pituitary hFSH (FSH 24 and FSH 18/18/21 , respectively), equine FSH [eqFSH, 90% hypo-glycosylated], and Chinese-hamster ovary cell-produced human recombinant FSH (recFSH or follitropin-; 80% tetra-glycosylated) (Fig. 1).These four preparations also differ to varying extent in microheterogeneity [15,23,24,27,[31][32][33][34][35][36][37][38].Our analysis was focused on the effects of each FSH preparation at three levels: gene expression, enriched biological processes, and perturbed pathways.

MATERIAL AND METHODS Hormones
Human pituitary FSH 24 and FSH 18/21 were purified after extraction of human pituitary tissue (a generous gift of Dr. James A. Dias, University at Albany, Albany, NY, USA), employing Superdex G75 chromatography of fractions obtained after Sephacryl S-100 and immunoaffinity chromatography, as described in detail previously [24].The fractions possessing largely 24 kDa-FSHβ were pooled to generate FSH 24 and those possessing largely 18kDa-and 21 kDa-FSHβ were pooled to obtain FSH 18/21 .Recombinant human FSH produced in Chinese hamster ovary cells (batches AU012310 and BA024393), was a gift of Merck Serono (Mexico City, Mexico); according to the manufacturer [34], the batch-to-batch consistency for this recombinant FSH compound exhibits coefficients of variation that ranges from 7% to 15% for its most and least abundant isoforms.Equine FSH (batch VB-I-171) was purified from horse pituitaries obtained from Animal Technologies, Inc., (Tyler, TX, USA), as previously described [39] Rat granulosa cells culture Granulosa cells (GC) from immature (21 days old) Wistar rats pretreated for 4 days with 10 mg diethylstilbestrol (DES) through an implanted 10 mm x 1.5 mm silastic capsule, were collected and cultured following the method described by Jia and Hsue [40] with some modifications [41].Briefly, granulosa cells were added to 6-well (34 mm diameter) polystyrene culture dishes (Nunc, Roskilde, Denmark) at a density of 1.0 x 10

Data bioinformatics analysis
The main objective of the computational pipeline was to identify differences among the four FSH glycoforms tested at three distinct levels: a. gene expression; b. active biological processes; and c. perturbed pathways.To achieve this goal, we applied three different computational approaches: differential g e n e expression, over-representation analysis and pathway perturbation analysis.Data pre-processing Prior to data analysis, a series of pre-processing procedures were applied to the data in order to eliminate/diminish technical errors, biases and other sources of undesired variation Quality parameters of the FASTQ files were checked using the quality control FASTQC tool (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/).FASTQ reads were filtered and trimmed using the AfterQC tool [43]; default parameters were used to discard low-quality reads, trim adaptor sequences, and to eliminate poor-quality bases.The reads passing quality control parameters were aligned and quantified to the Rattus norvegicus reference transcriptome Rnor_6.0 using Salmon [44] version 0.8.2.Transcript-level quantifications were imported using the tximport package [45] version 1.22.0.
Differential expression Differential gene expression analysis was performed with the edgeR package [46], by comparing each FSH treatment at 0, 6, and 12 h.Since the number of replicates was limited, it was decided to conserve the top-200 differentially expressed genes (FDR < 0.01).Additionally, the similarities and differences between groups of genes reported as differentially expressed were compared under conditions: a.All FSH glycoform preparations at the same incubation time; and b.Different incubation times (0, 6, and 12 h) for the same FSH glycoform (eg.eqFSH, recFSH, FSH 18/21 , and FSH 24 ).Each of these comparisons followed distinct scientific questions.
On the one hand, the first comparison was useful to analyze the behavior of granulosa cells under the action of four distinctly different FSH glycoforms, on the other hand the latter comparison served to analyze the dynamics of a given FSH glycoform at different times, providing a set of snapshots of the temporal effects of the molecules.
Over-representation analysis By taking separately the overexpressed and underexpressed gene sets, over-representation analysis to detect whether the differentially expressed genes for each contrast were associated to a specific set of biological processes was performed.We implemented independent analyses for over-and underexpressed genes, based on the idea that sets of overexpresed genes may exacerbate a given process.Conversely, underexpressed genes may indicate a depleted/diminished behavior in a particular biological event.We used the Gene Ontology and the Kyoto Encyclopedia of Genes and Genomes (KEGG) [47] as databases for biological categories.Significance threshold was set at p<0.01 to consider a biological category as significant.
Pathway perturbation analysis Pathway perturbation was assessed using the gene expression data previously described.A subset of KEGG pathways was selected based on the information previously described [29].
Pathways were collected using the Graphite package [48].Contrasts were made between the experimental conditions described above and the baseline cell culture using the Generally Applicable Gene-set Enrichment for Pathway Analysis (GAGE) package [49].Multiple test correction was performed using the Benjamini-Hochberg post hoc method [50].Top 4 pathways by q-value (corrected p value) for each contrast were selected, and their network representation was merged into a single metapathway.Betweenness centrality (a measure for the relative importance of molecules for the communication in the corresponding pathways) for each gene in the metapathway was calculated.For each contrast, key genes were identified as those that had above-median logfold change in the differential expression analysis, and above-median betweenness centrality in the metapathway network.

Validation of FSH glycoforms-sensitive transcripts by real time (RT)-PCR
To validate differentially expressed genes in the RNAseq data, we analysed the mRNA level by RT-PCR of five selected FSH glycoform-sensitive genes (Pld1, Npy1R, Amh, Vegf-B, and Bcl2l1).cDNA synthesis was performed in 2.5 μg total RNA with 260/280 nm ratio >1.8 following the manufacturer'as instructions for the Maxima First Strand cDNA Synthesis Kit for RT-qPCR (Thermo Fisher), employing random primers and oligo dT.RT-PCR was performed employing 10 ng cDNA in a 20 μl reaction mixture containing 1 μl TaqMan Gene Expression Assay ® , 1μl TaqMan Gene Expression Assay for Actin or GAPDH reference gene, and 10 μl Master Mix TaqMan Universal PCR ® (all from Applied Biosystems, Waltham, MA, USA).All reactions were performed under the following conditions: 2 min at 50°C (UNG incubation), 10 min at 95°C (AmpliTaq Gold activation), followed by 40 cycles at 95°C for 15 sec (denaturation) and at 60°C for 1 min (anealing and amplification), employing a Step One plus thermocycler (Applied Biosystems).The results were normalized and expressed as fold changes in gene expression levels relative to control cDNA (2 -ΔΔCT method).The TaqMan gene expression assays employed were: Rn01493709_m1, Rn02769337_s1, Rn00563731_g1, Rn01454585_g1, Rn06267811_g1 (for the above described test genes, respectively) and Rn01426628_g1 and Rn01775763_g1 (for actin and GAPDH, respectively).Differences in expression levels of selected genes expressed in response to FSH stimulation as assesed by RT-PCR, were analyzed employing the Student's t test at an alpha level for significance <0.05.

General contrast between 6 and 12 hours
By comparing the four glycoforms at 6 and 12 h vs control (no FSH added) it was possible to determine how many genes were shared among these glycoforms.Figures 2A and 2B show the top200 differentially expressed genes for the four glycoforms at 6 h (A) and 12 h (B).The whole set of differentially expressed genes for each glycoform and contrasts can be found in the supplementary information (Tables S1 to S4).This analysis allowed to identify particular changes in gene expression induced by each glycoform.At 6 h, the number of genes shared among the glycoforms (center of Fig. 2A) were only 4 genes, whereas at 12 h this number increased to 118 genes (center in Fig. 2B).Thus, at 12 h the four glycoforms behaved more similarly to each other regarding the induction of differential gene expression.

Figure 2. Venn diagrams of the top-200 differentially expressed genes between FSH glycoforms at 6 h (A) and 12 h (B)
. Numbers inside the figure represent the number of genes shared in the corresponding set.Gray scale is proportional to the number of genes in the subset compared to the total of genes.Note that the number of genes shared between the four groups (center of the figure) at 6 h is only 4 genes; meanwhile, at 12 h, the very same set is 118 genes, meaning that at 12 h, the four compounds behave much similar than at 6 h.

Contrasts among each FSH glycoform vs control at 6 and 12 hours
Figure 3 A-D shows the contrast between 6 h and 12 h for each FSH glycoform.In all cases, the number of top200 differentially expressed genes induced at 6 h was higher than at 12h. recFSH showed the strongest difference between 6 h and 12 h since the intersections between these incubation times were only 34 genes (Fig. 3D).The results with all comparisons, as well as the statistics of each contrast can be found at https://osf.io/57j3f/.Figure 3. Contrast between 6 and 12 hours for each isoform of FSH.In the four cases, the top200 more differentially expressed genes are more different between 6 hours and control than 12 hours vs control.The contrast of recombinant FSH shows that there is a strong difference at 6 and 12 hours, since the intersection between 12 and 6 hours is only 34 genes.

Gene ontology enrichment
Each contrast was enriched by means of the DAVID online tool [51].Gene Ontology (GO) terms were obtained with the DAVID-standard p-values of enrichment (at p<0.01).This method was applied to explore those processes in which differentially expressed genes participate.The four different treatments were analyzed at 6 h and 12 h time points, using the list of overexpressed and underexpressed genes separately.
Processes associated with overexpressed genes after 6 hours of FSH stimulation (Tables 1   and 2, and Fig. 4) Exposure of granulosa cells to recFSH and FSH 18/21 was followed by induction of several enriched processes that included response to estradiol and calcium, angiogenesis (in the case of recFSH), and response to the gonadotropin stimulus and steroid biosynthetic process (for FSH 18/21 ) (Tables 1 and 2).For these two FSH compounds the only category shared was reponse to drugs.In contrast, FSH 18/21 shared several processes with eqFSH, including response to estrogen, ovulation, and cAMP response, the latter also shared with FSH 24 .This finding suggests that eqFSH and FSH 18/21 behaved more similarly in the early (6 h) response.
Interestingly, and in contrast to FSH 18/21 , FSH 24 did not exhibit any exclusively enriched process after stimulation during 6 h and showed only one shared biological process (cellular response to cAMP, shared with eqFSH and FSH 18/21 ). Figure 4 shows the enriched processes of overexpressed genes for the four glycoforms at 6 h, with the nodes corresponding to each glycoform colored in yellow.Processes associated with overexpressed genes at 12 hours of FSH stimulation (Tables 3   and 4, and Fig. 5).
Similar to the data shown by the Venn diagrams shown in Fig. 2, all compounds exhibited more shared processes at 12 h than at 6 h, as illustrated in Fig. 5 and Tables 3 and 4. In fact, eleven processes (including response to estrogen and cholesterol biosynthetic process), were shared by all FSH phenotypes.At this time, there were fewer treatment-specific processes than after 6 h of FSH exposure (Tables 1 and 2 and Fig. 5), being recFSH and FSH 24 (both tetraglycosylated FSH molecules) the glycoforms exhibiting a larger number of exclusive processes.
FSH 18/21 and eqFSH showed only one exclusive process at this time.
Figure 5. Biological processes associated with overexpressed genes at 12 hours after FSH glycoforms addition.Color code is the same as in Fig. 4. It is evident that the number of shared processes was larger at 12 h than at 6 h of exposure.See also legend of Fig. 4.   and 6, and Fig. 6).
For the underexpressed genes, there was a lower number of significantly enriched processes (n=12) despite the gene sets of overexpressed and underexpressed genes were of the same size (i.e.200).FSH 18/21 was the preparation exhibiting more negatively regulated exclusive processes (n=7), including positive regulation of angiogenesis, whereas for eqFSH and recFSH there was only one exclusively enriched process, regulation of cell growth and thyroid hormone metabolic process, respectively.Positive regulation of cell proliferation and cell migration were negatively associated processes between FSH 18/21 and eqFSH at this time.For FSH 24 we only identified cell response to cAMP as a negatively regulated shared process (with recFSH).Processes associated with underexpressed genes at 12 hours of FSH stimulation (Fig. 7 and Tables 7 and 8).
For underexpressed genes, a larger number of significantly enriched biological processes were identified at 12 h than at 6 h (Fig. 7).At this time, eqFSH was the preparation showing the largest number of significantly underregulated processes (both exclusive and associated) (n= 18 processes vs 3 processes at 6 h), from which 9 were shared and 10 were exclusive; the latter were more abundant than those detected with the exclusively overexpressed genes at the same time, during which only 1 exclusive process (cell adhesion) was detected (see Fig. 5).In contrast to the findings at 6 h, FSH  7 and 8).Timeline for each FSH phenotype (Figs.S1 to S4 in Supplementary information).FSH 18/21   As shown in Fig. S1, overexpressed gene-associated processes exhibited exclusive biological processes at each time (6 and 12 h), with a subset of shared processes at both times.
At 6 h, all FSH 18/21 -stimulated exclusive processes were associated with cell responses, whereas the four shared processes present at 12 h included responses to steroid hormones and drugs, as well as biosynthetic steroid process.At 12 h, 13 processes, including responses to estradiol, cholesterol, hypoxia, nutrient, and organic cyclic compound, as well as male gonad development appeared as significantly enriched.Meanwhile, among underexpressed genesrelated processes, positive regulation of proliferation and migration, wound healing, and regulation of gene expression appeared as exclusive at 6 h, and DNA replication, chromosome segregation, microtubule-based movement, and response to substances were enriched at 12 h.Shared underexpressed processes between 6 and 12 h were not evident.Contrary to those corresponding to overexpressed genes, there were fewer underexpressed genes-associated processes at 12h than at 6h. FSH 24   Interestingly, exposure to FSH 24 for 6 h yielded the same associated process based on overexpressed and underexpressed genes, i.e. cellular response to cAMP (Fig. S2 bottom).
This was the only process enriched at 6 h for this particular compound.On the other hand, several processes associated with cellular response to different stimuli (in the case of overexpressed genes), and cell division, DNA replication, and cell movement (for underexpressed genes), among others, were observed after 12 h of exposure to this glycoform.
eqFSH For overexpressed genes-related processes stimulated after 6 h exposure to eqFSH, the response to cAMP and processes associated with response to FSH stimulus (e.g.ovulation, response to peptide hormone) were exclusive, whereas response to estrogen and cAMP were shared with those processes stimulated after exposure during 12 h (Fig. S3).Among the exclusive processes detected after 12 h of exposure were cholesterol metabolic and biosynthetic processes, extracellular matrix organization, and cell adhesion.For the case of underexpressed genes-related processes, regulation of cell growth was the only exclusive function, and cell migration and proliferation were shared with processes stimulated at 12 h by this glycoform.Biological processes associated with cell proliferation were among the 12 h exclusive ones.Interestingly enough was finding that the response to estradiol was an enriched process associated with overexpressed and underexpressed genes at 12 h, similar to the cellular response to cAMP observed for FSH 24 at 6 h.In both cases, stimulation followed by inhibition of these processes may be related to their selective (biased) desensitization following stimulation with FSH 24 and eqFSH.recFSH At 6 h, microtubule depolymerization, mitotic sister segregation, and angiogenesis were among the exclusive, overexpressed gene-related processes (Fig. S4).Elements shared between 6 h and 12 h were response to drug, hypoxia and estradiol, as well as extracellular matrix organization.Comparatively, more exclusive processes were detected after 12 h of recFSH stimulation, including female gonad development, steroid biosynthetic process, and response to estrogen, among others.
The only process significantly enriched from underexpressed genes at 6 h was thyroid hormone metabolic process (Fig. S4, bottom left).Microtubule-based movement, cell proliferation and intracellular signaling appeared as processes enriched after 12 h of recFSH exposure.The case of chromosome segregation was interesting in that this process was significantly enriched for overexpressed genes at 6 h and at 12 h for the underexpressed ones, perhaps pointing to a transitional process.Another interesting finding was that cellular response to cAMP appeared as an enriched process related to underexpressed genes at 6 h but also to overexpressed genes at 12 h.This could mean that this particular process changed its behavior from one pole to the other during a 12 h exposure to recFSH.

Pathway perturbation analysis
Pathway perturbation analyses allow linking the global effect of gene perturbation observed in a given set of experiments with known functional and molecular processes.Based on previous knowledge of the effects of FSH on FSHR-mediated intracellular signaling [29,30], we focused our analysis on cAMP-PKA, MAPK-, and PI3K/AKT-mediated signaling.Our analyses identified these pathways as altered in virtually all experimental comparisons.
Behavior of the cAMP signaling pathway in response to each FSH compound (Suppl.

Information Figs. S5 and S6).
A closer inspection on the cAMP pathway, allowed us to identify that said dysregulation is driven by different sets of perturbed genes.Figures S5 and S6 show the differentially expressed genes in KEGG-resolved cAMP signaling pathway for each treatment at 6 and 12 h.As shown, each compound exhibited a time-dependent individual molecular signature.For example, 6 h after FSH 18/21 and FSH 24 exposure, membrane-bound FSH-stimulated FSHR expression was found generally underexpressed, most probably due to internalization and subsequent degradation of the receptor.Simultaneously, phosphodiesterase was overexpressed to prevent prolonged activation of the cAMP/PKA pathway.In the case of FSH 18/21 , CREB was overexpressed at 6 h, whereas no changes in this transcription factor were observed in the case of FSH 24 at this time.After 6 h of FSH 24 exposure, G protein-coupled receptors (GPCRs) that responded to colinergic and GABAergic stimuli were modestly overexpressed whereas no change was observed for FSH 18/21 at this time.Changes similar to those provoked by FSH 18/21   exposure for 6 h were observed for eqFSH, whereas for recFSH no GPCR exhibited dysregulation.
After 12 h of FSH exposure, the effects of gene expression on pathway molecules were more consistent.Most molecules in the cAMP-PKA pathway (eg.CREB) exhibited underexpression at this time point, regardless of their cellular or pathway position.The exception was AMPA (αamino-3-hydroxy-5-methyl-4-isoxazolepropionic acid) receptor, which is stimulated by glutamic acid present in follicular fluid [52] and in turn regulates signaling molecules (presumably via calcium influx) involved in ovulation [53]; AMPAR showed overexpression with all treatments.
At this time, signaling was also very similar between FSH 18/21 and eqFSH with the exception of GPCRs for glycoprotein hormones and c-Jun N-terminal kinase (JNK; a kinase involved in cell cycle progression, mitosis and apoptosis [54,55]) signaling overexpression stimulated by eqFSH.It is interesting to note the contrast between recFSH-and FSH 24 -(both fullyglycosylated human FSH compounds) stimulated signaling.JNK-mediated signaling and CREB-regulated gene expression (at 6 h) as well as MEK and JNK signaling (at 12 h), were markedly overexpressed after exposure to recFSH but not to FSH 24 .In contrast, after exposure to FSH 24 both CREB-regulated gene expression and JNK signaling remained unmodified or were underexpressed at 6 h and 12 h, respectively, whereas MAPK signaling was slightly overexpressed at both times.

Behavior of the MAPK and PI3K-AKT signaling pathways in response each FSH compound (Suppl. Information Figs. S7 to S10).
The effects of the FSH compounds tested on the MAPK signaling pathway also differed among the glycoforms (Figs.S7 and S8).Notably, MEK1 and MEK2 pathways were upregulated by FSH 18/21 at 6 h but not by FSH 24 , whereas at 12 h, ERK signaling was modestly upregulated in response to both glycoforms.Meanwhile, in the case of eqFSH, MEK1 and 2 as well as ERKsignaling were modestly overexpressed at 6 h and 12 h, respectively, whereas after exposure to recFSH, both MEK2 and ERK were overexpressed at 12 h but not at 6 h.Worthy of note in the PI3-AKT signaling pathway (Figs.S9 and S10) was the overexpression of SGK (serum and glucocorticoid-regulated kinase, which is involved in survival signaling, growth, and proliferation [56]) after 12 h of cell exposure to FSH 18/21 , recFSH, and eqFSH, whereas this signaling pathway remained unaltered or markedly underexpressed at 6 h for all glycoforms.

Validation of FSH glycoforms-regulated gene expression by RT-PCR
To validate the results obtained from the RNA-seq experiment, we analyzed by RT-PCR the relative mRNA expression of four over-expressed and one under-expressed genes under the stimulation of distinct glycoforms (Fig. 8).Genes examined that were putatively up-regulated by FSH glycoforms included: a) Phospholipase D1 (Pld1), which encodes an enzyme implicated in a number of cellular pathways, including signal transduction mediated by GPCRs and receptor tyrosine kinases, subcellular trafficking, and regulation of mitosis [57]; b) Neuropeptide Y (NPY) receptor Y1 (Npy1R), which encodes for a GPCR involved in Npy-regulated granulosa cell proliferation and apoptosis [58,59]; c) Antimüllerian hormone gene (Amh), encoding a glycoprotein belonging to the transforming growth factor beta superfamily, which is involved in follicular development [60]; and d) Vascular endothelial growth factor B (Vegf-B), a gene that encodes the VEGF-B protein, a factor involved in angiogenesis [61].The FSH down-regulated gene examined was BCL2 like 1 (Bcl2l1) which acts as anti-or pro-apoptotic regulators that are involved in a wide variety of cellular activities [62], including prolactin signaling [63].The selection of those genes for validation was based on the role they play in FSH-stimulated follicular maturation, rather than on fold change variation, and their corresponding mRNAs were quantified in three independent biological replicates (samples coming from three different culture wells) from granulosa cells exposed to different FSH glycoforms during 6h or 12h.The mRNA expression of the 4 up-regulated genes in response to particular glycoforms vs zero hours [with statistically significant (p<0.01 or p<0.05) differences] was the following: Amh and Npy1R for recFSH and FSH 18/21 at 6 h, respectively; Pld1 for FSH 18/21 and eqFSH at 12 h; and Vegf-B for recFSH at 12 h.For the down-regulated gene (Bcl2l1) a statistically significant difference (p<0.05) was found for FSH 24 at 6 h.These results indicated the existence of a positive validation (Fig. 8).The use of only shown that hypoglycoslated FSH exhibits higher FSHR binding activity and receptor affinity as well as potency and efficiency to trigger intracellular signaling than fully glycosylated FSH [12,14,15,24,25,65,66]) despite its relatively shorter plasma half-life [66,67].In fact, when injected to Fshb null female mice, both hypo-and fully glycosylated FSH stimulated ovarian weight response and induced particular ovarian genes in a similar fashion, whereas in Fshb null male mice, hypoglycosylated FSH 18/21 was more active than fully-glycosylated FSH 24 in inducing FSH-responsive genes and Sertoli cell proliferation [25].More recently, short-term in vivo studies in prepubertal mice showed that FSH 18/21 was more effective than FSH 24 in promoting follicle development and health, as well as in rescuing granulosa cells from apoptosis [66].Although ovarian gene transcription in vivo was stimulated by both glycoforms, FSH 18/21   induced greater activation of Gs-mediated cAMP-PKA signaling as well as PI3/K and MAPK/ERK signaling pathways, with greater expression of early response genes upon administration of this particular glycoform [66].These data were concordant with findings from previous studies showing that in vivo administration of more basic, short-lived FSH isoforms to hypophysectomized rats, mantained equally or even more efficiently granulosa cell proliferation than their more acidic counterparts [20].
In the present study, we analyzed whether FSH glycoforms exhibiting distinct glycosylation patterns differentially induce gene transcription and activation of biological processes and signaling pathways in cultured rat granulosa cells exposed during 6h and 12h to equal doses of the glycoforms.Analysis of differentially expressed genes identified particular changes in gene expression induced by each glycoform, with a limited number of genes shared among the glycoforms at 6 h, followed by a sharp increase at 12 h, time at which the four glycoforms behaved more similarly for inducing differential gene expression.These data reflect the relatively early differential effects of the glycoforms on gene expression, with a platau reached at 12 hours.Analysis of individual FSH glycoforms revealed the expression of more FSHdependent genes at 6 h than at 12 hours, particularly for the hypo-glycosylated preparations (FSH 18/21 and eqFSH), which was expected given that FSH-mediated signaling occurs rapidly after receptor activation in in vitro conditions and that hypoglycosylated FSH exhibits a more favorable FSHR binding profile than its fully-glycosylated counterpart from kinetic and thermodynamic points of view, as disclosed by molecular dynamics simulations [12,26].
Nevertheless, the data also showed that even after 12 h of FSH exposure the effects of FSH stimulation on gene expression persisted, independenly of the kinetics by which each glycosylation variant activates the FSHR and stimulates intracellular signaling.
Gene ontology analysis revealed that exposure of granulosa cells to each FSH compound during 6 h or 12 h was differentially associated with induction or inhibition of one or several genes linked to distinct biological processes.Some of these processes were exclusive or unique for a given FSH compound whereas others where shared by or associated with several FSH glycoforms, demonstrating the distinctly different ability of each glycoform to activate/inhibit FSH-dependent biological processes.In this analysis, it was interesting to find that in contrast to hypo-glycosylated and recFSH, stimulation with FSH 24 and eqFSH for 6 h, was poorly associated with exclusive overexpressed processes.In fact, in the case of FSH 24 no exclusive biological processes and only one associated process (cAMP response) were clearly identified at this time.This observation contrasted with the findings after 12 hours of FSH exposure, where multiple associated processes where shared among all glycoforms, again emphasizing the long-term effects of FSH on biological effects in vitro.Here it is important to emphasize that at this time (12 h) the exclusive processes stimulated by FSH 18/21 where virtually absent, underlying the short and acute effects of this particular glycoform on granulosa cell responsiveness, an observation that is in agreement with a previous study on the effects of FSH 18/21 and FSH 24 on ovarian global trascriptomics in vivo [66].In this scenario it is possible to conclude that recFSH (tetra-glycosylated), was the most potent preparation to evoke particular/exclusive biological processes at 6 h and 12 h and that the effects of all FSH compounds converged in a number of similar processes in a time-dependent fashion, as it probably occurs in vivo during follicular maturation.Meanwhile, individual differences in time also were observed in processes associated to underegulated genes, with exposure to FSH 18/21 for 6 h resulting in more underegulated genes thatin those by other glycoforms but not later, at 12 h, time at which FSH 18/21 did not show any exclusive process, an observation that was in sharp contrast with the effects of the remaining glycoforms.These data indicates that different glycoforms also have distinct abilities to turn off selective genes in a time dependent manner, with the FSH 18/21 hypo-glycosylated glycoform acting faster than the other compounds.This observation is in line with the ability of this particular glycoform to efficiently trigger intracellular signaling and also to rapidly evoke down-regulation of stimulated signaling and also probaby resensitization of the FSH 18/21 /FSHR complex [12,14,24].In fact, it has been shown that FSH glycosylation can influence the turnon time of FSHR-mediated signaling [24,68,69].
Pathway perturbation analysis unveiled interesting differential effects of the FSH glycoforms on intracellular FSH-regulated signaling.We analyzed by this bioinformatics procedure three well-known pathways stimulated by gonadotropins, the canonical cAMP-PKA as well as the MAPK and PI3K-AKT signaling pathways [29,30,70].It is believed that these signaling cascades lead to fine-tuning regulation of the FSH stimulus, where activation and/or inhibition of their downstream activated components vary depending on the cell context, cell developmental stage, and concentration of the ligand and the receptor.Simultaneous activation of these signaling modules eventually converges on the ligand/receptor-integrated biological responses, which include cell proliferation, differentiation and survival of responsive cells and, at the molecular level, differential gene expression, as observed in the present study.In the cAMP pathway, it was interesting to find that for FSH 18/21 but not for the fully glycosylated FSH 24 , the transcription factor CREB was overactivated at 6 hours indicating a rapid activation of the cAMP-PKA pathway by this particular naturally occurring glycoform, whereas after 12 h stimulation of this pathway persisted active for both glycoforms, albeit to a lesser extent after FSH 24 exposure.The overall data underline the differential activation of the cAMP signaling pathway by these two glycoforms.
The MAPK and PI3K-AKT cascades are also part of the intertwined signals regulated by the FSH/FSHR complex and both are involved in controlling follicle development as well as in regulating gene transcription [29,[71][72][73].In these pathways, the effects of the FSH compounds tested on the MAPK signaling pathway also exhibited time-related differences.For example, MEK1 and MEK2 pathways were upregulated by FSH 18/21 at 6 h but not by FSH 24 , whereas at 12 h, ERK signaling was modestly upregulated in response to both glycoforms.Worthy of note is the fact that activation of FSHR involves cross-talk with pathways regulated by growth factors and the estrogen GPCR (GPER) [56,74,75] ).In the former, RTK plays an important role in folliculogenesis and ovulation [76,77], and it is known that this kinase may activate JNK in a Src and PI3K dependent fashion, whereas in the latter FSHR/GPER heterodimers trigger antiapoptotic/proliferative signaling through the Gβγ dimer [74,78].In this vein, it was interesting to find that all FSH glycoforms activated RTK and MEK after 6 h of FSH exposure, but not later, underlying the important role of this pathways on FSH-evoked cell proliferation and angiogenesis.
Regarding all these signaling pathway analyses, it is important to emphasize on that pathway perturbation captures a static representation of what is essentially a dynamic process involving several signaling cascades regulated by distinct factors either directly by the FSHR or via crosstalk with other membrane receptors.As such, unexpected behaviors such as under-and overexpression of certain molecules can be capturing dynamic feedback effects.With this in mind, the differences and similarities observed between different pathways can be thought to be capturing different kinetic behaviors by the distinct FSH treatments.In any case, these data underline the distinct time-and glycosylation-dependent fingerprints for each particular FSH compound.
In summary, the present transcriptomic analysis allowed dissection of some distinctly different glycosylation-dependent effects of the human FSH glycoforms after exposure of cultured rat granulosa to these compounds at different times, comparing these effects with those of eqFSH, which albeit differing in amino acid sequence and glycosylation pattern, exerts potent effects at the human FSHR [24].Of particular interest are the data on the transcriptomic effects of the naturally-occurring FSH glycoforms, which are consistent with previous in vivo and in vitro studies showing that hypoglycosylated FSH 18/21 exhibits greater biological activity at the target cell level and that both FSH 21/18 and FSH 24 initiate different FSHR-mediated intracellular signaling activation, eventually leading to differential impacts on follicle development.[66,79].
The overall findings are important to better understand the physiological significance of FSH glycoforms in follicular maturation and ovulation in naturally-occurring cycles in women.

Figure 1 .
Figure 1.Typical glycans attached to human pituitary FSH (hFSH), human recombinant FSH produced by Chinese hamster ovary cells (recFSH) glycoforms, and equine FSH (eFSH).The green and cyan bars and ribbons (3-D structures at right) indicate the common-α and hormonespecific FSHβ subunits, respectively.N-glycosylation sites are indicated by the numbers below the bars.The eFSHα subunit has 4 additional residues at the N-terminus, accounting for the difference in numbering.The glycan at position β7 in eqFSH (black dotted square) is absent in the bulk (90%) of the molecules contained in highly purified preparations[27].Glycans in recFSH were taken from Mastrangeli et al[28].Human FSH glycoform models created with the GLYCAM web tool are shown on the right with the subunits rendered as cartoons using PyMol.The same color scheme is employed for subunits and glycans, which are rendered as spheres.Partially obscured Asn24 glycan is colored dark blue to distinquish it from Asn 7 glycan (light blue).

Figure
Figure 4. Enriched processes of overexpressed genes for the four FSH glycoforms at 6 h, with the nodes corresponding to glycoforms in yellow.In this network representation, red squares represent significantly (p<0.01)enriched biological processes in the treatment in which there is a link.There are some shared processes i.e. biological functions that appear enriched in more than one phenotype, which is indicated by more than one link.It is important to consider the number of associated processes for each isoform: eqFSH, 8 processes (2 exclusive); recFSH, 11 processes (only 1 exclusive); FSH 18/21 13 processes, from which 6 are exclusive, and FSH24  only one associated process.In this figure, the intersection between FSH 18/21 at 6 h is depicted against the other 3 isoforms at the same time.Table1.Processes associated with overexpressed genes after exposure for 6 hours to FSH 18/21 and FSH24 .

Figure 6 .
Figure 6.Biological processes associated with underexpressed genes at 6 hours of exposure to each of the four FSH treatments with the nodes representing each glycoform in yellow.In this representation, blue color represents the enriched processes.

Figure 7 .
Figure 7. Biological processes associated with underexpressed genes at 12 hours of exposure to each of the four treatments.In this representation, blue color represents the enriched processes.
Invitrogen, Waltham, MA, USA), and 0.1% bovine serum albumin (Sigma), and cultured for 24 h at 37ºC in a humidified atmosphere of 95% air-5% CO 2 .At the end of the culture period, cells were washed, redissolved in insulin-and serum-free, 6viable cells/well in serum and insulin-free 2 ml McCoy's 5a medium (Life Technologies, Grand Island, NY, USA), pH 7.0, supplemented with 2 nM glutamine (Sigma-Aldrich Inc., St. MO, USA), antibiotic (penicillin plus streptomycin) reagent (GC to increasing doses of each FSH preparation, in which maximal estrogen production was achieved with the 100 ng/ml dose (not shown).RNA isolation and sequencingTotal RNA was isolated from individual GC wells using the TRIzol TM reagent (Thermo Fisher Scientific, Waltham, MA USA)) and the Direct-zol RNA kit (Zymo Research, Irving, Ca,

Table 1 .
4. Enriched processes of overexpressed genes for the four FSH glycoforms at 6 h, with the nodes corresponding to glycoforms in yellow.In this network representation, red squares represent significantly (p<0.01)enrichedbiologicalprocesses in the treatment in which there is a link.There are some shared processes i.e. biological functions that appear enriched in more than one phenotype, which is indicated by more than one link.It is important to consider the number of associated processes for each isoform: eqFSH, 8 processes (2 exclusive); recFSH, 11 processes (only 1 exclusive); FSH 18/21 13 processes, from which 6 are exclusive, and FSH24only one associated process.In this figure, the intersection between FSH 18/21 at 6 h is depicted against the other 3 isoforms at the same time.Processes associated with overexpressed genes after exposure for 6 hours to FSH 18/21 and FSH24.

Table 2 .
Processes associated with overexpressed genes after 6 hours of cell exposure to eqFSH and recFSH.

Table 4 .
Processes associated with overexpressed genes after 12 hours exposure to eqFSH and recFSH.

Table 6 .
Processes associated with underexpressed genes after 6 hours exposure to eqFSH and recFSH exposure.

Table 8 .
Processes associated to underexpressed genes after 12 hours exposure to eqFSH and recFSH.
3 replicates was justified by the difficulties in generating more material from such experiments.It is clear that using more replicates may better support the PCR vs RNAseq results where statistical power is increased by the measurement of thousand of targets at the same time.DISCUSSION Previous studies have demonstrated that hypoglycosylated FSH (FSH 18/21 ) is the predominant form secreted by women in reproductive age and that a shift to fully-glycosylated glycoforms occurs as age advances [15, 64, 65/217].In vitro and in vivo studies also have