Transcriptome profile of goat folliculogenesis reveals the interaction of oocyte and granulosa cell in correlation with different fertility population

To understand the molecular and genetic mechanisms related to the litter size in one species of two different populations (high litter size and low litter size), we performed RNA-seq for the oocytes and granulosa cells (GCs) at different developmental stages of follicle, and identified the interaction of genes from both sides of follicle (oocyte and GCs) and the ligand-receptor pairs from these two sides. Our data were very comprehensive to uncover the difference between these two populations regarding the folliculogenesis. First, we identified a set of potential genes in oocyte and GCs as the marker genes which can be used to determine the goat fertility capability and ovarian reserve ability. The data showed that GRHPR, GPR84, CYB5A and ERAL1 were highly expressed in oocyte while JUNB, SCN2A, MEGE8, ZEB2, EGR1and PRRC2A were highly expressed in GCs. We found more functional genes were expressed in oocytes and GCs in high fertility group (HL) than that in low fertility group (LL). We uncovered that ligand-receptor pairs in Notch signaling pathway and transforming growth factor-β (TGF-β) superfamily pathways played important roles in goat folliculogenesis for the different fertility population. Moreover, we discovered that the correlations of the gene expression in oocytes and GCs at different stages in the two populations HL and LL were different, too. All the data reflected the gene expression landscape in oocytes and GCs which was correlated well with the fertility capability.


Scientific Reports
| (2021) 11:15698 | https://doi.org/10.1038/s41598-021-95215-z www.nature.com/scientificreports/ Library preparation, sequencing and data analysis. Library preparation and sequencing. We prepared libraries for the oocytes and GCs separately. For GCs, total RNA was extracted as described in our early articles 30 . Then Ribo-Zero rRNA Removal Kit was used to remove rRNA to purify the RNA samples. Then Ribo-Zero rRNA Removal Kit and DNA polymerase I and RNaseH were used to make the double strands cDNA followed by the sequencing by Illumina HiSeq X Ten equipment following the manufacturer's instruction. For the oocytes, total RNA was isolated as described in our early article 30 . Then the mRNA was enriched by the Oligo (dT) beads. Then the mRNA samples were used to make the double strand cDNA templates for the sequencing by Illumina HiSeq X Ten equipment following the manufacturer's instruction.
Gene expression levels and principal component analysis. The libraries were aligned against the Capra_hircus genome 31 . Only reads showing a unique match to the genome and less than five mismatches were further filtered to eliminate duplicates with picard. The bam files containing non-duplicated reads were employed as input for Cufflinks (v.2.2.1), together with Ensembl gene annotation, to obtain FPKM (Fragments Per Kilobase of exon model per Million mapped fragments) data. Genes were subjected to analytical procedures if FPKM > 0.5. All statistical analyses were conducted in R software (Biase FH, Kimble 2018). To identify the main sources of variation in the dataset, we employed the FPKM values as the input for principal component analysis using the FactorMiner R package. The significance of the principal components was obtained with the Seurat package via a permutation test, after 1000 randomized samplings 24 .
Identification of putative ligands and receptors in cumulus-oocyte complexes. First, we produced a comprehensive protein-protein interaction (PPI) database as described in the articles 24,[32][33][34] . Briefly, the gene from our data were used to search the ligand-receptor pairs as reported in these articles. Then the ligand-receptor pairs from GCs or oocytes were mapped the gene identifier 24 .

