Construction of protein–protein interaction network based on transcriptome profiling of ovine granulosa cells during the sheep’s anestrus phase

Abstract

BACKGROUND: Small antral follicles as the final reserve of folliculogenesis, are existed throughout the reproductive life span of sheep. However, the ovarian cycles of ewe cease in the anestrus phase. This study was aimed to elucidate the ovarian small antral follicles transcriptome in the ewe’s anestrus phase.

METHODS: Granulosa cells of small antral follicles (≤3 mm) were collected from two groups of Mehraban ewes under long days of summer in the non-breeding season as anestrus phase. Transcriptome profiling of these granulosa cells were obtained using the RNA-Seq technology. An integrative analysis was utilized to identify key regulatory genes which may have potential impacts on intra-ovarian molecular activities.

RESULTS: Globally, 14506 genes were expressed whose higher expressions belonged to genes that encoded ribosomal proteins. Top significant terms of gene ontology were pertained to protein translational processes. Apart of this, most of highly significant terms were also relevant to apoptotic process through extracellular vesicles, including apoptotic bodies and exosomes. Regarding to node effect property, UBA52 (ubiquitin A-52 residue ribosomal protein fusion product 1) and RPS5 (ribosomal proteins S5) contained in highest out-degree and in-degree, respectively.

CONCLUSION: Our data suggest that ribosomal mRNA/proteins could make the granulosa cells undergo a lot of changes from the point of view that ovarian activities are ceased in the anestrus phase.

Keywords

RNA-sequencing / ovarian follicles / anestrus phase / protein-protein interaction network / functional modules

Introduction

Reproductive activities in ewe are characterized by a seasonality influenced by several factors such as photoperiod, latitude, temperature, nutrition and breed (^{Mbayahaga et al., 1998}). Photoperiod is the main factor that regulates the breeding season of wool sheep in high latitudes in the northern hemisphere (^{Hafez, 1952}). Patterns of reproductive function in the ewe are dominated by two distinct rhythms. The first of these is a 16- to 17-day estrous cycle. The other major rhythm in ovine reproduction is an annual cycle of ovarian activity. In most breeds of sheep, normal estrous cycles occur in the autumn and winter that this called breeding season. On the contrary, the ovarian cycles cease in the spring and summer that this called anestrus phase (^{Hafez, 1952}; ^{Bartlewski et al., 2011}; ^{Goodman and Inskeep, 2015}).

Each estrous cycle produces one or two mature oocytes for ovulation and fertilization. However, the oocyte cannot reach this point of development without the support of its follicle (^{Hennet and Combelles, 2012}). Follicle is composed of somatic cells and the developing oocyte. The two primary somatic cell types in the ovarian follicle are the theca cells and granulosa cells. These two somatic cell types are the site of action and synthesis of a number of hormones which promote a complex regulation of follicular development (^{Skinner et al., 2008}). Granulosa cells are the most important follicular somatic cells in the ovary that undergo serious changes morphologically and physiologically during the processes of proliferation, differentiation, ovulation, lutenization and atresia (Cuiling et al., 2005).

Cohort of small antral follicles enters to gonadotrophin-dependent development upon FSH stimulation, called recruited or ingrowing follicles (^{Shimizu, 2016}). Subsequently, during follicular selection, one or more of the recruited follicle(s) are selected to continue(s) maturation, ultimately being ovulated or becoming atretic (^{Noel et al., 1993}). Therefore, small antral follicles represent the fundamental developmental unit of the mammalian ovary and, as such, serve the needs of the entire reproductive life span. In spite, many small antral follicles are laid on the ovaries throughout the ovine reproductive activities (^{Fortune, 1994}), but ewes cannot mate during anestrus and no lambs can be born for several months (^{Di et al., 2014}). Importantly, the fate about cohort of small antral follicles has been located on the ovaries in the anestrus phases remained ambiguous.

In recent years, transcriptome analysis of ovaries has been leveraged to study questions related to reproductive biology of ovaries such as ovine early folliculogenesis (^{Bonnet et al., 2013}), bovine small atretic follicles (^{Hatzirodos et al., 2014a}), bovine small and large follicles (^{Hatzirodos et al., 2014b}), porcine with high and low litter size (^{Zhang et al., 2015}), distinct sizes of bovine follicles (^{Labrecque et al., 2016}), ovarian stages in pigeon (^{Xu et al., 2016a}), bovine follicular and luteal somatic cells (^{Romereim et al., 2017}), porcine small atretic follicles (^{Terenina et al., 2017}), and in our previous study on ovine granulosa cells of follicular and luteal phases during the estrous cycle (^{Talebi et al., 2018b}). Therefore, there is a gap in molecular knowledge regarding the recruitment of small antral follicles in the non-breeding season.

The purpose of this study was mainly focused to thoroughly investigate the global transcriptome (mRNA levels) of ovarian granulosa cells from the small antral follicles in anestrus phase of Mehraban sheep by RNA sequencing method. As well as, an integrative analysis was utilized to elucidate key regulatory genes whose may have potential impacts on intra-ovarian molecular activities of the granulosa cells in non-breeding season of sheep.

