Cumulus Cell Transcriptome after Cumulus-Oocyte Complex Exposure to Nanomolar Cadmium in an In Vitro Animal Model of Prepubertal and Adult Age

Simple Summary Cadmium (Cd), a highly toxic environmental contaminant, negatively affects human and animal fertility in females. In sheep, Cd displays age-dependent bioaccumulation at ovarian level. At environmental nanomolar concentrations, it reduces oocyte fertilization by inducing oxidative stress on cumulus cells (CCs), a cell population of the cumulus-oocyte complex (COC) supporting oocyte growth, functional maturation and fertilization. In this study, the modifications induced by Cd exposure to all the genes expressed in CCs of in vitro matured COCs, recovered from the ovaries of adult and prepubertal sheep, were analysed by RNA sequencing. A set of genes significantly dysregulated upon Cd exposure was identified. Effects of Cd were more relevant in CCs from adult than from prepubertal COCs. Most genes were upregulated while a minority of them were downregulated. Some genes were already known as involved in ovarian activity or Cd-induced effects, whereas others were completely new in these fields. These findings identify in the sheep, an important livestock species with translational relevance in human reproduction, the set of genes controlling oocyte functional competence, altered by Cd. These biomarkers will make it possible to identify oocytes that cannot be fertilized to evaluate whether they are to be discarded or recovered with detoxifying treatments. Abstract Cadmium (Cd), a highly toxic pollutant, impairs oocyte fertilization, through oxidative damage on cumulus cells (CCs). This study analysed the transcriptomic profile of CCs of cumulus-oocyte complexes (COCs) from adult and prepubertal sheep, exposed to Cd nanomolar concentration during in vitro maturation. In both age-groups, CCs of matured oocytes underwent RNA-seq, data analysis and validation. Differentially expressed genes (DEGs) were identified in adult (n = 99 DEGs) and prepubertal (n = 18 DEGs) CCs upon Cd exposure. Transcriptomes of adult CCs clustered separately between Cd-exposed and control samples, whereas prepubertal ones did not as observed by Principal Component Analysis. The transcriptomic signature of Cd-induced CC toxicity was identified by gene annotation and literature search. Genes associated with previous studies on ovarian functions and/or Cd effects were confirmed and new genes were identified, thus implementing the knowledge on their involvement in such processes. Enrichment and validation analysis showed that, in adult CCs, Cd acted as endocrine disruptor on DEGs involved in hormone biosynthesis, cumulus expansion, regulation of cell signalling, growth and differentiation and oocyte maturation, whereas in prepubertal CCs, Cd affected DEGs involved in CC development and viability and CC-oocyte communications. In conclusion, these DEGs could be used as valuable non-invasive biomarkers for oocyte competence.


Introduction
In recent years, a consistent body of studies reported that environmental contamination, due to chemicals of industrial and biological origin, affects female fertility in humans and animals [1,2]. Heavy metals are a major concern for reproductive health, due to their high global annual emission rate [3] and among them, Cadmium (Cd) is one of the most toxic.
In vitro exposure to Cd was reported to negatively affect oocyte maturation, subsequent fertilization and embryo development in different animal species [16][17][18][19][20][21][22][23] with most of the studies performed by analysing the effects of micromolar Cd concentrations. In a previous study in a sheep model, we found that environmental nanomolar Cd concentrations did not compromise oocyte nuclear maturation; rather, it adversely affected oocyte developmental competence. Indeed, exposure to nanomolar Cd concentrations during in vitro maturation (IVM) impaired oocyte in vitro fertilization by inducing oxidative stress, observed as increased ROS (reactive oxygen species) levels and lipid peroxidation, prevailingly on cumulus cells (CCs) rather than on the oocyte [16]. Based on these data, it emerged that nanomolar Cd induces essentially functional effects while the cumulusoocyte complex (COC) morphology remains unaltered. Furthermore, high percentages of oocytes exposed to Cd reach meiotic maturation. These observations prompted us to conduct studies aimed at identifying non-invasive biomarkers of nanomolar Cd exposure in CCs, allowing oocyte preservation for clinical use in IVF or ICSI programs followed by embryo culture and transfer.
Cumulus cells are a follicular granulosa cell subpopulation which represent an ideal cell substrate for the identification of non-invasive biomarkers of oocyte quality. These somatic cells surround growing oocytes, to which they are directly connected through cytoplasmic protrusions, and can be isolated from the COC without compromising oocyte viability. CCs support oocyte growth, maturation and acquisition of developmental competence [24][25][26][27][28]. Indeed, oocyte nuclear and cytoplasmic maturation depends on rapid transcriptional events, governed by paracrine and autocrine signalling before ovulation in which CCs play significant roles [28][29][30][31]. Once matured, the metaphase II (MII) oocyte is less transcriptionally active and relies on stored mRNA transcripts, acquired throughout maturation, to undergo successful fertilization and early embryo development until embryonic genome activation (EGA) [32][33][34]. Moreover, other than cytoplasmic stored transcripts, CCs provide additional transcripts to the oocyte by active transportation through trans-zonal projections [35,36]. This is why investigating transcriptional CC activity is a key strategy for detecting biomarkers of oocyte quality, improving the outcomes of assisted reproductive programs.
Studies on CCs transcriptomic profile have started in recent years and have developed rapidly in human as well as in some animal species, initially using microarrays [37] and then by RNA-seq [27]. In reproductive toxicology [38], this technology can open the way to the identification of new biomarkers of oocyte functional damage induced by specific contaminants. The aim of the present study was to perform CC transcriptome analysis for the identification of non-invasive biomarkers of oocyte functional impairment due to in vitro exposure to Cd at environmental nanomolar concentration. The IVM of sheep oocytes was used as a large animal in vitro model, with oocytes recovered from adult sheep or prepubertal lambs, representing those with high or low developmental competence, respectively.

Chemicals
All chemicals for in vitro cultures and analyses were purchased from Sigma-Aldrich (Milan, Italy) unless otherwise indicated.

