Seasonal effects on miRNA and transcriptomic profile of oocytes and follicular cells in buffalo (Bubalus bubalis)

Season clearly influences oocyte competence in buffalo (Bubalus bubalis); however, changes in the oocyte molecular status in relation to season are poorly understood. This study characterizes the microRNA (miRNA) and transcriptomic profiles of oocytes (OOs) and corresponding follicular cells (FCs) from buffalo ovaries collected in the breeding (BS) and non-breeding (NBS) seasons. In the BS, cleavage and blastocyst rates are significantly higher compared to NBS. Thirteen miRNAs and two mRNAs showed differential expression (DE) in FCs between BS and NBS. DE-miRNAs target gene analysis uncovered pathways associated with transforming growth factor β (TGFβ) and circadian clock photoperiod. Oocytes cluster in function of season for their miRNA content, showing 13 DE-miRNAs between BS and NBS. Between the two seasons, 22 differentially expressed genes were also observed. Gene Ontology (GO) analysis of miRNA target genes and differentially expressed genes (DEGs) in OOs highlights pathways related to triglyceride and sterol biosynthesis and storage. Co-expression analysis of miRNAs and mRNAs revealed a positive correlation between miR-296-3p and genes related to metabolism and hormone regulation. In conclusion, season significantly affects female fertility in buffalo and impacts on oocyte transcriptomic of genes related to folliculogenesis and acquisition of oocyte competence.