Calculation of the pairwise correlation of transcript levels between genes expressed in oocytes and surrounding
GCs. We calculated the Spearman's correlation coefficient 30 between the expression levels of genes expressed in oocytes and the surrounding GCs. Then the data were expressed by heatmap package in R stadio.
Functional annotation. The function of the different expressed genes was enriched by Metascape online.
q-RT-PCR. The procedure for mRNA q-RT-PCR was reported in our early publications using qPCR Master Mix (Roche, German) by the Roche LightCycler1 480 (Roche, German) 31 . The primers for qPCR analysis were synthesized by Invitrogen and present in Table S1. Briefly, total RNA was extracted using TRIzol Reagent (Invitrogen Corp., Carlsbad, CA, USA) and purified using an RT2 qPCR-Grade RNA Isolation Kit from SABiosciences Co., Ltd (MD, USA). Total RNA was quantified using a Nanodrop 3300 (ThermoScientific, DE, USA). The quality of RNA was controlled by the A260:A280 ratio being greater than 2.0 and confirmed by electrophoreses, with a fraction of each total RNA sample with sharp 18S and 28S ribosomal RNA (rRNA) bands as reported in our recent publication. One microgram of total RNA was used to make the first strand cDNA in 20 μl. The program for the reaction of miRNA and lncRNA was 25

Results
The transcript profiles of oocytes and granulosa cells during folliculogenesis. In order to discover the interactions of GCs and oocytes in the folliculogenesis, we analyzed the gene expression in different stages of oocytes and GCs of different fertility populations (HL, LL). We found that during the development, the gene expression profiles changed dramatically in the different stages of oocytes (Fig. 1a, Tables S2, S3). The gene expression profiles for different stages of GCs did not change much (Fig. 1a). In order to discover the molecular difference in oocytes and GCs, we used Seurat program to detect the marker genes and found the highly changed genes in oocytes and GCs, such as HSD17B2 and INHA (Fig. 1b). Moreover, the PCA analysis discovered that important marker genes in oocytes and GCs in PC1 component (Fig. 1c). The cluster analysis www.nature.com/scientificreports/ by UMAP program found that GCs and oocytes can be separated very well (Fig. 1d). However, there were two subclusters of oocytes ( Fig. 1d) which indicated that there was some difference for gene expression in different stages of oocytes. Based on the significant expressed marker genes, we found lots of dramatically different expressed genes in oocytes or GCs (Fig. 2), such as GDF9 and BMP15 which were highly expressed in oocytes while FST, FSHR, and LHCGR which were highly expressed in GCs (Fig. 2a,b). At the same time, we discovered new marker genes which can be used to separate the different cell types, such as GRHPR, GPR84, CYB5A, ERAL1 for oocytes ( Fig. 2c), while JUNB, SCN2A, MEGE8, ZEB2, EGR1, and PRRC2A for GCs (Fig. 2d). These genes can be potential marker genes of goat follicular cells. At the same time, the protein of the different genes for oocytes (GPD1) and GCs (TST) were analyzed by IHF staining (Fig. 2e,f).