Collection of Ovaries
Ovaries were collected from prepubertal lambs (under 2 months old) and adult sheep (3)(4)(5) years old) at a local slaughterhouse (Fin. Sud Import s.r.l.; Conversano, Bari). All animals were subjected to routine veterinary inspection in accordance with the specific health requirements stated in Council Directive 89/556/ECC and subsequent modifications. Ovaries were transported to the laboratory at room temperature within 2 to 4 h from slaughter.

In Vitro Maturation (IVM)
For the retrieval of cumulus-oocyte complexes (COCs), ovaries were processed differently in relation to the donor age. Ovaries of prepubertal lambs underwent the slicing procedure [39], whereas ovaries from adult sheep underwent follicular fluid aspiration from large developing follicles using a 18G needle. In both procedures, the follicular contents were released in sterile Petri dishes containing phosphate-buffered saline (PBS) supplemented with 1 mg/mL heparin. Only COCs with several intact cumulus cells layers and homogenous cytoplasm were selected for culture. In vitro maturation was performed as previously reported [40]. Briefly, IVM culture medium was composed by TCM-199 with Earle's salts, buffered with 5.87 mM 4-(2-hydroxyethyl)-1-piperazine ethane sulfonic acid (HEPES) and 33.09 mM sodium bicarbonate and supplemented with 0.1 g/L L-glutamine, 2.27 mM sodium pyruvate, calcium-l-lactate pentahydrate (1.62 mM calcium and 3.9 mM Lactate, 50 µg/mL gentamicin, 20% (v/v) Foetal Calf Serum (FCS), 10 µg/mL ovine follicle stimulating hormone (FSH) and 20 µg/mL ovine luteinizing hormone (LH), and 1 µg/mL 17 beta oestradiol. Before IVM, collected COCs were washed three times in TCM-199 with Hank's salts (Gibco ® , Life Technologies, Paisley, UK) supplemented with 10% FCS. COCs were individually cultured in single 10 µL microdrops of IVM medium placed in 60 mm petri dishes covered with pre-equilibrated paraffin oil. This individual IVM culture system allows separation of CCs from MII oocytes from those of immature ones. IVM culture was performed for 24 h at 38.5 • C under 5% CO 2 in air. COCs were exposed to Cd at 100 or 0 (controls) nM. This concentration was chosen on the basis of previous studies reporting that 100 nM CdCl 2 significantly reduced the fertilization rates of oocytes from prepubertal and adult sheep [16]. Cadmium working solution was prepared on the day of use, starting from a 10 mM CdCl 2 stock solution in distilled water.

Cumulus Cells Isolation from Matured Oocytes
After IVM, individual COCs underwent CC removal and oocyte meiotic stage assessment by polar body (PB) visualization under a Nikon SMZ-1500 stereomicroscope. For each experimental condition (100 nM CdCl 2 and controls), CCs isolated from COCs with matured oocytes, showing the 1st PB extruded, were pooled in groups of 20-25 cumuli/group, collected in RNAse-free tubes and washed twice in cold PBS by centrifugation at 300× g for 1 min. After supernatant removal, CC pellets were snap-frozen in liquid nitrogen and stored at −80 • C until molecular analyses. The meiotic stage of oocytes which did not show the 1st PB extruded was assessed by nuclear chromatin evaluation.

Oocyte Nuclear Chromatin Evaluation by Epifluorescence Microscopy
To evaluate nuclear chromatin, oocytes were fixed in 3.8% formaldehyde solution in PBS, stained with 2.5 µg/mL Hoechst 33258 in 3:1 (v/v) glycerol/PBS, mounted on microscope slides covered with cover slips, sealed with nail polish and kept at 4 • C in the dark until observation. Oocytes were evaluated under an epifluorescence microscope (Nikon Eclipse 600, Nikon Instruments, Firenze, Italy; 400× magnification) equipped with a 346 nm excitation/460 nm emission filter, as germinal vesicle (GV), metaphase to telophase I (MI to TI), MII with 1st PB extruded and abnormal [39,41].

RNA-Seq and Data Analysis
For transcriptomic analysis, for each experimental condition (prepubertal Cd-treated, prepubertal controls, adult Cd-treated and adult controls), 5 groups of CCs from matured COCs were used. In each group, CCs isolated from 20-25 matured oocytes were pooled and processed for RNA extraction and CCs from equal numbers of treated and control COCs were analysed. Total RNA was extracted and purified with the mirVana kit (Thermo Fisher Scientific, Waltham, MA, USA) following the manufacturer's protocol. During extraction procedures, RNAs were treated with DNase to exclude any contamination of genomic DNA. Total RNA integrity was evaluated using the 2100-Bioanalyzer (Agilent Technologies, Santa Clara, California, USA) with the RNA PicoLab Chip (Agilent Technologies) and RNA concentration and purity were evaluated using NanoDrop 2000c (Thermo Fisher Scientific). RNA-seq libraries were prepared using the Illumina's TruSeq Stranded Total RNA Sample Preparation Kit (Illumina, San Diego, CA, USA), according to the manufacturer's protocol. Sequencing was performed on Illumina NextSeq500 platform, generating 100 bp pairedend reads. Raw reads in FASTQ format were quality checked using the FastQC program (http://www.bioinformatics.babraham.ac.uk/projects/fastqc, accessed on 12 December 2022) and adaptors as well as low-quality regions (phred cutoff of 25) were trimmed using Trim_Galore. Cleaned reads were aligned onto the reference genome of the ovine species (Oar_v3.1 from UCSC) using STAR (version 020201) [42] with default parameters. Read counts per gene were performed by FeatureCounts (version 1.6.0). Normalized counts of 1000 most variable genes were used to perform a PCA analysis. Differentially expressed genes (DEGs) between Cd-treated and control samples were detected by using the DESeq2 R package [43], selecting only genes with a p value < 0.05 and |log2fc| > 1.2, and displayed as heatmaps (using the zscore (z = (X − µ)/σ) conversion).