Water buffalo (Bubalus bubalis) is an important livestock resource for both developing and developed countries. The major factor affecting buffalo farming profitability is reproductive seasonality, resulting in cycles of calving and milk production. Buffalo is a short-day breeder, with increased fertility in response to decreasing day length 1,2 . This photoperiod-dependent seasonality pattern is more pronounced as distance from the equator, together with variations in the light/dark ratio, increases. In Italy, in order to satisfy market demand, out of breeding mating strategy (OBMS), consisting in interrupting sexual promiscuity or the use of artificial insemination (AI) during the breeding season (BS), is commonly utilized 2 . The OBMS improves the distribution of calving throughout the year, but it reduces fertility 3 . Longer post-partum anoestrus periods as well as higher incidence of embryonic mortality are observed in months with increasing daylight length and particularly in mid-winter, which coincides with the transition to seasonal anoestrus at Italian latitudes 1,4 . The embryonic mortality is due to inadequate luteal growth and function, resulting in reduced progesterone secretion 5 . This has a negative impact on embryo growth, associated with alterations in transcriptomic and proteomic profiles of the embryos and chorioamnios/caruncles 6,7 , which ultimately impair embryo attachment to the uterine endometrium. About 10 million reads were sequenced for both OOs and FCs (see Supplementary File S1 for statistics). About 1% and 20% of them were assigned to miRNAs in OOs and FCs, respectively.
In total, 769 miRNAs were identified in at least three samples in all conditions (468 Bos taurus bta-miRNAs, 279 novels, and 22 novels homologous to related species). Among them, 467 were detected in at least three OO samples and 635 in at least three FC samples. Principal Component Analysis (PCA) clearly separates OOs and FCs according to their miRNA content, with 44% of the variance explained by component 1 (Supplementary file S2).
Few of the most expressed miRNAs showed a similar relative abundance in OOs and FCs (bta-miR-10b, bta-miR-148a and bta-miR-26a); on the contrary, expression rate of most miRNAs differed in the two cellular types (bta-miR-21-5p was highly expressed in FCs, whereas bta-miR-423-3p in OOs), (Supplementary file 3). In fact, there was a statistically significant difference in the expression of a high proportion of the miRNAs (n = 413) (False Discovery Rate FDR < 0.05) between OOs and FCs (Supplementary file 4). When the two seasons were considered, the PCA produced a good distinction between oocytes collected in BS and NBS, whereas FCs from BS and NBS could not be clearly distinguished (Fig. 1).
The number of differentially expressed miRNAs (DE-miRNAs, FDR < 0.05) between the two seasons was 13 for both OOs and FCs (Supplementary file 5 and Table 1). A view of the normalized expression of the most representative DE-miRNAs is shown in Fig. 2. Target prediction using human miRNAs homologous to buffalo DE-miRNAs led to the identification of 6,712 and 4,847 genes potentially regulated in OOs and FCs, respectively (P < 0.05). GO analysis using a subset of more significant target genes (n = 136 with P < 0.0005 for OOs, n = 139 with P < 0.001 for FCs) identified pathways related to triglyceride and cholesterol metabolism and transport, and mesoderm and epithelial cell morphogenesis differentiation for OOs, and related to photoperiodism, circadian clock regulation, and transforming growth factor beta signalling for FCs (Table2).
RNASeq. RNA-seq analysis was performed on the same samples used for miRNA profiling to evaluate the gene expression variation between the two cellular types and seasons. Approximately 23.5 ± 4.4 and 54.5 ± 10.5 millions of reads were obtained for OOs and FCs samples with a mapping rate of 93.5% and 92.4%, respectively (Supplementary file 6). A total of 22,013 unique genes present in at least three samples from both cellular types were identified (19, www.nature.com/scientificreports/ considering the relative expression of these genes showed a clear separation between the two cellular types (Supplementary file 7). There seems to be no seasonal effect in the overall transcript abundance for both OOs and FCs (Fig. 3).
The relative expression of mRNAs in the two cellular types was very different, with 14,680 (66.7%) DEGs between OOs and FCs. When DEGs were calculated between the two seasons, 22 mRNAs were found to differ between NBS and BS for OOs, whereas only two DEGs were present in FCs (Supplementary file 8 and Table 3).  www.nature.com/scientificreports/ Although a limited number of DEGs was found to differ in OOs between the two seasons, GO analysis revealed that some of them were related to lipid storage and localization and regulation of interleukin-8 (IL8) production (Table 4). miRNAs and mRNA interaction. In order to evaluate whether miRNAs could potentially regulate the expression of specific genes, the list of genes differentially expressed in oocytes between the two seasons was intersected with the list of DE-miRNA target genes observed in the same experimental condition. Six genes (CCL1, FOLR2, IGF2, HSPA1A, IL1B, CTSK) were found to be potentially regulated by specific DE-miRNAs. Interestingly, among the 8 DE-miRNAs (miR-143, miR-1468, miR-199a-3p, miR-199a-5p, miR-222, miR-25, miR-296-3p, miR-331-5p) targeting 6,712 genes, miR-296-3p targets 4 out of the 6 shared DEGs. In addition, all the miRNAs show a positive correlation with gene expression (Fig. 4).

Discussion
The present study aims to investigate the causes of the decreased oocyte competence during the NBS in buffalo. Many biological processes are required for developmental competence, with the exchange of information between oocyte and follicular environment promoting oocyte maturation 26 . Therefore, the focus of this study was to evaluate differences in the miRNA and transcriptomic profiles of OOs and corresponding FCs from buffalo www.nature.com/scientificreports/ www.nature.com/scientificreports/ ovaries collected in the BS and NBS. To our knowledge, this is the first study reporting together the miRNA and mRNA profiling from pools of low numbers of oocytes and corresponding FCs collected from abattoir-derived ovaries in livestock, as well as the first time that the seasonal effects on miRNA and mRNA profiling of oocytes and FCs are investigated in buffalo. Unfortunately, the amount of RNA obtained from such a limited number of oocytes was not sufficient to perform further experiments to validate our results.
In accordance with previous findings 10,11 , in the present study we observed a reduced oocyte developmental competence in the NBS, as indicated by decreased cleavage and blastocyst rates after in vitro fertilization (IVF), and this was associated to changes in miRNA and transcriptomic profiles both in OOs and FCs. Being miRNAs only one of the small RNA components, as expected those identified in buffalo represent only a fraction of the total small non-coding RNA present in both OOs and FCs 27,28 . The overall miRNA expression pattern was different between OOs obtained in the two seasons. However, FCs from the two seasons cluster together. Considering that different phenotypes of FCs can be observed within the follicle, it is likely that miRNA expression in FCs is mainly driven by cell position and function, thus masking seasonal effects 29 . Recently, Zhang et al. 30 reviewed miRNA profile studies related to ovarian development and function in mammals. Although many studies reported miRNA ovary profiling, only a few investigated miRNA variations in granulosa cells and only one in OOs 30 .
Interestingly, many of the differentially expressed miRNAs in the two seasons for both OOs and FCs are involved in follicular maturation and development regulation. In our study seasonal changes modified the expression of miR-143, miR-25, miR-222 and miR-199a in buffalo OOs. In the mouse ovary, miR-143 is highly expressed and related to oestradiol production and steroidogenesis gene expression 31 . MiR-143 and miR-25 were also shown to promote progesterone release in human ovarian granulosa cells 32 . Furthermore, cyclic variations in the expression of miR-222 and miR-199a were reported in cattle during follicle maturation, with expression increasing until the mid-luteal phase, and decreasing in the late follicular phase in the bovine dominant follicle 13 . Some of the differentially expressed miRNAs identified between the NBS and BS in buffalo FCs (miR-184, miR-2411 and miR-34c) were also reported to exhibit expression modulation during the cycle in cattle. In particular, temporal miRNA expression dynamics were observed for miR-184 in FCs between days 3 and 7 of the bovine oestrous cycle and for miR-2411 and miR-34c between subordinate and dominant follicles during the early luteal www.nature.com/scientificreports/ phase 33 . MiR-34c was shown to exert anti-proliferative and pro-apoptotic effects in porcine granulosa cells by targeting Forkhead box O3a (FoxO3a) 34 .
In addition, the expression of some of the DE-miRNAs detected in our study differs in several ovarian disorders. It was reported that mir-141 and miR-199a are respectively up and down regulated in human ovarian cancer 35 , miR-184 is a potential predictor of recurrence in human ovarian granulosa cell tumours 36 , and miR-486-5p is downregulated in cumulus cells collected from women affected by polycystic ovary syndrome 37 .
Interestingly, GO analysis of the predicted target genes for DE-microRNAs uncovered pathways associated with OOs and FCs physiology. Oocytes collected from the BS and NBS showed DE-miRNAs able to regulate genes for triglyceride and sterol biosynthesis essential for lipid metabolism, which provides a potent source of energy during oocyte maturation 38 . In FCs, the DE-miRNA target genes were related to pathways involved in transformation of growth factor β (TGFβ) and circadian clock photoperiod. TGFβ promotes granulosa cell proliferation regulating the expression of luteinizing hormone receptor (LH-R) [39][40][41] . Altered photoperiod can affect mRNA expression in ovaries 15 , in fact, transcriptome changes occurred between BS and NBS samples.
Considering DEGs between seasons, it is interesting to note that although only two DEGs were found in FCs, many of the DEGs in the OOs are known to be related to oocyte competence. In the NBS, decreased oocyte competence in buffalo was associated to change in the expression of secreted phosphoprotein 1 (SPP1), RUNX family transcription factor 2 (RUNX2) and Cathepsin K (CTSK) in OOs. Both SPP1 and RUNX2 expression was observed to change in oocyte and-granulosa cell complexes at various stages of follicle development in pigs 42 . In addition, variations in the expression of SPP1 were recorded in cumulus cells derived from Table 3. Differentially expressed gene DEGs (false discovery rate (FDR) < 0.05) calculated between the two seasons (non breeding season NBS, breeding season) for oocytes (OOs) and follicular cells (FCs).   43 , and RUNX2 expression was associated with controlled ovarian stimulation outcome in assisted reproductive technology treatment in women 44 . Furthermore, CTSK in cumulus cells was suggested as a predictive marker for oocyte competence in bovine COC 45 .

GENE ID Human and cattle ortholog logFC FDR GENE ID Human and cattle ortholog logFC FDR
In buffalo OOs during the NBS, a decreased expression of heat shock protein family A (Hsp70) member 1A (HSPA1A), known to be related to oocyte survival and apoptosis, was also observed. The HSPA1A plays a critical role through its protective action against apoptosis and its expression is reduced in poor, as compared to competent, ovine COCs 46 .
Another transcript down-regulated in the NBS OOs is interleukin-1 beta (IL-1β). Although IL-1β deficiency in mice prolongs ovarian lifespan 47 ), IL-1β stimulates the growth and sustains maturation in mare 48 and bovine oocytes 49 . In addition, IL-1β was postulated to be involved in different ovulation-associated events such as prostaglandin production and steroidogenesis 50 . Modulation of expression levels was observed in this study also for other genes related to gonadotropic hormone synthesis and metabolism, showing a reduced expression in buffalo OOs during NBS. Apolipoprotein E (APOE) is expressed in cultured ovarian granulosa cells, and is present in human follicular fluid where its relative levels are correlated with serum estrogen concentration 51 . In rats, APOE exerts a role in directing cholesterol during steroidogenesis and regulating follicular estrogenic production 52 . Insulin like growth factor 2 (IGF2) was observed to be expressed in bovine oocytes 53 . The expression of IGF2 is modulated by growth hormone (GH) in in vitro matured Rhesus macaque oocytes 54 . Furthermore, it is known that IGF2, in combination with follicle stimulating hormone (FSH), acts directly on oocyte competence in caprine follicles 55 . Finally, the folate receptor beta (FOLR2), also down-regulated in OOs during the NBS, is a key gene linked to methionine/folate cycles in bovine oocyte 56 , also involved in folate transport in mice oocytes during follicular development 57 .
In our study, a positive correlation between DEGs and DE-miRNA target genes was observed. MiRNAs usually mediate repression of their target mRNAs by inhibiting their translation, therefore reducing the abundance of their products 58 . However, several studies reported a positive miRNA-mRNA regulation with a feed-forward mechanism probably mediated by transcription factors 59 . Notably, among all DE-miRNAs, miR-296-3p is the most correlated with transcripts changing in OOs between BS and NBS. MicroR-296-3p was previously reported to be expressed in ovaries in mice 60 and to repress cell plasticity in different tumour lines 61 , promoting apoptosis in liver 62 and in mammalian pancreatic α cells 63 . Recently, altered expression of bta-miR-296-3p was detected in muscle, kidney, and liver, in bovine foetuses with large offspring syndrome (LOS) 64 . In addition, miR-296 was also observed to be epigenetically regulated as a part of the imprinted Gnas/GNAS clusters 65 .