Dynamic gene expression profiles in different stages of oocytes of the two fertility populations.
There were two population of goats used in this investigation: high fertility group (> 3/litter; HL) and low fertility group (≤ 2/litter; LL). The main purpose of this investigation was to explore the difference in the gene expression in different stages (large follicle > 7 mm in size; medium follicle 3-7 mm in size; small follicle < 3 mm in size) of oocytes and GCs in the two fertility population goats.
PCA analysis found that the gene expression pattern was closer for the oocytes in large follicles in high fertility group (HLO) than that in low fertility group (LLO) (Fig. 3a). Similarly, the gene expression pattern was closer for the oocytes in medium follicles in high fertility group (HMO) than that in low fertility group (LMO) (Fig. 3b). Moreover, the gene expression pattern was closer for the oocytes in small follicles in high fertility group (HSO) than that in low fertility group (LSO) (Fig. 3c). The data suggested that the oocytes in high fertility group were more homogeneous than that in low fertility group.
There were 161 genes highly expressed in HLO while 68 genes were highly expressed in LLO (Fig. 3d). The functions of the changed genes were enriched by Metascape online. Most of the enriched functional pathways were related to metabolism which was related to the growth of large follicle to prepare for ovulation (Fig. 3d).
There were 707 genes highly expressed in HMO while 119 genes were highly expressed in LMO (Fig. 3e). Most of the enriched functional pathways were related to transport which was related to the interaction of oocytes and GCs for the exchange of materials (Fig. 3e).
There were 218 genes highly expressed in HSO while 337 genes were highly expressed in LSO (Fig. 3f). The enriched functional pathways were in wide range of molecular functions. There were some pathways related to cell defense in the increased genes (Fig. 3f). The data in this section suggested that for the gene expression, the big difference was in the oocytes of the medium follicles. And the HMO has the highly developmental potential to HLO. Dynamic gene expression profiles in different stages of GCs of the two fertility populations. PCA analysis uncovered that the gene expression pattern was similar for the GCs in large follicles of high fertility group (HLG) and low fertility group (LLG) (Fig. 4a); in medium follicles of high fertility group (HMG) and low fertility group (LMG) (Fig. 4b); in small follicles of high fertility group (HSG) and low fertility group (LSG) (Fig. 4c). www.nature.com/scientificreports/ There were 5664 genes highly expressed in HLG while 2944 genes were highly expressed in LLG (Fig. 4d). There were two pathways related to differentiation and embryo development were enriched in the increased genes in HLG (Fig. 4d).
There were 4323 genes were highly expressed in HMG while 5907 genes were highly expressed in LMG (Fig. 4e). There were some pathways related to cell division, and embryonic development enriched in the increased genes in HMG (Fig. 4e) which further indicated that GCs in the large follicles of high fertility group pose more potential for the oocyte development.
There were 3925 genes were highly expressed in HSG while 5992 genes were highly expressed in LSG (Fig. 4f). The enriched functional pathways were diverse in this comparation (Fig. 4f) which suggested that at small follicle stage, the GCs in high fertility and low fertility groups are close.
Some of the gene expression levels were confirmed by q-RT-PCR (Fig. 4g).  NOTCH signaling pathway and TGF-β signaling pathway are the two most important signaling pathways in the folliculogenesis, and we found that many interesting unique protein-protein interactions (PPIs) of ligands and receptors for the oocytes and GCs. In NOTCH signaling pathway, the ligands DLL3 and JAG1 were specifically highly expressed in oocytes, while the ligand JAG2 was highly expressed in early stages of GCs and all stages of oocytes (Fig. 5a). The receptors NOTCH1, NOTCH2 were highly expressed in GCs while NOTCH3 was expressed in oocytes and GCs. Especially, NOTCH1 and NOTCH2 were very highly expressed in early stages of GCs (Fig. 5a), while NOTCH3 was very highly expressed in late stages of GCs (Fig. 5a). And the target gene HES1 was specifically expressed in oocytes (Fig. 5a).
In TGF-β signaling pathway, the ligands GDF9 and BMP15 were specifically expressed in oocytes, while their receptors BMPR1A, BMPR1B and BMPR2 were expressed in oocytes and GCs while the expression was time dependent (developmental stages) (Fig. 5b). The target genes ID3 was expressed in oocytes and ID4 was expressed in GCs and oocytes (Fig. 5b). The junction protein GJC1 and GJA4 were highly expressed in oocytes while GJA1 and GJA9 were highly expressed in GCs (Fig. S1).
Moreover, based on the gene expression the ligands and receptors pairs have been identified in different stages of oocytes and GCs of the two population goats (Figs. 6, 7, 8). 44 ligand-receptor pairs have been identified for the ligands in GCs and receptors in oocytes, while 65 ligand-receptor pairs have been found for the ligands in oocytes and receptors in GCs in the large follicles of high fertility group (HL) (Fig. 6). 27 ligand-receptor pairs have been identified for the ligands in GCs and receptors in oocytes, while 42 ligand-receptor pairs have been found for the ligands in oocytes and receptors in GCs in the large follicles of low fertility group (LL) (Fig. 6). 47 ligand-receptor pairs have been identified for the ligands in GCs and receptors in oocytes, while 68 ligandreceptor pairs have been found for the ligands in oocytes and receptors in GCs in the medium follicles of high fertility group (HM) (Fig. 7). 50 ligand-receptor pairs have been identified for the ligands in GCs and receptors in oocytes, while 97 ligand-receptor pairs have been found for the ligands in oocytes and receptors in GCs in the medium follicles of low fertility group (LM) (Fig. 7). 50 ligand-receptor pairs have been identified for the ligands in GCs and receptors in oocytes, while 69 ligand-receptor pairs have been found for the ligands in oocytes and receptors in GCs in the small follicles of high fertility group (HS) (Fig. 8). 55 ligand-receptor pairs have been identified for the ligands in GCs and receptors in oocytes, while 82 ligand-receptor pairs have been found for the ligands in oocytes and receptors in GCs in the small follicles of low fertility group (LS) (Fig. 8).
The functions of these PPI have been enriched by Metascape online. There were some specific functional pathways in the different stages of follicles of high or low fertility population. The interesting pathways "Signaling pathways regulating pluripotency of stem cells" and "TGF-beta signaling pathway" were enriched in the HL group not in LL population (Fig. S2). The pathways "developmental growth" and "epithelial cell differentiatio" were enriched in the HM group not in LM population (Fig. S3). The pathways "developmental growth", "MAPK cascade", "ECM-receptor interaction" and "GPCR ligand binding" were enriched in the HS group not in LS population (Fig. S4). The data suggested that the ligand-receptor interactions were different in the two fertility population which may be related to the fertility capability in goat.

