Signalling pathways and mechanistic cues highlighted by transcriptomic analysis of primordial, primary, and secondary ovarian follicles in domestic cat

In vitro growth (IVG) of dormant primordial ovarian follicles aims to produce mature competent oocytes for assisted reproduction. Success is dependent on optimal in vitro conditions complemented with an understanding of oocyte and ovarian follicle development in vivo. Complete IVG has not been achieved in any other mammalian species besides mice. Furthermore, ovarian folliculogenesis remains sparsely understood overall. Here, gene expression patterns were characterised by RNA-sequencing in primordial (PrF), primary (PF), and secondary (SF) ovarian follicles from Felis catus (domestic cat) ovaries. Two major transitions were investigated: PrF-PF and PF-SF. Transcriptional analysis revealed a higher proportion in gene expression changes during the PrF-PF transition. Key influencing factors during this transition included the interaction between the extracellular matrix (ECM) and matrix metalloproteinase (MMPs) along with nuclear components such as, histone HIST1H1T (H1.6). Conserved signalling factors and expression patterns previously described during mammalian ovarian folliculogenesis were observed. Species-specific features during domestic cat ovarian folliculogenesis were also found. The signalling pathway terms “PI3K-Akt”, “transforming growth factor-β receptor”, “ErbB”, and “HIF-1” from the functional annotation analysis were studied. Some results highlighted mechanistic cues potentially involved in PrF development in the domestic cat. Overall, this study provides an insight into regulatory factors and pathways during preantral ovarian folliculogenesis in domestic cat.