Gene Network Analysis
The GeneMANIA (http://www.genemania.org, accessed on 12 December 2022); [44] prediction server was used to analyse the DEG functions and to find neighbouring genes associated with DEGs by constructing the gene network.

RNA-Seq Data Validation by Quantitative RT-PCR
Total RNA was isolated with RNeasy Micro Kit (Qiagen, Hilden, Germany) following manufacturer's instructions. During the procedure, RNA was treated with DNase I to exclude any potential genomic DNA contamination. The isolated RNA was measured using a Nanodrop spectrophotometer (Thermofisher) and used for reverse transcriptionpolymerase chain reaction (RT-PCR). The High-Capacity Complementary DNA (cDNA) Reverse Transcription kit (Life Technologies) was used to convert RNA to cDNA. Each RNA sample (1 µg) was added to 2 µL 10× RT buffer, 0.8 µL 25× dNTP mix, 2 µL RT random primers, 1 µL M-MLV RT, 1 µL RNase inhibitor and nuclease-free H 2 O for a total volume of 20 µL and then mixed gently and centrifuged briefly. Reaction tubes were incubated at 10 • C for 10 min, then at 37 • C for 120 min, and finally at 85 • C for 5 min. The relative quantification of the transcripts was carried out by Real-Time RT-PCR with the StepOne instrument (Applied Biosystems, Foster City, CA, USA). Specific ovine cD-NAs were amplified by PCR using the primers shown in Table 1. PCR was performed in 20 µL reaction volume containing: 10 µL PowerUp SYBR Green PCR Master Mix (Applied Biosystems, 2×), 200 nM of each primer, 2 µL of diluted cDNA (1:10) and nuclease-free water up to 20 µL. Cycling parameters were: 95 • C for 2 min, 40 cycles of denaturation at 94 • C (45 s), annealing at 60 • C (45 s) and extension at 72 • C (45 s), final extension at 72 • C for 5 min. The analysis was carried out in triplicate. Data were collected by using the StepOne Software (Applied Biosystem) and relative quantification was performed by using a comparative method after determining the Ct (threshold cycle) values for the reference endogenous control (beta actin) and the target gene in each sample sets, according to the 2 −∆∆Ct method. Changes in mRNA expression levels were calculated after normalization to Beta actin. The program calculates the ∆Ct and the ∆∆Ct with the formulas below: ∆Ct = Ct_Mean (beta actin) − Ct_Mean (target gene); ∆∆Ct = ∆Ct − ∆Ct_Mean, so that the gene expression level = 2 −∆∆Ct . Changes in gene expression were reported as percentage changes relative to controls.

Statistical Analysis
Oocyte nuclear maturation rates were compared between treated and control groups by the Chi-square test. The evaluation of differences in gene expression between control and treated cells was performed by Student's t-test. Differences with p < 0.05 were considered as statistically significant.

Gene Annotation and Literature Search
To determine the biological significance of our bioinformatic findings, DEGs were cross referenced with available datasets. DEGs were reviewed using the GeneCards database (http://www.genecards.org/, accessed on 12 December 2022) to retrieve general information on their structure and function and to correlate our bioinformatic findings with hallmark physiological and pathological processes in the ovary, when available (Entrez gene summary), and to search for information on gene expression at the mRNA (GTEx and Illumina BodyMap) and protein (Moped and ProteomicsDB) levels in the human normal ovary. Then, the PubMed literature database was used to search for previous studies assessing the expression and functional role of each gene in the ovary of human or other animal species, by associating the gene symbol with the following key words called, alternatively: cumulus cells, granulosa cells, oocyte, ovary.

Nanomolar Cd Does Not Affect Meiotic Progression of Sheep Oocytes
A total of 742 COCs were cultured and analysed, 401 of which were from adult sheep and 341 from prepubertal lambs. In both age groups, exposure to nanomolar Cd during IVM did not affect oocyte meiotic progression and maturation, as no significant differences were identified in the percentages of oocytes found at the examined meiotic stages. Furthermore, no differences were observed in the MII rates of adult versus prepubertal sheep, also considering that prepubertal COCs used in this study were recovered from small antral follicles but displayed the oocyte diameter, number of CC layers and the ability to respond to gonadotropin in vitro stimulation, similarly to their adult counterparts. Overall, in all sample types, oocyte maturation rate reached values around 90% (87.6-95.0) ( Table 2). Table 2. Effects of in vitro exposure to nanomolar Cd during IVM on meiotic progression and maturation of oocytes from adult sheep and prepubertal lambs.

PCA Shows Relevant Cd-Induced Effects on Transcriptomic Profile of CCs from Adult Sheep
Transcriptomes of CCs isolated from oocytes selected after IVM as belonging to matured MII oocytes of adult sheep and prepubertal lambs were sequenced employing the Illumina technology on the NextSeq500 platform. The principal component analysis (PCA) based on gene expression data was performed by grouping sample categories according to donor age (adult versus prepubertal) and culture conditions (Cd-exposed versus controls). In adult samples, Cd exposure induced clear transcriptome separation along the x-axis (main component) ( Figure 1A). This effect was not induced in CCs from prepubertal samples, in which an overlapping gene expression pattern between Cd-exposed and control samples was found ( Figure 1B).

Cadmium Induces Higher Number of DEGs in CCs from Adult versus Prepubertal Sheep
From CCs isolated from matured oocytes of adult sheep, on average, 23.5 million pairs of stranded reads per sample were obtained. About 83% of these reads were uniquely aligned to the reference genome (Ovis_aries.Oar_V3.1) and assigned to known mRNA transcripts. A total of 17,606 gene loci out of 20,028 Ensembl annotations were

Cadmium Induces Higher Number of DEGs in CCs from Adult versus Prepubertal Sheep
From CCs isolated from matured oocytes of adult sheep, on average, 23.5 million pairs of stranded reads per sample were obtained. About 83% of these reads were uniquely aligned to the reference genome (Ovis_aries.Oar_V3.1) and assigned to known mRNA transcripts. A total of 17,606 gene loci out of 20,028 Ensembl annotations were found expressed in CCs cells. By using unique and concordant RNA-seq reads, ninety-nine genes appeared differentially expressed (DEG) upon Cd exposure and are shown in the heatmap of Figure 2. Of these, 73 genes were upregulated (Table S1) whereas 26 genes were downregulated (Table S2) after in vitro exposure to Cd. Heatmaps of DEGs in Cd-exposed versus control CCs from adult sheep. Blue to red colouration denotes low to high expression levels, respectively. Abbreviations: Cd, cadmium; CTRL, controls.
RNA-seq of CCs isolated from matured oocytes of prepubertal lambs produced, on average, 64 million pairs of stranded reads per sample. After the removal of low-quality regions, about 80% of cleaned reads were uniquely aligned to the sheep reference genome and assigned to known gene loci. On the whole, 12,201 genes over the 20,028 known annotations were found expressed. Differential gene expression analysis revealed a lower number of DEGs (n = 18) upon Cd exposure in CCs of oocytes from prepubertal lambs compared to adults, as shown in the heatmap of Figure 3. Compared to controls, 17 genes appeared upregulated and only one was downregulated after Cd exposure (Table S3). To RNA-seq of CCs isolated from matured oocytes of prepubertal lambs produced, on average, 64 million pairs of stranded reads per sample. After the removal of low-quality regions, about 80% of cleaned reads were uniquely aligned to the sheep reference genome and assigned to known gene loci. On the whole, 12,201 genes over the 20,028 known annota-tions were found expressed. Differential gene expression analysis revealed a lower number of DEGs (n = 18) upon Cd exposure in CCs of oocytes from prepubertal lambs compared to adults, as shown in the heatmap of Figure 3. Compared to controls, 17 genes appeared upregulated and only one was downregulated after Cd exposure (Table S3). To better characterize the functional role of DEG genes in both age groups, our attention was focused on protein-coding genes having an orthologue in humans through the GeneCards database. Among adult upregulated genes, 46 were found in the GeneCards database. Of them, six DEGs met all three search criteria (IGFBP2, ITGAX, NOS2, PLA2G4D, SLC27A3 and YBX2) as they were found to be expressed in the human normal ovary, both at mRNA and protein level, and were found to be associated with PUBMED literature data on expression/function at ovarian level. As well, 17 DEGs (APOA2, CCNE2, CDC6, CENPK, CXCL14, CYP19A1, DDIT4L, GDF3, IHH, IL15, INHBE, LOXL1, MEI1, MT1A, MYO5C, NLRP14 and RSAD2) met two of the three criteria as they were found to be expressed in human normal ovary at mRNA level and were associated with protein data from Gene-Cards (Moped or ProteomicsDB) or with RNA/protein data from PUBMED. The remaining 23 DEGs (ABCC8, ACCSL, ADGRF4, ATAD5, CD274, CRISPLD1, DLGAP1, DPP10, EF21, E2F7, ETNPPL, FHAD1, KCNE3, PAX6, PDEA4B, PEAR1, SERTAD4, SLC3A1, SLC6A15, SLC24A5, SPTLC3, RDH5, TRAM1L1) displayed only mRNA expression in human normal ovary in the GeneCards database (Table S1). Among adult upregulated genes, 46 were found in the GeneCards database. Of them, six DEGs met all three search criteria (IGFBP2, ITGAX, NOS2, PLA2G4D, SLC27A3 and YBX2) as they were found to be expressed in the human normal ovary, both at mRNA and protein level, and were found to be associated with PUBMED literature data on expression/function at ovarian level. As well, 17 DEGs (APOA2, CCNE2, CDC6, CENPK, CXCL14, CYP19A1, DDIT4L, GDF3, IHH, IL15, INHBE, LOXL1, MEI1, MT1A, MYO5C, NLRP14 and RSAD2) met two of the three criteria as they were found to be expressed in human normal ovary at mRNA level and were associated with protein data from GeneCards (Moped or ProteomicsDB) or with RNA/protein data from PUBMED. The remaining 23 DEGs (ABCC8, ACCSL, ADGRF4, ATAD5, CD274, CRISPLD1, DLGAP1, DPP10, EF21,  E2F7, ETNPPL, FHAD1, KCNE3, PAX6, PDEA4B, PEAR1, SERTAD4, SLC3A1, SLC6A15,  SLC24A5, SPTLC3, RDH5, TRAM1L1) displayed only mRNA expression in human normal ovary in the GeneCards database (Table S1).
Among prepubertal DEGs, 14 genes were found in the GeneCards database and all of them were upregulated. Of them, three DEGs (CULLIN 2, DSG2 and SLC30A2) met all three search criteria (mRNA and protein expression in normal human ovary on GeneCards and PUBMED literature) and four DEGs (ASTL, BMP15, NLRP5, WEE2) were included in two of the three criteria, being expressed in human normal ovary at mRNA level and associated with protein data from GeneCards (Moped or ProteomicsDB) or with RNA/protein data from PUBMED literature. The remaining seven DEGs (CENPU, GABRA3, HSPA6, MICAL2, MT1A, MT2A, NACHT) displayed only mRNA expression in human normal ovary (Table S3).
Some adult upregulated genes known as being involved in various pathways of CC expansion and oocyte maturation, such as CYP19A1 (aromatase), IGFBP2, (insulin-like growth factor binding protein 2), NOS2 (nitric oxide inducible synthase) and IHH (Indian hedgehog signalling molecule) coding for proteins involved in the regulation of metabolic processes, cell signalling, growth and differentiation, and MT1A (metallothionein 1A), involved in ion transport, were chosen for subsequent enrichment and validation analysis.
In addition, some prepubertal overexpressed genes, known to be involved in various pathways of CC development and viability and CC-oocyte communications, such as SLC30A2 (Solute Carrier family 30 member 2) coding for a protein involved in the regulation of metal ions transport; HSPA6 (Heat Shock Protein family A member 6, Hsp70) having a role in chaperone-mediated protein folding; DSG2 (desmoglein 2) involved in cell-to-cell contact and BMP15 (Bone Morphogenetic Protein 15) involved in the regulation of cumulusoocyte interaction crucial for fertilization, together with MT1A and MT2A (metallothionein 1A and 2A) involved in heavy metal detoxification, were used in subsequent enrichment and validation analysis.

Gene Networks Identifies Different Cd-Induced Pathways in Adult versus Prepubertal CCs
Different gene networks were identified in CCs of adult and prepubertal lambs, indicating that different pathways were activated upon Cd exposure in the two age groups. The gene network of adult upregulated genes CYP19A1, IHH, IGFBP2, MT1A and NOS2 is shown in Figure 4. The network revealed complex interactions among the target genes and between them and other genes (Table S4). Considering the functions, some of these genes were reported to be associated with hormone and steroid hormone biosynthetic processes (CYP19A1, CYP17A1, HSD17B1 and HSD17B3), insulin-like growth factor receptor signalling pathway (IGF1 and IGFBP2), detoxification of inorganic compound (MT1A, MT1B, MT1E, MT1G, MT1H), cell-cell adhesion (IGF1, IGF2, IGFBP2, IHH and SHH) and embryo development (IGF2, IHH, TGFB1 and SHH). Considering the network, it can be seen that: (1) CYP19A1 has physical interactions with CYP17A1, HDS17B1 and HSD17B3 and it is coexpressed and/or has genetic interaction with several proteins of this network; (2) IGFBP2 has physical interactions and pathway interactions with IGF1 and IGF2; (3) NOS2 has various interactions, such as with TGFB1 (genetic interaction), GUCY1A2, GUCY1A1, GUCY1B1 and UCHL5 (physical interactions and pathways), ACTN4 (predicted); (4) IHH has different interactions, such as with PTCH, GAS1 and HHIP (pathways), DHH and SHH (predicted); (5) MT1A shows co-expression with four further metallothioneins. The gene network of prepubertal upregulated genes BMP15, DSG2, HSPA6, MT1A and SLC30A2 is shown in Figure 5. As for adult samples, in prepubertal samples, the network revealed the same complex interactions among the target genes and between them and other genes (Table S4). Considering the functions, some of these genes were reported to be associated with chaperone-mediated protein folding (HSPA6, HSPA8, HPSA1A,  The gene network of prepubertal upregulated genes BMP15, DSG2, HSPA6, MT1A and SLC30A2 is shown in Figure 5. As for adult samples, in prepubertal samples, the network revealed the same complex interactions among the target genes and between them and other genes (Table S4). Considering the functions, some of these genes were reported to be associated with chaperone-mediated protein folding (HSPA6, HSPA8, HPSA1A, STUB1 and DNAJB4), detoxification of inorganic compounds and stress response to metal ion

Cd-Induced CC Expression Pattern Is Consistent with the Sequencing Results
The validation of RNA-seq results was performed by quantitative RT-PCR using three independent biological replicates for each experimental condition (adult Cd-treated, adult controls, prepubertal Cd-treated, prepubertal controls). In adult samples, significantly increased expression of CYP19A1 (p < 0.01), NOS2 (p < 0.05), IGFBP2 (p < 0.05), MT1A (p < 0.05) and IHH (p < 0.05) was found in Cd-exposed CCs compared with controls ( Figure 6, panels A-E). The expression pattern of all examined genes was in agreement with the DEGs profile. As well, in prepubertal samples, significantly increased expression of MT1A, DSG2 and BMP15 was found in Cd-exposed CCs compared with controls, in accordance with transcriptomic profile obtained after RNA-seq (p < 0.05; Figure 7, panels A-D). No significant differences were found for SLC30A2.

Cd-Induced CC Expression Pattern Is Consistent with the Sequencing Results
The validation of RNA-seq results was performed by quantitative RT-PCR using three independent biological replicates for each experimental condition (adult Cd-treated, adult controls, prepubertal Cd-treated, prepubertal controls). In adult samples, significantly increased expression of CYP19A1 (p < 0.01), NOS2 (p < 0.05), IGFBP2 (p < 0.05), MT1A (p < 0.05) and IHH (p < 0.05) was found in Cd-exposed CCs compared with controls ( Figure 6, panels A-E). The expression pattern of all examined genes was in agreement with the DEGs profile. As well, in prepubertal samples, significantly increased expression of MT1A, DSG2 and BMP15 was found in Cd-exposed CCs compared with controls, in accordance with transcriptomic profile obtained after RNA-seq (p < 0.05; Figure 7, panels A-D). No significant differences were found for SLC30A2.

Discussion
This study analysed effects of in vitro exposure of COCs of the sheep model to environmental nanomolar Cd levels on transcriptomic profile of CCs in order to identify noninvasive biomarkers of Cd-induced oocyte dysfunction. To the best of our knowledge, this is the first study reporting the effects of this specific endocrine disruptor on CC transcriptome of a mammalian animal species. Previous studies analysed its effects on ovarian transcriptome in invertebrate species [45][46][47]. Methodological strengths of this reproductive toxicology study include: (1) the use of COCs of two donor age-groups, adult and prepubertal, recovered from two phases of female reproductive life in which oocytes are characterized by higher and lower developmental competence for adult and prepubertal

Discussion
This study analysed effects of in vitro exposure of COCs of the sheep model to environmental nanomolar Cd levels on transcriptomic profile of CCs in order to identify noninvasive biomarkers of Cd-induced oocyte dysfunction. To the best of our knowledge, this is the first study reporting the effects of this specific endocrine disruptor on CC transcriptome of a mammalian animal species. Previous studies analysed its effects on ovarian transcriptome in invertebrate species [45][46][47]. Methodological strengths of this reproductive toxicology study include: (1) the use of COCs of two donor age-groups, adult and prepubertal, recovered from two phases of female reproductive life in which oocytes are characterized by higher and lower developmental competence for adult and prepubertal respectively and, possibly, by different ability to respond to Cd-induced stress; (2) the exposure of COCs at environmental nanomolar Cd concentration, chosen within the range

Discussion
This study analysed effects of in vitro exposure of COCs of the sheep model to environmental nanomolar Cd levels on transcriptomic profile of CCs in order to identify non-invasive biomarkers of Cd-induced oocyte dysfunction. To the best of our knowledge, this is the first study reporting the effects of this specific endocrine disruptor on CC transcriptome of a mammalian animal species. Previous studies analysed its effects on ovarian transcriptome in invertebrate species [45][46][47]. Methodological strengths of this reproductive toxicology study include: (1) the use of COCs of two donor age-groups, adult and prepubertal, recovered from two phases of female reproductive life in which oocytes are characterized by higher and lower developmental competence for adult and prepubertal respectively and, possibly, by different ability to respond to Cd-induced stress; (2) the exposure of COCs at environmental nanomolar Cd concentration, chosen within the range detected in ovarian tissues of prepubertal (0-2.7 ng/g ∼ = 0-24.02 nM) and adult (7.89-24.39 ng/g ∼ = 70.19-216.99 nM) sheep [16] and (3) the choice of exploring the transcriptomic dynamics in CCs which are considered as valuable non-invasive markers for MII oocyte quality [24][25][26][27][28][29][30][31][32][33][34][35][36].
Regarding IVM results, our data are in agreement and confirm those reported in our previous study in which 100 nM Cd did not affect oocyte maturation rate but decreased fertilization [16]. Similarly, in another study [22], exposure of bovine oocytes, during IVM, to 200 nM Cd did not affect oocyte nuclear maturation rate but reduced the cleavage and blastocyst rate compared with controls, indicating functional damage to the matured oocyte. As well, in porcine oocytes, exposure to 400 nM Cd did not affect oocyte maturation rate [23]. To the best of our knowledge, all other studies published to date that have found Cd deleterious effects of Cd on oocyte maturation have examined the effects of higher micromolar Cd concentrations [17][18][19][20][21][22][23]48]. In humans, in agreement with these studies, negative association of follicular fluid Cd nanomolar environmental levels with the probabilities of pregnancy and live birth was found [15].
Based on this large body of previous studies reporting that exposure to nanomolar Cd does not affect maturation, but rather fertilization, we decided to perform a transcriptomic analysis for exploring the molecular signature underlying Cd-induced impaired fertilization. At first glance it emerged that, from the quantitative point of view, in vitro exposure to environmental nanomolar Cd did not induce an excessively significant response in the number of DEGs involved. This is a positive finding considering that this environmental pollutant significantly affects animal and human fertility (See Introduction for References). Indeed, in comparison with recent, even non toxicological, transcriptomic studies on CCs by RNA-seq [27,49], the involvement of a hundred DEGs in adult and 18 DEGs in prepubertal CCs could be considered as indicative of a mild molecular response. However, the observed response should not be underestimated considering the variety and types of altered gene functions playing important roles in the development of oocyte ability to be fertilized. Some of these genes are confirmed by previous studies whereas others are found to be entirely new, as they were never previously associated with COC development and functions. Thus, they represent novel findings from this study that allow knowledge advancement and implementation of the available literature on molecular mechanisms underlying the acquisition of synchronized COC developmental competence.
Moreover, the different intensity of response between the two analysed age groups (higher number of DEGs in adult and lower number in prepubertal; see heatmaps and PCA analysis) led us to consider that these altered functions could be responsible for a full-blown condition in adults (Cd-induced infertility) whereas they could lead to a predisposition in prepubertal subjects (susceptibility to Cd-induced toxicity) in the used animal model. Furthermore, we can observe that basically the response has been mainly of upregulation both in adults and in prepubertal CCs. In our opinion, this could mean that the cells mainly tried to activate defence mechanisms (upregulated genes) and that a true inhibitory effect, at the tested Cd concentration, occurred on a lower number of genes (downregulated genes). Therefore, it could be hypothesized that exposure to this Cd concentration does not definitively damage the COC. Rather, its effects could be reversible, and detoxification strategies could be tested, as reported in a previous study from our group [48]. On the other hand, it cannot be excluded that some of upregulated genes could play inhibitory functions.
Among 46 adult upregulated genes found in the GeneCards database, 17 DEGs were previously associated with studies on cumulus/granulosa/oocyte functions reported in PUBMED. They include regulators of the cell cycle (CCNE2, CDC6); inflammation processes crucial for ovulation (CXCL14, IL15, ITGAX); steroid synthesis (CYP19A1); follicle/COC growth and differentiation (GDF3, IGFBP2, IHH, PLA2G4D); meiotic spindle organization (MEI1); periovulatory processes (MT1A, NLRP14, RSAD2); energetic homeostasis (SLC27A3); oxidative stress and apoptosis (NOS2); and stability and/or translation of germ cell mRNAs (YBX2). Given that some of these genes code for proteins that perform two or more of the aforementioned functions, the results of the present study confirm their expression and functional role in CCs of adult COCs. Moreover, the fact that these genes were overexpressed upon Cd exposure led us to hypothesize that they code for functions useful to prevent COC degeneration, thus preserving it for fertilization. Additional 6 DEGs were reported in GeneCards as expressed both at RNA and protein level in the ovary but, to our knowledge, no studies are available on PUBMED on their functional role in the COC. In light of their ovarian expression, their possible role and involvement in COC preparation to the fertilization process can be hypothesized as they are reported to code for several functions, such as: cholesterol transport (APOA2), component of the centromere (CENPK), signal transduction regulation (DDIT4L), cell proliferation, apoptosis, immune response and hormone secretion regulation (INHBE), extracellular matrix remodelling (LOXL1), actin filament organization and vesicle transport (MYO5C). The remaining 23 upregulated genes have not been annotated previously in studies on oocyte maturation and fertilization (Table S1 ).
Among 21 adult downregulated genes found in the GeneCards database, seven DEGs were associated with previous studies on cumulus/granulosa/oocyte functions reported in PUBMED. They include regulators of the cell cycle (MYCN), apoptosis (CIDEC), ovarian inflammatory processes (LY6G6C, MCP1), G protein signalling (RGS4), sperm-egg fusion (TSPAN18) and post-transcriptional modifications in oocyte target genes (U6). Additionally, six DEGs were reported in GeneCards as expressed both at RNA and protein level in the ovary but, to our knowledge, not in previous PUBMED studies. Based on their ovarian expression, it could be hypothesized that they play a role in COC functions as they are reported to code for: cell cycle progression, signal transduction and apoptosis (CORO1A, PSMA8), proteolysis (PRSS50), thioesterase activity (THEM6), cholesterol homeostasis (TTC39B), metal ion binding activity (ZMYND12). The remaining 8 DEGs (CNTN3, CPNE5,  ISL2, RASGRP3, SLC12A5, TARS3, TTC9, TXK), displaying only mRNA expression in human normal ovary, have never been annotated in previous studies on COC competence (Table S2 [ (Table S3 [ 28,79,).
Overall, the 23 upregulated and eight downregulated genes from adult CCs and the seven upregulated genes from prepubertal CCs, not yet annotated in previous study on COC functions, represent a novel group of genes captured by this study design implicated in its developmental competence. These DEGs should be further explored as they have never been previously reported in such processes.
DEGs used for enrichment and or validation analysis were chosen among upregulated ones and having known role in ovarian function or known response to Cd. They have been analysed to identify their possible role in the modulation of the COC response to Cd.
Concerning DEGs overexpressed in adult CCs, in the present study, significant CYP19A1 upregulation was found after Cd exposure, with an increase of about 250 times greater in Cd-exposed samples than in controls. This finding could mean that, in CCs, Cd may have acted as an endocrine disruptor, as the upregulation of CYP19A1, coding for aromatase, the most important enzyme involved in oestrogen biosynthesis, may have led to excess oestrogen levels, which is known to be related to deleterious effects on follicular development, oocyte and embryo quality [138]. This finding is in agreement with that of a previous study in mice, in which high levels of steroid hormones were found upon Cd exposure contributing to female sexual dysfunction with premature puberty onset [139]. However, CYP19A1 downregulation was found in response to Cd in other studies. In carp, it was showed that Cd affected steroidogenesis in a dose and time-dependent manner by reducing aromatase activity and expression [140]. Similarly, in another study, downregulation of aromatase expression in zebrafish female ovaries was found upon Cd exposure [46]. This apparent discrepancy on effects of Cd on aromatase expression could be linked to the different organism/cell regulation systems or Cd tested concentration.
According with deep sequencing results, Cd exposure increased the gene expression of IGFBP2 in CCs of adult sheep. Considering its function in ovarian follicles [141], the observed upregulation of IGFBP2 could lead to negative modulation of steroidogenesis in GCs and follicular atresia. To the best of our knowledge, this is the first study reporting the modulation of IGFBP2 expression following Cd exposure; therefore, our data can be discussed in the light of studies conducted in other cell systems or under different experimental conditions. IGFBP2 has been shown to block FSH-dependent E2 production in the ovarian follicle [64]. In our study, it cannot be excluded that the IGFBP2 overexpression could be a feedback mechanism induced by the excessive expression of CYP19A1 and the subsequent synthesis of oestrogens after Cd exposure.
The NOS2 gene, encoding for inducible nitric oxide synthase, was also upregulated by exposure to Cd in our study. It could be assumed that Cd damage on CCs is associated to inflammatory and oxidative state, as this gene is usually activated by inflammatory molecules or toxic elements. To the best of our knowledge, there are no studies that have analysed the effects of Cd on the expression of NOS2 in CCs; therefore, our data can be discussed in the light of studies conducted in other cell systems or under different experimental conditions. It has been demonstrated that Cd upregulated NOS2, resulting in increased nitric oxide (NO) production implicated in Cd-mediated cytotoxicity in mouse macrophages [142] and arteria [143], and in rat liver [144] and testis [145]. In human CCs, NOS mRNA expression has been negatively correlated with oocyte receptibility to fertilisation, as increased NOS2 expression was found in non-fertilized oocyte after ICSI compared to fertilized oocyte, so it could be considered a negative marker of oocyte competence [83]. On steroidogenesis, NO synthesis can directly inhibit aromatase both in human GCs and luteal cells [146] and this inhibitory effect on steroidogenesis can be considered a feedback mechanism induced by the excessive expression of CYP19A1 and the subsequent synthesis of oestrogen after Cd exposure.
The IHH gene, encoding for the Indian Hedgehog Signalling Molecule, was upregulated by exposure to Cd in our study. This signalling pathway has been shown to be involved in the regulation of several cellular activities, such as protein metabolic processes in bovine CCs [69], apoptosis in swine GCs [70] and GC differentiation in mice [73]. In our study, overexpression upon Cd exposure might indicate a response with a double meaning. On one side, IHH may have triggered an apoptotic process, since in a previous study, IHH was upregulated during the time of in vitro culture of the GCs, demonstrating an association with processes of aging and programmed GC death [70]. On the other side, Cd may have acted as an endocrine disruptor. In fact, in a previous study on women with polycystic ovary syndrome (PCOS), one of the most common endocrine disorders in women, IHH was abnormally highly expressed in the PCOS tissue [72]. To the best of our knowledge, the only study investigating IHH levels in relation to Cd effects was performed on zebrafish swim bladder [147]. Further studies are needed to clarify Cd effects on Hedgehog signalling pathway in CCs.
Among upregulated genes upon Cd exposure, there was MT1A involved in the regulation of metal ions transport. This finding is in agreement with previous studies in non-vertebrate organisms in which altered expression of the same gene or genes belonging to the same family were found, using transcriptomic approaches, in ovarian tissues [47]. Indeed, increased expression of MT1A was found which could be considered as a defensive cellular response activated by CCs after Cd exposure. In fact, due to the high affinity with zinc ions, MT1 plays an important role in the detoxification from heavy metals and free radicals produced during oxidative stress conditions [148]. In our experimental conditions, MT1A gene could have been overexpressed to provide increased amount of this protein to sequester excessive Cd for detoxification. Our results are in agreement with those obtained in previous studies in the mouse model in which increased expression of MT1 was observed in response to Cd absorption in liver [149], kidney, spleen, lung and heart, and this mechanism protected these organs from Cd toxicity [150]. In another study, uterine cells displayed increased expression of MT1A gene after Cd exposure, and this overexpression has been related to steroid estrogenic-dependent regulation [151]. Considering our findings, we could hypothesize that the promoted expression of MT1A in Cd-exposed CCs is a protective response of CCs against Cd-toxicity. Indeed, MT1A is a zinc-binding protein. When Cd is present, it binds to MT1A thus competing with zinc ions (ionic mimicry, namely the ability of a cationic form of a toxic metal to mimic an essential element or cationic species of an element at the site of a transporter of that element; [152] and removing Cd from the cell. Interestingly, this DEG was upregulated in adult and in prepubertal CCs. This finding led us to identify it as the most reliable biomarker of Cd-induced CC toxicity. Furthermore, MT2A, another member of metallothionein family, was found to be overexpressed. This gene also codes for a protein involved in homeostatic control and detoxification of heavy metals. Moreover, thanks to enrichment analysis, in adult and in prepubertal CCs, the involvement of many other members of MT family has been evidenced upon Cd exposure. Continuing among overexpressed genes in prepubertal CCs, we focused our attention on HSPA6, BMP15, DSG2 and SLC30A2. Although HSPA6 has been observed as deregulated in several studies upon Cd exposure [153][154][155][156], another study reported its overexpression in Cd-exposed human iPSC-derived renal cells, suggesting that other stressors (e.g., Cd) in addition to the heat stress, could regulate HSPA6 expression levels [157]. Similarly, HSPA6 expression was significantly upregulated in response to Cd exposure in human trophoblast cells [158], as a clear indication of activation of a defence system in response to a stressful factor.
Different studies have demonstrated that BMP15 is essential for normal follicular development [111][112][113] and for the acquisition of oocyte competence [159], through an efficient communication between the oocyte and its surrounding somatic cells [28,114]. Recently, BMP15 was also found to be overexpressed in the CCs of PCOS patients following drug treatments [160]. In our study, we identified a never shown before expression in CCs upon Cd exposure. The enrichment analysis identified co-expression with SLC30A3 and genetic interaction with SLC30A1; therefore, we could hypothesize that BMP15 overexpression in CCs was a stress-mediated response to the metal ion Cd, in addition to the other well-known detoxification systems mentioned already. To the best of our knowledge, to date, no studies have investigated the expression of this gene in relation to Cd effects, even in other cell systems, so further studies are needed to clarify these events.
DSG2 code for a desmoglein, a component of desmosomes, cell-cell junctions which have been reported as involved in apoptosis, upregulated genes in GCs from atretic follicles [116]. To the best of our knowledge, no previous studies have reported to date on effects of Cd on the expression of this gene, but its role in Cd-mediated CC toxicity can be hypothesized as it is well known that Cd affects intercellular cumulus communications and regular expansion.
According to deep sequencing results, Cd exposure increased the gene expression of SLC30A2 in CCs of prepubertal lambs. The Solute Carrier Family 30 Member 2 (SLC30A2) is a Zn transporter involved in Zn homeostasis within the cells [161]. Several studies have suggested that the toxic effects of Cd may be mediated by altered metabolism of essential elements, including Zn [162]. Zn homeostasis plays a pivotal role in the functioning of female germ cells. In murine COCs, there is a strict control of Zn levels within the oocyte supported by CCs [131], which actively regulate the uptake of Zn into the oocyte during maturation, which is necessary for the completion of meiosis and for MII arrest. In a previous study, Cd-exposure of rat mother caused depletion of Zn and upregulation of SLC30A2 gene expression in gastrointestinal tract and plasma of suckling puppets [162]. The validation procedure did not identify differences in the expression levels of SLC30A in CCs treated with Cd compared to controls. This result could be related to the use of independent CC samples for validation and transcriptomic analysis, which, although conducted under the same culture conditions, were derived from ovaries of different animals.
Interestingly, for all these genes, the enrichment analysis, in adult CCs, revealed interactions among all aforementioned genes involved in COC growth/differentiation, functional and detoxification pathways in response to Cd. Instead, in prepubertal CCs, independent clusters of genes, involved essentially in detoxification functions (MT family and SLC30) were identified. This difference could be explained considering that in prepubertal CCs, some functions could be not yet expressed or active (Table S4).
As a final consideration, downregulation was a lower intense response in this study, as it accounted for 23% of dysregulated genes. However, no less importance should be attributed to the findings of significantly downregulated genes upon Cd exposure among which, in particular, the presence of MYNC stands out as associated with oocyte maturation pathways [103], along with RGS4, which is downregulated in women with diminished ovarian reserves [105], and TSPAN18, which is involved in facilitating sperm inner acrosomal membrane apposition with the oolemma [107] (See Table S2 for further details).

Conclusions
In conclusion, our findings identified genes over-and under-expressed upon in vitro exposure to nanomolar Cd in CCs of adult and prepubertal oocytes in sheep. For genes previously associated with COC molecular developmental competence, the results of this study confirm and implement the knowledge on their involvement in these processes. Moreover, the novel identified genes can serve as a springboard to future studies. We plan to perform parallel studies aimed to identify, in humans and in the sheep model, shared biomarkers of Cd-induced infertility, to correlate them with clinical outcomes and develop appropriate detoxification strategies.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/biology12020249/s1, Table S1: Upregulated genes in CCs of adult sheep exposed in vitro to nanomolar cadmium, Table S2: Downregulated genes in CCs isolated from adult sheep and exposed in vitro to nanomolar cadmium, Table S3: Upregulated genes in CCs isolated from prepubertal lambs and exposed in vitro to nanomolar cadmium, Table S4: Detailed results of the gene network analysis.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.