Materials and methods

Animals and experimental design

Two non-pregnant ewes of the Mehraban breed were selected for this experiment from the research farm of Bu-Ali Sina University, Hamedan, Iran. Each ewe was slaughtered in the anestrus phase during the non-breeding season on 21 June. Ovaries were dissected from the oviducts by creaser and then, two separated ovaries were suspended into normal saline solution (NSS) fortified with antibiotics and carried to the genetic laboratory of Bu-Ali Sina University at 35-37°C within 3 h of collection.

Collection of granulosa cells

Ovaries were washed 3-4 times in NSS and rinsed in 70% ethyl alcohol for 2-3 min to eliminate surface organisms. Follicular fluid of small antral follicles (≤3 mm) was finely aspirated via suction pump of vacuum to collection of mural granulosa cells of each ovary. For each group 6 follicles roughly were pooled by animal and kept for granulosa cell isolation into medium culture of Fetal Bovine Serum (10%) with heparin (10 mg/mL). The granulosa cells were removed from the remainder of the follicles by gentle rubbing with a glass Pasteur pipette, previously modified by heat sealing the tip into a rounded smooth surface. The medium culture containing the granulosa cells was centrifuged at 500 × g for 7 min at 4°C, the medium was removed by aspiration and the cells washed twice in phosphate-buffered saline (PBS) cold, pH 7.4. Finally the cells were resuspended in RNA Stabilization Reagent (QIAGEN, Cat No./ID: 76154), and stored at - 20°C for subsequent RNA isolation.

RNA extraction and assessment of quantity, purity and quality

Total RNA was isolated from granulosa cells with RNeasy Protect Mini Kit (QIAGEN, Cat No./ID: 74126) including a DNAse1 treatment step and following the manufacturer instruction. The RNA purity was determined via the OD₂₆₀/OD₂₈₀ ratio and OD₂₆₀/OD₂₃₀ ratio by the NanoDrop® ND-8000. For such verification, RNA samples which extracted from two groups of granulosa cells were also resolved with 0.7% agarose gel electrophoresis. The quality of the RNA was assessed using a Fragment Analyzer™ (Advanced Analytical Technologies) with a RNA Quality Number (RQN).

Library preparation and high-throughput RNA sequencing

All RNA sequencing libraries were constructed from 1 mg of total RNA using Illumina TruSeq RNA Sample Prep (Illumina) according to the manufacturer’s instructions. Libraries were quantified using the qPCR NGS Library Quantification Kit (Agilent). Each cDNA library was tagged with a specific Illumina adapter sequence and the two libraries were subsequently pooled and loaded on one sequencing lane. Sequencing was performed on an Illumina HiSeq 3000 apparatus using the Illumina TruSeq Standard mRNA kit v2 to obtain paired-end reads (2 × 150 bp).

Bioinformatics analysis

The barcoding tags were trimmed in order to separate the two samples. Quality control of Illumina reads in FastQ format was checked using Fastqc v0.11.2 software. The valid Illumina paired reads in FastQ format, which can be used for subsequent analysis, were obtained by removing adapter sequences and low quality raw reads. Four cleaned fastq files and metadata were released in publicly available database: Sequence Read Archive (SRA): SRP148599: https://www.ncbi.nlm.nih.gov/sra/SRP148599. The read mapping, counting and statistical analysis were performed using the local instance of Galaxy (https://galaxyproject.org) at the Toulouse Midi-Pyrénées bioinformatics platform (http://sigenae-workbench.toulouse.inra.fr/), France. The cleaned paired reads were merged and mapped against the ovine genome assembly (Oar_v3.1.84) using STAR2 (STAR_2.4.0i, Galaxy tool v2.0). Obtained BAM files were sorted using Samtools_sort (Galaxy tool V1.0). Rmdup (v1.0.0) (^{Li et al., 2009}) was used to remove potential PCR duplicates. A GFF3 annotation file was obtained from Ensembl (Ovis_aries.Oar_V3.1.84) and converted to GFF using the GFF-to-GTF tool (version 2.0.0). Reads counting for annotated transcripts, genes and isoforms (FPKM values, fragments per kilobase of transcript per million fragments mapped) for each sample was done using SigCufflinks (version 1.0.0) adapted from the Cufflink software (^{Trapnell et al., 2010}).

Enrichment analysis of Gene Ontology (GO) pathways