Co-regulated gene expression between the oocyte and GCs of the two fertility populations.
The analysis of the gene co-expression between oocytes and GCs in different stages of follicles in two different fertility populations was performed in current investigation to explore the difference in the interactions of oocytes and GCs at different stages and in different populations.
The patterns of correlation of gene expression for different stage of follicles and different populations were present in Fig. 9. The X-axis was for oocytes and the Y-axis was for GCs. The very positively expressed genes in oocytes and GCs at different stages in the two populations were enriched to identify the functions. It was interesting that similar functional pathways were identified for oocytes and GCs with the strong positive correlation in all the stages of follicle in both high and low fertility populations. These functional pathways were related to  www.nature.com/scientificreports/ RNA processing, protein modification (Fig. 9). And the cell cycle related pathways were enriched in the medium and small follicles not large follicles which indicated that in the early stage oocytes and GCs cell proliferation, division, or differentiation (Fig. 9b,c,e,f). There were some pathways enriched only in high fertility population, not in low fertility population: such as "microtubule-based process", "positive regulation of telomere maintenance via telomerase" (Fig. 9a,d); "chromosome segregation", "cytoplasmic sequestering of protein" (Fig. 9b,e); "posttranscriptional regulation of gene expression", "cell cycle checkpoints" (Fig. 9c,f). The data suggested that the oocytes and GCs were different with the correlation in gene expression for the two populations.

Discussion
As the technology development, RNA-seq is a very useful tool to uncover the transcriptomic changes of oocytes and granulosa cell and to discover the interaction of different types of cells and even the mechanisms of folliculogenesis 23 . In current investigation, we explored the gene expression landscape of oocyte and GCs during different developmental stages by the tube based single cell RNA-seq. We aimed to discover the gene expression patterns in these two types of cells oocytes and GCs, the dynamic gene expression changes during development, and the interactions of genes in both sides to deeply understand the mechanism of folliculogenesis in two different populations. The data will be used to help the regulation of goat fertility (litter size).
First, we identified a set of potential genes in oocyte and GCs as the marker genes which can be used to determine the goat fertility capability and ovarian reserve ability. The data showed that GRHPR, GPR84, CYB5A and ERAL1 were highly expressed in oocyte while JUNB, SCN2A, MEGE8, ZEB2, EGR1and PRRC2A were highly expressed in GCs. It has been reported that ZEB2 has been identified as a new marker gene in human follicular GCs during the folliculogenesis 23 . Moreover, estrogen is mainly produced by aromatase in GCs. Jun, highly expressed in GCs, regulated aromatase through cAMP signaling pathway to control the estrogen level 35 . EGR1 has been reported to be one of the target genes of fibroblast growth factors (FGF) in cow GCs 36 and it can promote the apoptosis of GCs in aging mouse ovary 37 . In our another study (unpublished data) MEGE8 can be used as a marker gene in goat GCs. Many of these marker genes need to be identified in the future.
Next, we tried to explore the functional gene expression in oocytes and GCs in the different populations (high fertility group, HL; low fertility group, LL). We found much more genes are highly expressed in large and medium oocytes of HL population than that in LL population, while, more gene expressed in small oocytes of LL population than that in HL population. The enriched functional pathways in the high expressed gene in the large and medium oocytes of HL population are related to protein translation, metabolism, and transportation which are very close to the oocyte function at these stages. While the enriched functional pathways for the genes highly expressed in LL population are related to other functions. The data indicated that oocytes in HL population pose the high potential for folliculogenesis.
The gene expression profile poses the developmental trend in GCs, too. And many more genes were differentially expressed in GCs of the two populations (HL vs. LL) than that in oocytes. More genes were highly expressed in GCs of large follicles in HL group than that in LL group, while more genes were highly expressed in GCs of medium follicles in LL group than that in HL group. The enriched functional pathways in the highly expressed genes in large and medium follicular GCs were related to cell cycle and proliferation which is close to the function of GCs at these stages. The data suggested that the expression of genes in GCs were correlated to the oocyte for the folliculogenesis which may reflect the fertility difference in these populations.
It has been reported recently ligands and receptors mediate the regulatory signaling between follicular oocytes and GCs during folliculogenesis (Biase and Kimble 2018). And the transcriptome profile of oocytes 18-21 and GCs 7-17 can reflect the developmental potential of oocyte to be successfully fertilized to develop to be embryo. Moreover, ligand-receptor pairs can tune the signaling for the regulation of gene expression (Biase and Kimble 2018). In current investigation, we found that ligand-receptor pairs played important roles in goat folliculogenesis. Notch signaling pathway has been reported to be one of the most important pathways in folliculogenesis 38 . DLL3, a ligand in Notch pathway, was highly expressed in oocytes, while another ligand in Notch pathways JAG2 was highly expressed in GCs compare with DLL3 and JAG1. The receptors in Notch pathway NOTCH2 and NOTCH3 were highly expressed in GCs while the target gene HES1 was highly expressed in oocytes. The transforming growth factor-β (TGF-β) superfamily plays vital roles in the tuning folliculogenesis 39 . The important ligands in TGF-β pathway were highly expressed in oocytes while the receptor BMPR1b was highly expressed in GCs. The target gene ID3 was highly expressed in oocytes while ID4 was highly expressed in GCs. Moreover, we found ligand-receptor pairs were differentially expressed in the oocytes and GCs of HL and LL at different developmental stages. The data suggested that the gene expression levels are differentially regulated by ligandreceptors in HL and LL population during folliculogenesis.
Furthermore, the correlations of the gene expression in oocytes and GCs at different stages in the two populations HL and LL were different, too. It has been reported that the association between the transcript levels of genes expressed in oocytes and GCs play important roles in follicogenesis, too 24 . Our data matched this notion very well.

Conclusions
In summary, in the current investigation we found many marker genes differentially expressed in oocytes or GCs during folliculogenesis. The ligand-receptor pairs in oocyte-GCs are involved in the regulation of the gene expression profile in two fertility population HL and LL. At the same time, the gene expression pattern was correlated well in oocyte and GCs in HL or LL. All the data reflect the gene expression landscape in oocytes and GCs which is correlated well with the fertility capability.