Conclusion
In conclusion, the reduced oocyte developmental competence recorded during the NBS in buffalo is associated with changes in miRNA and mRNA content in OOs and corresponding FCs. The GO analysis showed overrepresentation of key genes related to lipid and sterol biosynthesis and hormone regulation, crucial for folliculogenesis and acquisition of oocyte competence. These observations might help to explain the seasonal difference in the potential of buffalo oocytes, thus providing the basis for the development of strategies to improve oocyte competence in the NBS. Nevertheless, further efforts are still needed to validate expression modulation of miR-NAs and key genes identified in our study and deeply investigate their role in seasonal reproduction in buffalo. www.nature.com/scientificreports/

Materials and methods
Collection of oocytes and granulosa cells. The study was carried out in Southern Italy (latitude 40.5°-41.5° N and longitude 13.5-15.5) in October, i.e. autumn (BS) and January, i.e. mid-winter (NBS). Buffalo ovaries were collected at a local slaughterhouse (Real Beef s.r.l., Flumeri (AV), Italy under national food hygiene regulations, and transported to the laboratory in physiological saline supplemented with 150 mg/L kanamycin at 30-35 °C within 4 h after slaughter. In order to reduce variability, the ovaries were collected from a homogeneous population of buffaloes, i.e. 134 cyclic multiparous Italian Mediterranean Buffalo cows with a mean weight and age of 552.6. ± 12.1 kg and 5.3 ± 0.4 years, over a total of 10 replicates (5/season). Cyclic ovarian activity was assessed by two clinical examinations carried out 12 days apart before slaughter, to detect the presence of a follicle greater than 1 cm and/or corpus luteum on the ovary. For each day of collection (n = 10), 2-8 mm follicles were aspirated under controlled pressure to collect both OOs and FCs for molecular analyses, while a group of cumulus oocyte complexes (COCs) were in vitro matured, fertilized and cultured up to the blastocyst stage (n = 238 and 234, respectively in the BS and NBS).
Follicular fluid was aspirated using an 18 G needle under vacuum (40)(41)(42)(43)(44)(45)(46)(47)(48)(49)(50) in Falcon tubes and poured into a petri dish for COC recovery. The COCs were evaluated according to morphology and classified according to Di Francesco et al. 10 . Grade A and B COCs, considered suitable for in vitro embryo production (IVEP), were quickly selected from the dish and washed thoroughly in medium H199.
For each replicate, COCs were denuded of their cumulus cells by gentle pipetting and denuded oocytes were washed in phosphate buffer solution (PBS) + 0.1% polyvinyl alcohol (PVA), pooled (20/pool), snap frozen in liquid nitrogen and stored at − 80 °C until RNA isolation.
The follicular fluid was centrifuged at 300×g for 10 min at 4 °C to separate the follicular fluid and the FCs. After centrifugation, the supernatant was centrifuged again at 2000g for 10 min and the pellet containing FCs was snap frozen in liquid nitrogen and stored at − 80 °C until RNA isolation.
In vitro embryo production. Unless otherwise stated, reagents were purchased from Sigma Chemical Company (Milano, Italy). The methods for in vitro maturation (IVM) described below have been reproduced in part from Gasparrini et al. 66 . For each replicate, Grade A and B COCs recovered by follicular aspiration were rinsed in HEPES-buffered TCM199 supplemented with 10% fetal calf serum (FCS) and in vitro matured, fertilized and cultured to the blastocyst stage. Briefly, COCs were allocated to 50 µL drops (10 per drop) of IVM medium, i.e. in TCM199 buffered with 25 mM sodium bicarbonate and supplemented with 10% FCS, 0.2 mM sodium pyruvate, 0.5 µg/mL FSH, 5 µg/mL LH, 1 µg/mL 17 β-estradiol and 50 µg/mL kanamycin, and incubated at 38.5° C for 21 h in a controlled gas atmosphere of 5% CO 2 in humidified air 66 .
The methods for in vitro fertilization (IVF) and culture (IVC) described below have been reproduced from Di Francesco et al. 2012 11 . Frozen straw from a bull previously tested for IVF were thawed at 37 °C for 40 s and sperm was selected by centrifugation (25 min at 300g) on a Percoll discontinuous gradient (45% and 80%). The sperm pellet was re-suspended to a final concentration of 2 × 10 6 mL −1 in the IVF medium, consisting of Tyrode albumin lactate pyruvate 67 supplemented with 0.2 mM penicillamine, 0.1 mM hypotaurine and 0.01 mM heparin. Insemination was performed in 50 µL drops of IVF medium under mineral oil (5 oocytes per drop) at 38.5° C under humidified 5% CO 2 in air. Twenty hours after IVF, putative zygotes were removed from the IVF medium, stripped of cumulus cells by gentle pipetting and allocated to 20 µ@@L drops of IVC medium, i.e. synthetic oviduct fluid (SOF) including essential and non-essential amino acids and 8 mg/mL bovine serum albumin 68 . Culture was carried out under humidified air with 5% CO 2 , 7% O 2 and 88% N 2 at 38.5 °C. On day 5 post-insemination (pi) the cleavage rate was assessed and the embryos transferred into fresh medium for further 2 days of IVC, when blastocyst rates were recorded. RNA isolation. Samples for RNA isolation were obtained from pools (n = 20) of OOs and FCs for both conditions (BS and NBS). The methods described below have been reproduced in part from Lange-Consiglio et al. 69 . Total RNA was isolated by NucleoSpin miRNA kit (Macherey-Nagel, Germany), following the protocol in combination with TRIzol (Invitrogen, Carlsbad, CA, USA) lysis with small and large RNA in one fraction (total RNA). Concentration and quality of RNA were determined by Agilent 2,100 Bioanalyzer (RIN ≥ 6.5 and 7.5 for OOs and FCs, respectively) (Santa Clara, CA, USA). The isolated RNAs were stored at − 80 °C until use. Data analysis. miRNA analysis. Illumina raw sequences were quality checked with FastQC (https ://www. bioin forma tics.babra ham.ac.uk/proje cts/fastq c/) and trimmed with Trimmomatic (version 0.32) 71 , then miR-Deep2 (miRDeep2 (version 2.0.0.5) 72 was used for miRNA detection and discovery. Known miRNAs available at MirBase (https ://www.mirba se.org/) were used to support miRNA identification. In particular, Bos taurus miRNAs were input to support known miRNA detection and miRNAs from related species (sheep, goat and human) were input to support novel miRNA identification. All the identified miRNAs were quantified using the Scientific RepoRtS | (2020) 10:13557 | https://doi.org/10.1038/s41598-020-70546-5 www.nature.com/scientificreports/ miRDeep2 quantifier module. The Bioconductor edgeR package (version 2.4) was used to identify statistically significant differential expression between groups of samples (false discovery rate [FDR] < 0.05) 73 . Predicted miRNA gene targeting of differentially expressed Bos taurus miRNAs (DEmiRNAs) was performed with miR-Walk2.0 74 , using homologous human miRNAs as input identifiers. Target genes were submitted to GO analysis. GO classification of the DEGs was performed according to canonical GO categories, using the Cytoscape (version.3.2.1) plug-in ClueGO (version 2.3.5) which integrates GO and enhances biological interpretation of large lists of genes 75 . MicroRNA cluster analysis was performed with Genesis (version1.8.1) 76 .
In vitro embryo production. Differences in cleavage and blastocyst rates between seasons were analyzed by Chi square test. The level of significance was set at P < 0.05.

Data availability
RNA-Seq data are available in the Sequence Reads Archive (SRA), BioProject accession number, PRJNA599337. Novel miRNA precursors and novel miRNA mature sequences are reported in Supplementary files S10 and S11.