RNA-sequencing.
Total RNA was extracted with NucleoSpin RNA Plus XS (Macherey-Nagel GmbH & Co) from the 9 samples following the manufacturer's instructions. Total RNA sample integrity and concentration was measured by Agilent High Sensitivity RNA ScreenTape Assay (Agilent 2200 TapeStation system, Agilent Technologies, Inc., Santa Clara, US). Sample concentrations for quantification were determined by Qubit 2.0 Fluorometer (Thermo Fisher Scientific Inc., Waltham, USA). During extraction, genomic DNA was removed by gDNase during on-column DNA digestion. The purified RNA was dissolved in 20 µL RNase-free water and stored at − 80 °C. The cDNA was generated and amplified in 9 cycles with the SMART-Seq v4 Ultra Low Input RNA Kit for Sequencing (Takara Bio, Inc. Mountain View, USA). Libraries were generated (Nextera XT DNA Library Prep Kit, Illumina, San Diego, USA) and sample quality was checked by Agilent 2100 Bioanalyzer (Agilent, Santa Clara, USA) and by Agilent High Sensitivity D1000 ScreenTape Assay (Agilent 2200 TapeStation system). Custom index primers were tagged to the libraries to allow for multiplexing during RNA-sequencing. Libraries were quantified, normalised based on measurements determined by Qubit 2.0 Fluorometer, and sequenced on the Illumina NextSeq 500 system (150 cycles). The raw RNA-Seq data was deposited and released with the BioProject accession number of PRJNA635095.
Raw data consisted of paired-end, double-indexed cDNA library sequencing reads. The PhiX adapter-ligated control was sequenced at a standard concentration of 5%. The de-multiplexed data received in fastq file format was quality checked with FastQC 23 and summarised with MultiQC 24 . Salmon 25 quantified transcripts from a FASTA file containing the reference transcriptome and FASTQ files containing the sequence reads. Transcript abundance, counts, and length were summarised by tximport 26 . The DESeqDataSet under three factor levels "type" determined differential gene expression for PrFs, PFs, and SFs with DESeq2 27 . Ovarian follicle type contrasts, PrF versus PF (PrF-PF) and PF versus SF (PF-SF), were designed. Differentially expressed genes (DEGs) were estimated from the un-normalised, paired-end fragments by the Independent Hypothesis Weighting 28 method with an alpha = 0.05, an adjusted P value < 0.05, and absolute log 2 fold-change of 1 (Supplementary Data S1). For quality control, the log fold-change shrinkage estimates of normalised data were visualised with heatmaps of Euclidean distances and with a principal component analysis (PCA). To assign Entrez gene identifiers (IDs) BioMart 29,30 was employed. The Entrez IDs were input into the web-based portal "the database for annotation, visualisation and integrated discovery" (DAVID) Bioinformatics Resources 6.8 (https:// david. ncifc rf. gov/) 31,32 . The tool DAVID was chosen to: (1) identify major gene groups (functional classification); (2) to elucidate enriched annotation terms (functional annotation chart and clustering, respectively); (3) and to provide an overview of gene annotations (functional annotation table) 31 . Importantly, the last update, DAVID 6.8, occurred in October 2016 and the data analysed here was submitted to DAVID 6.8 at the beginning of 2019. Researchers interested in enrichment analysis should consider a combination of other tools such as gProfiler 33 , Enrichr 34,35 , and/or Metascape 36,37 . However, these portals and others such as, PANTHER 38 , InterMine 39,40 , and GeneTrail2 41 may be relying on old knowledge bases 37 . Here, the functional annotation and enrichment analysis workflow is summarised as follows: analysis with DAVID 6.8 output gene ontology (GO) 42 and KEGG (Kyoto Encyclopedia of Genes and Genomes) orthology (KO) 43 terms which were categorised into functional annotation clusters. The functional annotation tool (FAT) was applied to the GO categories BP, CC, and MF which filtered out broad GO terms based on measured specificity. Default categories were unselected and "GOTERM_BP_FAT", "GOTERM_CC_FAT", "GOTERM_MF_FAT", and "KEGG_PATHWAY" were selected. An over-representative hypergeometric test on the resulting GO and KO terms was performed and the clustering options within DAVID were modified as follows: Similarity Term Overlap 3, Similarity Threshold 0.60, Initial Group Membership 3, Final Group Membership 3, Multiple Linkage Threshold 0.50, and an Enrichment Threshold EASE 0.2. The GO and KO terms with an arbitrary enrichment score (ES) > 0.1 were considered for further investigation. A false discovery rate (FDR), Q value (Q), corrected for multiple testing of P values of the enriched terms. A Q < 0.05 was considered as significant. For visualisation purposes the web-based portal Metascape (http:// metas cape. org) was utilised 36 . Metascape currently does not support the domestic cat but it does support human. As gene annotation databases are primarily compiled for human genes it is useful to use human orthologs for model organisms prior to analysis 36 . As a model organism, the domestic cat shares several female germ cell features to humans: (1) the oocyte proper and the germinal vesicle diameter is equivalent; (2) oocytes reach metaphase II (MII) stage of meiosis after 24 h in culture; and (3) the nuclear configuration is comparable 44,45 . Thus, g:Profiler g:Orth orthology (https:// biit. cs. ut. ee/ gprofi ler/ orth) 33 converted the Felis catus Entrez IDs into Homo sapiens Ensembl IDs (an orthologous species recognisable by Metascape). The "Express Analysis" in Metascape was selected which automatically removed ontology terms that did not satisfy the minimal statistical criteria. The analysis report produced functional annotation network cluster graphs and ontology enrichment bar plots which Gene expression by qRT-PCR. Total RNA extraction for qRT-PCR was performed as described previously. Reverse transcription of RNA into cDNA was performed as described by Hryciuk et al. 46 . Intron-spanning primers were designed according to gene sequences listed in GenBank (Supplementary Table S1). The cDNA was diluted 1:10 and analysed with the CFX96 Real-Time PCR detection system utilising the SsoFast EvaGreen Supermix (Bio-Rad Laboratories GmbH, Munich), at 98 °C and 8 s at different annealing temperatures. Quantification of qRT-PCR products was performed with CFX Manager Software 3.1 (Bio-Rad Laboratories GmbH). Serial dilutions of plasmid DNA carrying genes of interest sequences or of qRT-PCR products were utilised for calibration with β-actin (BACT ) as a reference gene 47 . Statistical analysis was performed in R using R 3.4.4 22  Ethics declarations. Domestic cat ovaries were obtained from animal shelters in Berlin after routine ovariectomy for the purpose of permanent contraception primarily from stray domestic cat individuals. The surgical procedures were not related to the purpose of the experiment. The animal shelter agreed to donate the excised ovaries which we utilised as samples. Castrations are compliant with the "Protection of Animals Act" in Germany; no further guidelines had to be considered.  Supplementary Fig. S1a). Trimming of reads was not performed. Log count distributions between samples before and after normalisation were visualised ( Supplementary Fig. S1b). The PCA of the nine samples explained sample variance by the first two principal components (PC1 versus PC2) which accounted for 37% and 16% variation, respectively whereby, PC1 separated PrF from PF and SF, whereas PC2 separated PF from SF. This emphasised the sample-to-sample relationships based on ovarian follicle type ( Supplementary Fig. S2a). The PC1 versus PC3 for all genes accounted for 37% and 12% variation, respectively ( Supplementary Fig. S2b). The inclusion of PC3 explained a reduced proportion of variation between the PF and SF samples. Collectively, the first three PCs explained approximately 65% of total variation. The PC1 versus PC2 for 500 genes accounted for 39% and 19% variation, respectively ( Supplementary Fig. S2c). A hierarchical heatmap demonstrated that samples clustered together based on ovarian follicle type ( Supplementary Fig. S3). Clustering of all PrF samples was identified with a high degree of correlation. An outlier PF sample was observed (sample_4_S4) but the remaining two samples were closely associated. All SF samples were highly correlated. Overall, the degree of correlation demonstrated a relationship between PF and SF follicle samples. These samples were either less correlated or not correlated at all with the PrF samples ( Supplementary Fig. S3). Additionally, sample-to-sample relationships were visualised with hierarchical cluster dendograms ( Supplementary Fig. S4).

RNA
Differential gene expression analysis. The highest number of DEGs was identified during the PrF-PF transition. This included 2,226 DEGs comprised of 1,206 down-regulated and 1,020 up-regulated genes. The lowest number of DEGs was found in the PF-SF transition whereby, the total 154 DEGs estimate included 122 down-regulated and 32 up-regulated genes (Fig. 2a, c). This revealed that the number of DEGs decreased when comparing the PrF-PF and PF-SF transitions. The patterns of gene expression changes across all samples were graphically depicted in heatmap format (Fig. 2b) and summarised with a diagram of ovarian follicle DEGs (Fig. 2c). All differential gene expression lists were combined into one DEG database (Supplementary Data S1). Collectively, 2,380 DEG transcripts including 1,328 down-regulated and 1,052 up-regulated genes were detected in early preantral ovarian follicles in domestic cat. . The KO terms such as, "ECM receptor interaction" (Q = 1.60E−06) and "focal adhesion" (Q = 1.80E−03) were significantly enriched also. Other enriched terms and pathways were involved in "cellular response to growth factor stimulus" and the "phosphatidylinositol 3-kinase (PI3K)-Akt/Protein Kinase B (Akt)" and "transmembrane receptor protein serine/threonine kinase" signalling pathways (Table 1). Additionally, "platelet derived growth factor", "integrin-mediated", and "transforming growth factor (TGF)-β" signalling were observed during the PrF-PF transition (Supplementary Data S3). Interestingly, "axonemal dynein complex assembly" was identified in Annotation Cluster 32 and was coupled with "ciliary plasm" and "axoneme cellular components" (Supplementary Data S2). Functional annotation clustering of the 154 PF-SF DEGs resulted in 128 DAVID IDs in 27 clusters (Supplementary Data S2). Identified terms included "reproductive process", "cellular response to chemical stimulus", "ovulation cycle process", and "epithelial cell proliferation" ( Table 2). The KOs included terms involved in "focal adhesion" and "inositol phosphate metabolism" along with "sphingolipid", "phosphatidylinositol", "thyroid The estimated number of DEGs between ovarian follicle contrasts: primordial versus primary and primary versus secondary with an adjusted P value of < 0.05 and log 2 fold-change ≥ 1 were visualised utilising MA plots whereby, non-significant genes (values not satisfying the P value and log 2 fold-change threshold) are shown in black and significant genes (values satisfying the P value and log 2 fold-change threshold) are shown in red, log 2 foldchange is mapped to y and the normalised mean is mapped to x (transformed to log 10 scale); (b) Patterns of gene expression changes are shown on a heatmap with an adjusted P value of < 0.01 and log 2 fold-change ≥ 1. The colour red represents relative increase in expression, blue represents relative decrease, and pale orange (in the middle of the colour scale) represents no change. The columns are labelled with follicle type, respectively; (c) Pie charts summarise the estimated number of DEGs between the primordial versus primary and primary versus secondary contrasts with an adjusted P value of < 0.05 and log 2 fold-change ≥ 1. www.nature.com/scientificreports/ hormone", "calcium", "vascular endothelial growth factor (VEGF)", "ErbB" and "hypoxia-inducible factor 1 (HIF-1)" signalling pathways (Table 2). Additionally, "β-catenin", "apelin", "extra-nuclear estrogen", "IL8 CXCR1", "RAS", and "TXA2" signalling were identified. Lastly, "E2F", "forkhead box O (FOXO)", "Hedgehog", "IL6 JAK STAT3", "insulin-like growth factor (IGF)", "KIT", "Notch", and "Wnt" were identified though not significantly enriched (Supplementary Data S3). As described in the methods section, Metascape generated functional annotation network clusters and bar plots of enriched ontology terms for visualisation. The results from Metascape were extracted for PrF-PF ( Fig. 3 and Supplementary Fig. S5) and PF-SF ( Fig. 4 and Supplementary Fig. S6), respectively. Based on Metascape analysis, the functional annotation cluster and enrichment PrF-PF and PF-SF results were summarised in Excel format, respectively (Supplementary Data S3 and S4, respectively). The following pathways, processes, and regulators were selected for further analysis: PI3K-Akt; TGF-β; erythroblastoma (ErbB); HIF-1 signalling pathways; over-represented biological processes associated with the extracellular matrix (ECM); and potential structural regulators of ovarian folliculogenesis such as, the dynein complex.
Selected RNA-sequencing transcript expression. Supplementary Data S3 and S4 identified which genes were associated with the aforementioned selection, respectively. Several genes were selected (investigating both the significant and non-significant differential expression gene lists), transcript expression levels were plotted, and DESeq2 results were extracted for each ( Fig. 5 and Table 3). Some genes were selected due to functional implications related to ovarian folliculogenesis in other species. Whereas, others such as, bone morphogenetic proteins (BMPs) and MMPs were selected due to the piling evidence describing the essential role of these factors during mammalian ovarian folliculogenesis. For each numbered section below the following genes are listed in order of decreasing "baseMean" values: (1) PI3K-Akt signalling pathway: TSC complex subunit 2 (TSC2); pyruvate dehydrogenase kinase 4 (PDK4); AKT serine/threonine kinase 2 (AKT2); epidermal growth factor (EGF); regulatory associated protein of mTOR complex 1 (RPTOR); and PI3K subunit delta (PIK3CD) (Fig. 5a). Transcript expression increased for each of these genes except EGF during the PrF-PF transition ( Table 3, Fig. 5a). In contrast, the non-significantly differentially expressed gene EGF decreased sharply during the PrF-PF transition only to increase sharply during the PF-SF transition (Fig. 5a). (2) TGF-β signalling pathway: growth differentiation factor 9 (GDF-9); BMP and activin membrane bound inhibitor (BAMBI); BMP15; twisted gastrulation BMP signalling modulator 1 (TWSG1); SMAD specific E3 ubiquitin protein ligase 2 (SMURF2); BMP4; bone morphogenetic protein receptor type 1B (BMPR1B); and BMP binding endothelial regulator (BMPER) ( Table 3, Table 1. Functional annotation clusters identified during the PrF-PF transition in domestic cat using DAVID. The top 20 categories grouped by similar GO and KO terms are listed. The full list containing 137 identified clusters is found in Supplementary Data 2, "DAVID PrF-PF". 1 enrichment score, geometric mean of member's P values of the corresponding annotation cluster in -log10 scale of the annotation cluster; 2 number of gene counts, 3 Fold-change, 4 Benjamini-Hochberg value, and 5 false discovery rate (P value adjusted). *internal gene ID defines a unique gene cluster belonging to a gene entry. www.nature.com/scientificreports/ Fig. 5b). Overall, the genes associated with the TGF-β signalling pathway demonstrated the most dynamic transcript expression. The highest "baseMean" value was estimated for GDF-9 however, this gene was not significantly differentially expressed (Table 3). Transcript expression of GDF-9 was in a comparable range throughout each ovarian follicle stage (Fig. 5b). BAMBI, BMP15, TWSG1, BMP4 and BMPR1B showed a significant increase from PrF to PF (Fig. 5b). There was an oppositional pattern observed for SMURF2 and BMPER (Fig. 5b). Transcript counts of BMP4 and BMP15 increases towards SF too but this was not significant (Fig. 5b). (3) ErbB signalling pathway: Erb-B2 receptor tyrosine kinase 2 (ERBB2); SH3 domain containing kinase binding protein 1 (SH3KBP1); and ErbB factor neuregulin 2 (NRG2) (Fig. 5c). A sequential decrease all the way until the SF stage was observed for NRG2. In contrast, expression of ERBB2 and SH3KBP1 transcripts increases from PrF to PF stage (Fig. 5c, Table 3). (4) HIF-1 signalling pathway: phosphofructokinase, liver type (PFKL); hypoxia inducible factor 1 subunit alpha inhibitor (HIF1AN); egl-9 family hypoxia inducible factor 2 (EGLN2); basic www.nature.com/scientificreports/ helix-loop-helix family member e40 (BHLHE40); and furin, paired basic amino acid cleaving enzyme (FURIN) (Fig. 5d). Transcript counts for three of them (PFKL, HIFNAN, EGLN2) increased significantly during PrF-PF transition whereas the other two increased strongly in the PF-SF transition (Fig. 5d, Table 3). (5) Variable transcript expression patterns were found in the MMP group: MMP12, MMP21 and MMP7 expression decreased   www.nature.com/scientificreports/ 63 (CCDC63); armadillo repeat containing 4 (ARMC4); dynein axonemal heavy chain 5 (DNAH5); and sperm associated antigen 1 (SPAG1) (Fig. 5f). All counts were highest in PrFs in comparison to PFs. Values between PF and SF were comparable (Table 3).
qRT-PCR comparative analysis. The genes, BMP15 and histone H1t (HIST1H1T), were compared to the RNA-sequencing data by evaluating relative gene expression trends with qRT-PCR data. Diagnostic plots with RNA-sequencing data of selected genes demonstrated an increase in BMP15 gene expression as ovarian follicles developed from PrFs toward SFs ( Supplementary Fig. S7a). A decrease in HISTH1T gene expression was observed as ovarian follicles progressed from the PrF stage onward (Supplementary Fig. S7a). The negative binomial general linearised model of normalised RNA-sequencing data showed a similar increase in BMP15 gene expression as the ovarian follicles developed ( Supplementary Fig. S7b, above). Applying the same approach to HISTH1T, gene expression was detected in PrF only ( Supplementary Fig. S7b, below). Analysis of BMP15 by qRT-PCR revealed that gene expression was nearly undetectable in PrFs and PFs however, BMP15 gene expression was observed in SFs ( Supplementary Fig. S7b, above right). The Kruskal-Wallis test detected significant changes of gene expression throughout ovarian follicular development (P value = 0.04627), but the subsequent post-hoc test reported only a tendency for higher expression in SF compared to PrF and PF. In contrast, a high relative gene expression of HIST1H1T was found in PrFs which decreased considerably as ovarian follicles progressed toward PF and SF stages, the differences of PrF and PF as well as PrF and SF were confirmed as statistically significant ( Supplementary Fig. S7b, below right).

Discussion
As far as we are aware, this is the first whole ovarian follicle transcriptomic analysis during primordial ovarian follicle (PrF) development and early follicular growth during the primary (PF) to secondary (SF) ovarian follicle transition in a non-rodent animal. In contrast to most gene expression studies performed on adult ovaries, in this study the ovarian follicle is considered as a functional unit. There was no differentiation made between oocytes and the surrounding ovarian follicular cells. This is primarily due to the fact that after isolation from the ovarian cortex, dissection of oocytes from early ovarian follicles is impossible without any severe damage. However, laser dissection microscopy (LCM) is an alternative method which may circumvent such problems and has thus far, been executed in human 20,49 and sheep 16,17 . Nevertheless, the use of whole ovarian follicles as a source for transcriptome analysis is essential to monitor follicular growth in vitro. The initial RNA-sequencing analysis revealed a higher proportion of differentially expressed genes (DEGs) during the activation of dormant ovarian follicles (PrF-PF: 2,226 DEGs) in comparison to the growth phase (PF-SF: 156 DEGs) (Fig. 2). Analysing transcripts of early stages of human oocytes, 223 and 268 genes were detected as significantly expressed in oocytes from PrFs and PFs, respectively 50 . The differences in differential gene expression levels may be due to the use of whole ovarian follicles in comparison to the latter oocyte-specific approach. Another group analysed oocytes and granulosa cells (GCs) of ovine ovarian follicles separately 17 . In this study, less than 200 DEGs in oocytes and GCs during the PrF-PF transition were detected. In contrast to this data, a higher number of gene expression changes during the PF-SF transition were observed-around 500 DEGs in oocytes and 400 DEGs in GCs 17 . So far it is not clear if the discrepancies in the number of DEGs are species-specific or due to the particular analytical approach such as, intact ovarian follicles versus laser dissected ovarian follicles from tissue slices. Additionally, we must consider that in this study the ratio of expressed genes of oocyte origin is decreasing with ovarian follicular growth. An intact PrF contains < 20 follicular cells per oocyte, whereas, this ratio is shifted to 1:50 for PF and 1:200 SFs, respectively 51 . To account for the number of ovarian follicular cells the number of follicles per pool were adapted (PrF n = 180, PF n = 45, and SF n = 9). Therefore, differential gene expression in PF-SF is mirroring more GC expression than oocyte expression, maybe leading to a lower number of detectable DEGs than in PrF-PF. In the future, it will be worthwhile to consider single cell transcriptomics to investigate gene expression levels individually. Figure 5. Transcript expression plots of normalised counts for genes of interest from preantral ovarian follicles in domestic cat (RNA-sequencing data). (a) Genes associated with the PI3K-Akt pathway from left-to-right: TSC complex subunit 2 (TSC2); pyruvate dehydrogenase kinase 4 (PDK4); AKT serine/threonine kinase 2 (AKT2); epidermal growth factor (EGF); regulatory associated protein of mTOR complex 1 (RPTOR); and PI3K subunit delta (PIK3CD); (b) TGF-β pathway from left-to-right: growth differentiation factor 9 (GDF-9); BMP and activin membrane bound inhibitor (BAMBI); bone morphogenetic protein 15 (BMP15); twisted gastrulation BMP signalling modulator 1 (TWSG1); SMAD specific E3 ubiquitin protein ligase 2 (SMURF2); BMP4; bone morphogenetic protein receptor type 1B (BMPR1B); and BMP binding endothelial regulator (BMPER); (c) erythroblastoma (ErbB) pathway: Erb-B2 receptor tyrosine kinase 2 (ERBB2); SH3 domain containing kinase binding protein 1 (SH3KBP1); and ErbB factor neuregulin 2 (NRG2); (d) hypoxia-inducible factor 1 (HIF-1) pathway: phosphofructokinase, liver type (PFKL); hypoxia inducible factor 1 subunit alpha inhibitor (HIF1AN); egl-9 family hypoxia inducible factor 2 (EGLN2); basic helix-loop-helix family member e40 (BHLHE40); and furin, paired basic amino acid cleaving enzyme (FURIN)  www.nature.com/scientificreports/ Table 3. Transcript expression DESeq2 data from preantral ovarian follicles in domestic cat (RNA-sequencing data). The transcript expression data was extracted from DESeq2 analysis and were summarised in order of decreasing "baseMean" values. "PI3K-Akt pathway" genes: TSC2; PDK4; AKT2; EGF; RPTOR; and PIK3CD. "TGF-β pathway" genes: GDF-9; BAMBI; BMP15; TWSG1; SMURF2; BMP4; BMPR1B; BMPER. "ErbB pathway" genes: ERBB2, SH3KBP1, and NRG2. "HIF-1 pathway" genes: PFKL; HIF1AN; EGLN2; BHLHE40; FURIN. "Extracellular matrix" genes: MMP12; MMP2; MMP1; MMP21; MMP7; and MMP9. "Axonemal dynein complex" genes: DNAH7; DRC7; CCDC63; ARMC4; DNAH5; and SPAG1. "baseMean" is the average of the normalised count values, dividing by size factors, taken over all samples; "FC" fold-change is log2FoldChange, the effect size estimate. This value indicates how much the gene or transcript's expression seems to have changed between the contrasts reported on a logarithmic scale to base 2; "padj" is the adjusted P value for multiple testing for the gene or transcript; PrF-PF and PF-SF denote primordial versus primary and primary versus secondary DESeq2 contrasts; and * indicates differential gene expression significance with an adjusted P value of < 0.05 and log2 fold-change ≥ 1. www.nature.com/scientificreports/ After differential gene expression analysis the gene lists were functionally annotated resulting in the identification of over-represented GO and KO terms. As described previously, conserved signalling pathways and biological processes (BPs) that were over-represented were considered. This included the two signalling pathways predominantly studied during ovarian folliculogenesis in other mammals: the PI3K-Akt and TGF-β pathways. Additionally, signalling pathways more ambiguously described though implicated in ovarian folliculogenesis such as, ErbB and HIF-1 pathways were given attention too. The PI3K-Akt signalling pathway within the oocyte is a key regulator of PrF activation and has been described during ovarian folliculogenesis in bovine, human, ovine, and porcine 52 . Currently, it is one of the major non-gonadotrophic growth factor pathways regulating ovarian follicles. As expected, the PI3K-Akt signalling pathway was identified during the PrF-PF transition in domestic cat ( Table 2). The TGF-β signalling pathway was identified in the PrF-PF data under the parent GO term "transmembrane receptor protein serine/threonine kinase signalling pathway" (Table 2). Similarly to the PI3K-Akt pathway, the TGF-β signalling pathway has been studied during ovarian folliculogenesis in bovine, human, ovine, porcine, rodents, and rhesus monkeys 53 . For the domestic cat, the development of ovarian follicles in vitro with the supplementation of PI3K-Akt and TGF-β-associated factors has been investigated 9,54 . The studies focused on the effect of epidermal growth factor (EGF), its receptor (EGFR), and the growth differentiation 9 factor (GDF-9) supplementation either in combination or alone. Signalling via EGF and EGFR up-regulates the PI3K pathway 55 . In goat, EGF stimulates in vitro growth during the PrF-PF transition 56 and in rat, it is implicated in the growth of PrFs toward the SF stage 57 . The TGF-β factor, GDF-9, has been shown to have a similar influence on preantral ovarian follicles in several mammals. In vitro studies have implicated GDF-9 supplementation in PrF activation and ovarian follicle viability in human 58 , goat 59,60 , bovine 60 , and hamster 61 . Additionally, GDF-9 has an over-arching function within the PI3K pathway in rat preantral ovarian follicles 62 . Interestingly, for the domestic cat, the culture of ovarian cortical slices in medium supplemented with EGF and/or GDF-9 have shown that EGF but not GDF-9 improved follicle viability 54 . Medium containing GDF-9 not only had no beneficial influence on ovarian tissues but negligibly impacted ovarian follicle viability 54 . Interestingly, although the transcriptomic data for GDF-9 was determined with a high "baseMean" value it was not significantly differentially expressed (Table 3). Therefore, for the domestic cat it is likely to be more essential in later stages instead. This is the case in the rat where GDF-9 supplementation in vitro promotes ovarian follicular growth only after the PF stage 63,64 . Other functionally annotated factors have also been studied in other models during ovarian folliculogenesis. This includes BMP4 which is also described in rats to function as an ovarian follicle survival factor promoting PrF development 65 and in humans, where BMP4 is implicated in regulating ovulation via remodelling of the ovarian extracellular matrix (ECM) by facilitating in cumulus-oocyte-complex and mural GC separation 66 . A comparative analysis was performed on BMP15 with the qRT-PCR method. The RNA-sequencing transcript counts of BMP15 were shown to increase from the PF stage onward (Supplementary Fig. S7). Similarly, BMP15 mRNA increased in expression with ongoing follicular development ( Supplementary Fig. S7). This supports previous studies in other species such as, human and rat where increased BMP15 expression promotes GC proliferation and theca layer development 67 .

Gene name symbol baseMean
As mentioned, the ErbB and HIF-1 signalling pathways were included as mechanisms of interest. Currently, there are no studies into the role of these pathways during ovarian folliculogenesis in the domestic cat. However, there is evidence in other species of its involvement in ovarian follicle development. Briefly, the ErbB has recently been found to be down-regulated in human oocytes during the PrF-PF transition 20 . Previously, the ErbB factor neuregulin 1 (NRG1) has been investigated in mural GCs and theca cells from rodent periovulatory follicles 68-70 . In the domestic cat, NRG2 transcript expression was observed along with other ErbB-associated factors such as, ERBB2 and SH3KBP1 (Table 3). In contrast to the rodent model, not only was a different variant found in the domestic cat but a sequential decrease of NRG2 expression was observed all the way until the SF stage (Fig. 5c). Investigations into the HIF-1 pathway may also provide an alternative insight into preantral ovarian folliculogenesis in the domestic cat. Some factors identified had either similar or higher "baseMean" values than those identified for the TGF-β pathway (Table 3). For example, PFKL had the highest "baseMean" overall the factors (excluding GDF-9) (Table 3). Interestingly, each HIF-1 signalling factor demonstrated a tendency to increase in transcript expression during ovarian follicle development (Fig. 5d). In other species, increased HIF-1α expression in GCs is implicated in regulating follicular growth and development in rats 71 . Additionally, HIF-1 signalling in bovine primary GCs is suggested to be involved in steroidogenesis and cell proliferation during follicular development 72 . Although not extensively described here it will be interesting to study these pathways in relation to domestic cat ovarian folliculgenesis in the future.
Noteworthy, "ECM-receptor interaction" and "focal adhesion" were identified within the same functional annotation cluster as "PI3K-Akt signalling pathway" (Table 1). During the development of ovarian follicles the ECM undergoes significant compositional remodelling. Upon the initial activation of PrFs, GCs become cuboidal, forming a PF, and the ovarian follicle moves toward the medullar region and away from the cortical region of the ovary 73 . Previously, in vivo imaging windows in mouse ovaries have revealed that collagen fibers can function as migration tracks for infiltrating tumor cells 74 . Additionally, 3D imaging analyses have revealed that the inward movement of ovarian follicles is likely due to the stiffer cortical enclosure as compared with the softer inner medulla 75 . Furthermore, this analysis demonstrated that ovarian follicles are in close contact to each other 75 . Thus, the migration of activated PrFs and developing preantral ovarian follicles may occur through a collective collaboration based on collagen fiber tracking, morphological reactions, and interactions between neighbouring ovarian follicles. Regarding the ECM, degradation occurs via matrix MMPs. This allows for continued enlargement of the developing ovarian follicle. In the ovary, MMP2 and MMP9 expression has been observed in rat, sheep, mouse, rhesus macaque, horse, cow, and human; MMP1 and MMP13 in rabbit, rat, horse, and rhesus macaque; and MMP7 mRNA in macaque pre-ovulatory GCs 76 . In the domestic cat, MMP1, MMP2, MMP3, MMP7, MMP9, and MMP13 mRNA are detectable at every developmental stage, and the abundance and expression of these enzymes was consistently dynamic 76 (Fig. 5e). Additionally, two MMP transcripts, MMP12 and MMP21, not previously described for the domestic cat were identified. A higher level of transcript expression was observed for MMP12 than MMP21 (Fig. 5e). Additionally, the highest "baseMean" value overall was observed for MMP12 (Table 3) therefore, it may be interesting to perform downstream mRNA analysis on this factor in the domestic cat. Subsequently, the changes that occur within the ECM subsequently affect nuclear dynamics through changes in lamina composition, membrane tension, and nuclear pore size which regulates chromatin configuration and gene expression 77 .
Currently, these changes are not well understood for all species in respect to ovarian folliculogenesis. Screening for factors involved in oocyte dormancy with RNA-sequencing revealed that microenvironmental compression elicited a nuclear rotation response which was mediated by a motor protein called dynein in mice 73 . The inhibition of dynein significantly increased the number of growing follicles demonstrating its essential role in follicle dormancy in cultured murine ovaries 73 . Differentially expressed genes of the axonemal dynein complex were identified in the PrF-PF transition only (Supplementary Data S2). The heavy chain regions of the dynein complex contain the motor domain which is capable of producing movement in vitro 78 . In the transcriptomic data, the DNAH5 and DNAH7 heavy chains transcripts were identified along with the regulatory DRC7 and structural factors CCDC63 and ARMC4 ( Fig. 5f and Table 3). Thus, the two heavy chain transcripts, DNAH5 and DNAH7, from the nine major phylogenetically classed dynein heavy chains may be involved in PrF dormancy and/or activation in the domestic cat. This, of course, needs further investigation in the future. Additionally, very high expression of histone HIST1H1T (H1.6) in PrFs was observed ( Supplementary Fig. S7). H1 histones belong to the so-called linker histones, meaning that they are not part of the histone-octomer-centre of nucleosomes but connect to the ends of DNA that coils around such a centre leading to chromatin compaction. Furthermore, H1 also functions through interactions with other proteins that will in turn modify chromatin or take part in DNA-based processes like transcription 79 . Besides versions of H1 (H1.1 to H1.5) which are present in most somatic cells, other subtypes exist 79 , some of which are germ-line variants 80 . HIST1H1T was so far described as a testis-specific variant 81 , but could be detected in tumour cell lines, mouse embryonic stem cells and some normal somatic cells also 82 . To our knowledge, this is the first description of its presence in ovarian follicles, meaning that it could be present in female germ cells. Male H1 germ variants like H1T appear after spermatogonia stop proliferating and before transition proteins are detected and histones are replaced by protamines at late spermatid stages. H1T is the first variant to be expressed in meiotic spermatocytes, the presence in later germ cell stage up to elongated spermatids seems to be species-specific 80 . The dormant PrFs are arrested in meiosis I 83 and primary spermatocytes are in prophase of meiosis I 84 . It could be possible that H1T fulfils a specific function during this cell division phase, for example, on chromatin structure. It is discussed that H1T supports the highly decondensed chromatin stage in early spermatocyte 80 . However, is has been shown that H1T repressed at least ribosomal DNA transcription by condensing chromatin structure 82 .

Conclusion
RNA-sequencing analysis of ovarian follicles in domestic cat contributed to the increasing knowledge on factors and processes which regulate the recruitment and growth of ovarian follicles. Many biological processes were comparable to known data obtained for species such as; human however, species-specific features for the domestic cat may be present. In this study, the analysis focused on signalling factors and pathways along with mechanistic cues associated with the ECM and nuclear dynamics. Overall, the results are relevant to fundamental ovarian follicle developmental biology with an outlook toward developing techniques such as, IVG of ovarian follicles in the future.