ClueGO plugin v2.1.7 (http://www.ici.upmc.fr/cluego/) (^{Bindea et al., 2009}) of Cytoscape V. 3.1.1 was used to create and visualize a functionally grouped network of biological processes (BP), cellular components (CC), and molecular functions (MF) from the genes whose high expressed (100≤FPKM) in ovarian granulosa cells in the ewe's anestrus phase. These were accumulated from the gene and protein symbols of HGNC (HUGO Gene Nomenclature Committee). ClueGO uses precompiled annotation files including GO, KEGG and BioCarta for a wide range of organisms (^{Bindea et al., 2009}). Additionally, ClueGO can easily integrate new annotation sources in a plugin like way (^{Bindea et al., 2009}). Pathway-based analysis further informed the biological functions suggesting significantly enriched metabolic or signal transduction pathways associated to the highest expressed genes, compared to the whole genomic background (^{Zhang et al., 2014}). Term enrichment was tested with a right-sided hyper-geometric test that was corrected for multiple testing by the Bonferroni stepdown method (^{Holm, 1979}). Only GO/pathway terms that were significantly enriched (corrected P-value<0.01) were included in the analysis. Kappa statistics were used to link and grouping of the enriched terms and functional grouping of them as described in (^{Bindea et al., 2009}). The minimum connectivity of the pathway network (kappa score) was set to 0.4 units.

Protein–protein interaction (PPI) network

A large PPI (protein–protein interaction) network was reconstructed through high expressed genes (100≤FPKM) at anestrus phase. PPI network was consisted of nodes (protein symbols of HGNC) and edges (interaction between identified nodes). To map pairwise interactions, all computational methods in STRING database Version. 10.0 containing neighborhoods, gene fusion, co-occurrence, homology, coexpression, experiments, databases and text mining were utilized with the highest confidence score≥0.9 (^{Szklarczyk et al., 2015}). This PPI network was visualized in Cytoscape Version 3.1.1 (^{Lotia et al., 2013}), which furnishes diverse plugins for multiple analyses.

Network analysis and modules selection

Cytoscape Version 3.1.1 (^{Lotia et al., 2013}) was used to plot and analyze the centralities, clustering and modularity of the network. The MCODE (Molecular Complex Detection) v1.4.0-beta2 was performed to screen modules of PPI network with degree cutoff= 2, node score cutoff= 0.2, k-core= 2, and max. depth= 100. This is a well-known automated algorithm to find highly interconnected subgraphs that detects densely connected regions in large PPI networks that may represent molecular complexes (^{Bader and Hogue, 2003}). In thus MCODE detects genes/proteins complexes whose are with the highest quality, in terms of the function and localization similarity of proteins within predicted complexes. Also, the functional modules were chosen with the number of nodes and scores (^{Xu et al., 2016b}). In scale-free biological networks (^{Albert, 2005}), nodes having large number of interacting partners represent hubs in the network. Hubs were detected by calculating the node degree distribution (^{Malik et al., 2015}) using the Network Analyzer (http://apps.cytoscape.org/apps/networkanalyzer) plugin of Cytoscape (^{Bindea et al., 2009}). The other topological parameters like network centrality options such as stress, betweenness and closeness centralities were taken as subsequent indexes for collection of hubs.

Result

Global analysis of transcriptome sequencing

The transcriptome of ovine granulosa cells of small antral follicles (≤3 mm) were analyzed by global RNA sequencing of samples collected from two Mehraban ewes (ovaries revealed without corpus lutea) during anestrus phase at the non-breeding season. The RNA quantity obtained from the two groups of granulosa cells as G2 and G3 were 357.8 and 366.22, respectively, and RNA purity with an absorbance ratio (A260/280 nm) were determined 2.04 and 2.1, respectively. This revealed adequate amounts of RNA for library construction and the sequencing. For the present next generation sequencing (NGS) assay, RQN scores for G2 and G3 were determined in 7.2 and 6.1, respectively (Fig. 1, part a,b). RNA integrity number (RIN) scores≥7 are considered acceptable for the most NGS assays (^{Copois et al., 2007}). Also, separated 28s and 18s rRNAs (rRNAs), were clearly observed after 20 min migration on agarose gel electrophoresis (Fig. 1, part c).

More than 94% of the clean-read data achieved the Q30 level of the Phred-like quality score (sequencing error of approximately 0.001). Following quality filtering, RNA sequencing experiment generated around 145 million of 150 bp reads. Approximately 24 million reads (81% of the total reads) were uniquely mapped onto the version 3.1 of the ovine genome (Ovis_aries.Oar_V3.1.84), and aligned with 27057 covering exons, introns, and intergenic regions (Table 1). Gene expression was quantified through FPKM conversion of raw reads. According to the 0.3 FPKM threshold suggested by (^{Ramsköld et al., 2009}), 14506 genes (mean FPKM of the 2 samples>0.3) were considered as expressed in ovine granulosa cells that, 1777 genes were identified as unknown over the gene annotations provided by Ensembl (Ovis_aries.Oar_V3.1.84). As shown in Fig. 2, in separately, 13794, 530, 251, 95 and 16 genes were categorized as autosome, X-linked (sexual), JH, AMGL and MT (mitochondrial), respectively (Fig. 2). The highest expression were belonged to genes that encoding ribosomal proteins (Table 2). To ensure that the cells isolated in anestrus phase were not contaminated with any thecal cells, no follicles with level of expression of CYP17A1 found in thecal samples were included. In fact, CYP17A1 is expressed exclusively in thecal cells (^{Rodgers et al., 1986}).

Functional enrichment of gene ontology pathways

In this study, 436 genes with higher expression level of FPKM (100≤FPKM) in ovine granulosa cells of small antral follicles from the anestrus phase, were analyzed for the gene ontology significant pathways (corrected with Bonferroni P-value<0.01). These genes were categorized into 246, 70 and 35 terms for biological processes (BP), cellular components (CC) and molecular functions (MF), respectively (Additional file 1). The top significant terms for BP, CC and MF, were belonged to protein localization to endoplasmic reticulum, extracellular vesicle and RNA binding, respectively (Table 3). Approximately, most of highly significant terms (P-value<0.01) of BP were relevant to apoptotic process and cell death (Additional file 1). Interestingly, these process were observed with most genes were enriched in cellular components (CC) of membrane-bounded vesicle (n = 229), extracellular organelle (n = 217) and extracellular vesicle (n = 217) (Additional file 1).

Construction of PPI network and topological analysis

Among 436 genes (100≤FPKM) from HUGO symbols, 292 proteins as nodes (highest confidence score≥0.9) were annotated on the PPI networks, and 3511 edges containing neighborhoods, gene fusion, co-occurrence, homology, coexpression, experiments, databases and text mining, were interacted between such genes/proteins (Additional file 2). The statistics of this PPI network including network density, network diameter and clustering coefficient, were 0.083, 9 and 0.643, respectively. The power law of degree distribution was followed with an R² = 0.946 for the PPI network. Meanwhile, R² is computed on logarithmized values. As shown in Fig. 3, a total PPI network was created from the 292 genes/proteins, that genes were clustered in different categories after using layout option of edge weighted spring embedded (Fig. 3). Top ten of genes/proteins with highest centrality degree in the PPI network were belonged to RPS27A, UBA52, RPL15, RPS15, RPS6, RPS8, RPL11, RBM8A, RNPS1 and RPS3, respectively (Fig. 3 and Additional file 2). Among these, RPS27A, UBA52, RBM8A RNPS1 and RPL15, had high values in all degree, stress, betweenness and closeness centralities (Additional file 2). These are crucial to identify some candidate genes among hubs as major regulators of small antral follicles in the sheep’s anestrus phase.

Identification of functional modules

Total PPI network was deciphered for identification of functional modules. In totally 12 modules were extracted among the total PPI network (Table 4). These modules were obtained using the MCODE plugin implementing the MCODE algorithm with the k-core value of 2.0, node score cutoff of 0.2, maximum depth from the seed node of 100 and graphics-processing- unit-based parallelization to find modules efficiently (^{Kotni et al., 2016}). Which three functional modules were detected with the intra-connection nodes≥16 and node score>10 (Table 4 and Fig. 3). These functional modules whose shown in Fig. 3 and Table 4, including module 1 (n = 73, score= 68.2), module 2 (n = 21, score= 21) and module 3 (n = 16, score= 12), were pertained to genes encoded ribosomal proteins, ribonucleoproteins and mitochondrial NADH dehydrogenases, respectively (data has not been shown).

Central-directed PPI network of hub genes/proteins

To identify candidate genes among discovered hubs with highest degree centrality, top 30 genes/proteins were extracted from the PPI network. Interestingly, these hubs were thoroughly involved in the module 1, indicating the reliability of the modules model. According to central-directed PPI network, UBA52 (Ubiquitin A-52 Residue Ribosomal Protein Fusion Product 1) and RPS5 (ribosomal proteins S5) contained in highest out-degree and in-degree, respectively (Additional file 3). Therefore, UBA52 is known as an upstream regulator can regulate 79 genes/proteins in the PPI network. By contrast, RPS5 is known as a downstream regulator has been influenced by 74 other genes/proteins (Additional file 3). These two candidate genes were also become a highest degree centrality, which belonged to module 1 in red rectangles (Fig. 3).

Discussion

In the present study, the global transcriptome profile of ovarian granulosa cells in sheep’s anestrus phase was determined for the first time. This study documented the expression of 14506 genes in granulosa cells of small antral follicles in sheep’s anestrus phase. These high numbers of expressed genes were close to the higher number of genes expressed in different RNA-seq studies on oocyte and granulosa cells (14172 and 14087 genes, respectively) by (^{Bonnet et al., 2013}). Cytochrome c oxidase subunit 1 (COX1) revealed a top expression level (FPKM= 3219.08) in small antral follicles of anestrus phase (Table 2). This is the terminal component and one of the three genomic components of the mitochondrial respiratory chain, plays an important role in oogenesis and follicle maturation (^{Zhen et al., 2015}). In addition to the above, the expression of follicle stimulating hormone receptor (FSHR) was interestingly observed into granulosa cells of small antral follicles in sheep’s anestrus phase (FPKM= 62.97, data has not been shown). One of the main roles of FSH is to promote antral follicular development but also to protect ovarian follicles from atresia (^{Peluso and Steger, 1978}). FSHR gene has also been known as potential biomarker of granulosa cells of antral follicles (^{Xu et al., 2011}; ^{Hatzirodos et al., 2014a},^b; ^{Romereim et al., 2017}). However, its expression is not notable in comparison with higher expressed genes as shown in Table 2. Because, many top significant pathways of apoptotic process and cell death for higher expressed genes were identified in sheep’s anestrus phase (Additional file 1). Therefore, this reveals that the granulosa cells of small antral follicles are nominated in response to the follicle stimulating hormone (FSH) no luteinizing hormone (LH) (^{Douville and Sirard, 2014}).

In previous study, we surprisingly have shown differences in transcriptome profiles between small antral follicles collected during the follicular and luteal phases. Furthermore, we revealed the possible initiation of early follicular atresia in small antral follicles during the follicular phase in interaction with the presence of immune cells (^{Talebi et al., 2018b}). In this study, the enrichment analyses have been demonstrated that the translational processes were overrepresented amid genes with higher expression in the anestrus phase. These genes had functions into cellular components of extracellular vesicles or exosomes (Table 3). There are other types of microvesicle, including apoptotic bodies and ectosomes during apoptosis has been long known (^{Cocucci et al., 2009}; ^{György et al., 2011}; ^{Edgar, 2016}). In this study, 18 biological pathways of higher expressed genes that has been enriched in apoptotic signaling pathway and cell death (Additional file 1), likely attributed to the apoptotic bodies and extracellular vesicles. This has been reported the reason of regression and elimination of pigeon ovulated follicles by (^{Xu et al. 2016a}). Exosomes can carry the MHC–peptide complexes that are recognized by T lymphocytes (^{Raposo et al., 1996}) and that secretion of such exosomes could promote anti-tumor immune responses in mice in vivo (^{Zitvogel et al., 1998}). Therefore, these are some potential reasons reveal despite expression of FSHR but many higher significant terms are attributed to the apoptosis of granulosa cells in the sheep’s anestrus phase.

Protein–protein interaction (PPI) networks are one the most known approaches for representation of candidate biomarkers beyond of high-throughput studies (^{Malik et al., 2015}; ^{Talebi et al., 2016}; ^{Xu et al., 2016b}; ^{Bahrami et al., 2017}; ^{Talebi et al., 2018a}). Following the PPI network, top 70 genes/proteins with the high centrality degree were belonged to ribosomal proteins encoded (Additional file 2). These also have been highly expressed (level of FPKM) into granulosa cells of small antral follicles in the sheep’s anestrus phase (Additional file 1). In this study, MCODE clustering algorithm revealed 12 modules have been derived from the PPI network, which three were highly interconnected (Table 4). A module is a group of closely related proteins that act in concert to perform specific biological functions in time and space (^{Lin et al., 2015}). As shown in Table 4, the most nodes (73 genes/proteins) with higher score (68.2) belonged to module 1. The PPI evolution scores were utilized to describe the module elements, including core and ring components. Functions of core components were highly correlated with those of essential genes (^{Lin et al., 2015}). Additionally, we utilized several centrality parameters including stress, betweenness and closeness instead of using degree centrality itself to identify the hubs genes. Among 292 genes/proteins, 30 were chosen as hubs in highest degree and some centrality parameters (Additional file 2), that thoroughly were belonged to the module 1 as known a subnetwork among 73 genes encoded ribosomal proteins (Fig. 3). Therefore, in agreement with (^{Bonnet et al., 2008}), the first subnetwork brought together the genes coding for ribosomal proteins. A number of ribosomal proteins (RPs) have been found to have a role in inducing the apoptosis and cell death (^{Chen and Ioannou, 1999}; ^{Naora et al., 1995}; ^{Lai and Xu, 2007}; ^{Zhou et al., 2015}; ^{Qi et al., 2017}).

Following the topological parameters in the central-directed PPI network, UBA52 has been identified as an upstream regulator with highest out-degree. By contrast, RPS5 with the highest in-degree well known as downstream regulator (Additional file 3). Possibly, these two genes can be potential biomarkers whose play critical roles in ovine granulosa cells of small antral follicles in the anestrus phase. It has been reported UBA52 has roles in regulation of FSHR for normal folliculogenesis of female ovaries (^{Cohen et al., 2003}). It has been shown that the UBA52 mRNA is preferentially overexpressed during hepatoma cell apoptosis (^{Han et al., 2012}), and in human colorectal carcinoma (^{Barnard et al., 1995}). Interestingly, Kobayashi et al. have discovered UBA52 as a dual regulator of ribosomal protein complex whose regulates cell fate and embryonic development (^{Kobayashi et al., 2016}). Terminal follicular growth of swine was notably accompanied with decreased expression of genes implicated in protein translation including subunits of ribosomal proteins (^{Bonnet et al., 2008}). Furthermore, it has been demonstrated the potential role of RPS5 gene expression and maybe RPS5 protein on the initiation of erythroid differentiation of murine erythroleukemia cells and cell cycle arrest (^{Matragkou et al., 2008}; ^{Vizirianakis et al., 2015}). Therefore, evidences existed to indicate that some ribosomal proteins besides being structural components of the ribosomal subunits are involved in the regulation of cell differentiation and apoptosis.

Conclusion

Despite the expression of FSHR into ovine granulosa cells of small antral follicles in the anestrus phase, many higher expressed genes have been enriched in the apoptosis terms. These events roughly ascribed to genes which encoded the ribosomal proteins. Among 12 subnetworks deviated from the PPI network, three which were identified as functional modules. Two genes of UBA52 and RPS5 were identified as potential regulators of ribosomal proteins among genes in sheep’s anestrus phase. Possibly, ribosomal mRNA/proteins could make the granulosa cells undergo a lot of changes during the ovarian activities are ceased in the anestrus phase. Our findings provided unique molecular insights into mural granulosa cells of small antral follicles in the sheep’s anestrus phase.

References

[1] Albert R (2005). Scale-free networks in cell biology. J Cell Sci, 118(21): 4947–4957

[2] Bader G D, Hogue C W (2003). An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinformatics, 4(1): 2

[3] Bahrami A, Miraie-Ashtiani S R, Sadeghi M, Najafi A (2017). miRNA-mRNA Network involved in folliculogenesis interactome: Systems biology approach. Reproduction, 154(1): 51–65

[4] Barnard G F, Mori M, Staniunas R J, Begum N A, Bao S, Puder M, Cobb J, Redman K L, Steele G D Jr, Chen L B (1995). Ubiquitin fusion proteins are overexpressed in colon cancer but not in gastric cancer. Biochim Biophys Acta, 1272(3): 147–153

[5] Bartlewski P M, Baby T E, Giffin J L (2011). Reproductive cycles in sheep. Anim Reprod Sci, 124(3–4): 259–268

[6] Bindea G, Mlecnik B, Hackl H, Charoentong P, Tosolini M, Kirilovsky A, Fridman W H, Pagès F, Trajanoski Z, Galon J (2009). ClueGO: a Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics, 15(8): 1091–1093

[7] Bonnet A, Cabau C, Bouchez O, Sarry J, Marsaud N, Foissac S, Woloszyn F, Mulsant P, Mandon-Pepin B (2013). An overview of gene expression dynamics during early ovarian folliculogenesis: specificity of follicular compartments and bi-directional dialog. BMC Genomics, 14(1): 904

[8] Bonnet A, Lê Cao K A, Sancristobal M, Benne F, Robert-Granié C, Law-So G, Fabre S, Besse P, De Billy E, Quesnel H, Hatey F, Tosser-Klopp G (2008). In vivo gene expression in granulosa cells during pig terminal follicular development. Reproduction, 136(2): 211–224

[9] Chen F W, Ioannou Y A (1999). Ribosomal proteins in cell proliferation and apoptosis. Int Rev Immunol, 18(5-6): 429–448

[10] Cocucci E, Racchetti G, Meldolesi J (2009). Shedding microvesicles: artefacts no more. Trends Cell Biol, 19(2): 43–51

[11] Cohen B D, Bariteau J T, Magenis L M, Dias J A (2003). Regulation of follitropin receptor cell surface residency by the ubiquitin-proteasome pathway. Endocrinology, 144(10): 4393–4402

[12] Copois V, Bibeau F, Bascoul-Mollevi C, Salvetat N, Chalbos P, Bareil C, Candeil L, Fraslon C, Conseiller E, Granci V, Mazière P, Kramar A, Ychou M, Pau B, Martineau P, Molina F, Del Rio M (2007). Impact of RNA degradation on gene expression profiles: assessment of different methods to reliably determine RNA quality. J Biotechnol, 127(4): 549–559

[13] Cuiling L, Wei Y, Zhaoyuan H, Yixun L (2005). Granulosa cell proliferation differentiation and its role in follicular development. Chin Sci Bull, 50(23): 2665–2671

[14] Di R, He J, Song S, Tian D, Liu Q, Liang X, Ma Q, Sun M, Wang J, Zhao W, Cao G, Wang J, Yang Z, Ge Y, Chu M (2014). Characterization and comparative profiling of ovarian microRNAs during ovine anestrus and the breeding season. BMC Genomics, 15(1): 899

[15] Douville G, Sirard M A (2014). Changes in granulosa cells gene expression associated with growth, plateau and atretic phases in medium bovine follicles. J Ovarian Res, 7(1): 50

[16] Edgar J R (2016). Q&A: What are exosomes, exactly? BMC Biol, 14(1): 46

[17] Fortune J E (1994). Ovarian follicular growth and development in mammals. Biol Reprod, 50(2): 225–232

[18] Goodman R L, Inskeep E K (2015). Control of the ovarian cycle of the sheep. In: Plant T M, Zeleznik A J, editors. Knobil and Neill’s Physiology of Reproduction 4. San Diego, CA, USA: Elsevier Inc, pp. 1259–1305

[19] György B, Szabó T G, Pásztói M, Pál Z, Misják P, Aradi B, László V, Pállinger E, Pap E, Kittel A, Nagy G, Falus A, Buzás E I (2011). Membrane vesicles, current state-of-the-art: emerging role of extracellular vesicles. Cell Mol Life Sci, 68(16): 2667–2688

[20] Hafez E (1952). Studies on the breeding season and reproduction of the ewe Part I. The breeding season in different environments Part II. The breeding season in one locality. J Agric Sci, 42(03): 189–231

[21] Han X J, Lee M J, Yu G R, Lee Z W, Bae J Y, Bae Y C, Kang S H, Kim D G (2012). Altered dynamics of ubiquitin hybrid proteins during tumor cell apoptosis. Cell Death Dis, 3(1): e255

[22] Hatzirodos N, Hummitzsch K, Irving-Rodgers H F, Harland M L, Morris S E, Rodgers R J (2014a). Transcriptome profiling of granulosa cells from bovine ovarian follicles during atresia. BMC Genomics, 15(1): 40

[23] Hatzirodos N, Irving-Rodgers H F, Hummitzsch K, Harland M L, Morris S E, Rodgers R J (2014b). Transcriptome profiling of granulosa cells of bovine ovarian follicles during growth from small to large antral sizes. BMC Genomics, 15(1): 24

[24] Hennet M L, Combelles C M (2012). The antral follicle: a microenvironment for oocyte differentiation. Int J Dev Biol, 56(10-11-12): 819–831

[25] Holm S (1979). A simple sequentially rejective multiple test procedure. Scand J Stat, 6: 65–70

[26] Kobayashi M, Oshima S, Maeyashiki C, Nibe Y, Otsubo K, Matsuzawa Y, Nemoto Y, Nagaishi T, Okamoto R, Tsuchiya K, Nakamura T, Watanabe M (2016). The ubiquitin hybrid gene UBA52 regulates ubiquitination of ribosome and sustains embryonic development. Sci Rep, 6(1): 36780

[27] Kotni M K, Zhao M, Wei D Q (2016). Gene expression profiles and protein-protein interaction networks in amyotrophic lateral sclerosis patients with C9orf72 mutation. Orphanet J Rare Dis, 11(1): 148

[28] Labrecque R, Fournier E, Sirard M A (2016). Transcriptome analysis of bovine oocytes from distinct follicle sizes: Insights from correlation network analysis. Mol Reprod Dev, 83(6): 558–569

[29] Lai M D, Xu J (2007). Ribosomal proteins and colorectal cancer. Curr Genomics, 8(1): 43–49

[30] Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R (2009). The sequence alignment/map format and SAMtools. Bioinformatics, 25(16): 2078–2079

[31] Lin C Y, Lee T L, Chiu Y Y, Lin Y W, Lo Y S, Lin C T, Yang J M (2015). Module organization and variance in protein-protein interaction networks. Sci Rep, 5(1): 9386

[32] Lotia S, Montojo J, Dong Y, Bader G D, Pico A R (2013). Cytoscape app store. Bioinformatics, 29(10): 1350–1351

[33] Malik A, Lee E J, Jan A T, Ahmad S, Cho K H, Kim J, Choi I (2015). Network Analysis for the Identification of Differentially Expressed Hub Genes Using Myogenin Knock-down Muscle Satellite Cells. PLoS One, 10(7): e0133597

[34] Matragkou C N, Papachristou E T, Tezias S S, Tsiftsoglou A S, Choli-Papadopoulou T, Vizirianakis I S (2008). The potential role of ribosomal protein S5 on cell cycle arrest and initiation of murine erythroleukemia cell differentiation. J Cell Biochem, 104(4): 1477–1490

[35] Mbayahaga J, Mandiki S N, Bister J L, Paquay R (1998). Body weight, oestrous and ovarian activity in local Burundian ewes and goats after parturition in the dry season. Anim Reprod Sci, 51(4): 289–300

[36] Naora H, Nishida T, Shindo Y, Adachi M, Naora H (1995). Association of nbl gene expression and glucocorticoid-induced apoptosis in mouse thymus in vivo. Immunology, 85: 63–68

[37] Noel B, Bister J L, Paquay R (1993). Ovarian follicular dynamics in Suffolk ewes at different periods of the year. J Reprod Fertil, 99(2): 695–700

[38] Peluso J J, Steger R W (1978). Role of FSH in regulating granulosa cell division and follicular atresia in rats. J Reprod Fertil, 54(2): 275–278

[39] Qi Y, Li X, Chang C, Xu F, He Q, Zhao Y, Wu L (2017). Ribosomal protein L23 negatively regulates cellular apoptosis via the RPL23/Miz-1/c-Myc circuit in higherrisk myelodysplastic syndrome. Sci Rep, 7(1): 2323

[40] Ramsköld D, Wang E T, Burge C B, Sandberg R (2009). An abundance of ubiquitously expressed genes revealed by tissue transcriptome sequence data. PLOS Comput Biol, 5(12): e1000598

[41] Raposo G, Nijman H W, Stoorvogel W, Liejendekker R, Harding C V, Melief C J, Geuze H J (1996). B lymphocytes secrete antigen-pre-senting vesicles. J Exp Med, 183(3): 1161–1172

[42] Rodgers R J, Rodgers H F, Hall P F, Waterman M R, Simpson E R (1986). Immunolocalization of cholesterol side-chain-cleavage cytochrome P-450 and 17 alpha-hydroxylase cytochrome P-450 in bovine ovarian follicles. J Reprod Fertil, 78(2): 627–638

[43] Romereim S M, Summers A F, Pohlmeier W E, Zhang P, Hou X, Talbott H A, Cushman R A, Wood J R, Davis J S, Cupp A S (2017). Gene expression profiling of bovine ovarian follicular and luteal cells provides insight into cellular identities and functions. Mol Cell Endocrinol, 439: 379–394

[44] Shimizu T (2016). Molecular and cellular mechanisms for the regulation of ovarian follicular function in cows. J Reprod Dev, 62(4): 323–329

[45] Skinner M K, Schmidt M, Savenkova M I, Sadler-Riggleman I, Nilsson E E (2008). Regulation of granulosa and theca cell transcriptomes during ovarian antral follicle development. Mol Reprod Dev, 75(9): 1457–1472

[46] Szklarczyk D, Franceschini A, Wyder S, Forslund K, Heller D, Huerta-Cepas J, Simonovic M, Roth A, Santos A, Tsafou K P, Kuhn M, Bork P, Jensen L J, von Mering C (2015). STRING v10: protein–protein interaction networks, integrated over the tree of life. Nucleic Acids Res, 43(Database issue): D447–D452

[47] Talebi R, Ahmadi A, Afraz F, Abdoli R (2016). Parkinson’s disease and lactoferrin: Analysis of dependent protein networks. Gene Rep, 4: 177–183

[48] Talebi R, Ahmadi A, Afraz F (2018a). Analysis of protein-protein interaction network based on transcriptome profiling of ovine granulosa cells identifies candidate genes in cyclic recruitment of ovarian follicles. J Anim Sci Technol, 60: 11

[49] Talebi R, Ahmadi A, Afraz F, Sarry J, Plisson-Petit F, Genêt C, Fabre S (2018b). Transcriptome analysis of ovine granulosa cells reveals differences between small antral follicles collected during the follicular and luteal phases. Theriogenology, 108: 103–117

[50] Terenina E, Fabre S, Bonnet A, Monniaux D, Robert-Granié C, SanCristobal M, Sarry J, Vignoles F, Gondret F, Monget P, Tosser-Klopp G (2017). Differentially expressed genes and gene networks involved in pig ovarian follicular atresia. Physiol Genomics, 49(2): 67–80

[51] Trapnell C, Williams B A, Pertea G, Mortazavi A, Kwan G, van Baren M J, Salzberg S L, Wold B J, Pachter L (2010). Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol, 28(5): 511–515

[52] Vizirianakis I S, Papachristou E T, Andreadis P, Zopounidou E, Matragkou C N, Tsiftsoglou A S (2015). Genetic manipulation of RPS5 gene expression modulates the initiation of commitment of MEL cells to erythroid maturation: Implications in understanding ribosomopathies. Int J Oncol, 47(1): 303–314

[53] Xu F, Stouffer R L, Müller J, Hennebold J D, Wright J W, Bahar A, Leder G, Peters M, Thorne M, Sims M, Wintermantel T, Lindenthal B (2011). Dynamics of the transcriptome in the primate ovulatory follicle. Mol Hum Reprod, 17(3): 152–165

[54] Xu X, Zhao X, Lu L, Duan X, Qin H, Du X, Li G, Tao Z, Zhong S, Wang G (2016a). Transcriptomic analysis of different stages of pigeon ovaries by RNA-Sequencing. Mol Reprod Dev, 83(7): 640–648

[55] Xu Z, Zhou Y, Cao Y, Dinh T L, Wan J, Zhao M (2016b). Identification of candidate biomarkers and analysis of prognostic values in ovarian cancer by integrated bioinformatics analysis. Med Oncol, 33(11): 130

[56] Zhang J, Pan Z, Moloney S, Sheppard A (2014). RNA-Seq Analysis Implicates Detoxification Pathways in Ovine Mycotoxin Resistance. PLoS ONE9, 6: e99975

[57] Zhang X, Huang L, Wu T, Feng Y, Ding Y, Ye P, Yin Z (2015). Transcriptomic analysis of ovaries from pigs with high and low litter size. PLoS One, 10(10): e0139514

[58] Zhen X, Wu B, Wang J, Lu C, Gao H, Qiao J (2015). Increased Incidence of Mitochondrial Cytochrome C Oxidase 1 Gene Mutations in Patients with Primary Ovarian Insufficiency. PLoS One, 10(7): e0132610

[59] Zhou X, Hao Q, Zhang Q, Liao J M, Ke J W, Liao P, Cao B, Lu H (2015). Ribosomal proteins L11 and L5 activate TAp73 by overcoming MDM2 inhibition. Cell Death Differ, 22(5): 755–766

[60] Zitvogel L, Regnault A, Lozier A, Wolfers J, Flament C, Tenza D, Ricciardi-Castagnoli P, Raposo G, Amigorena S (1998). Eradication of established murine tumors using a novel cell-free vaccine: dendritic cell-derived exosomes. Nat Med, 4(5): 594–600

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature