This is the first study evaluating the transcriptomic changes in growing (G), dominant (D) and preovulatory (P) follicles in southern white rhinoceros (SWR). This information is needed to advance assisted reproductive technologies (ARTs) by determining appropriate methods for manipulating ovarian function prior to oocyte collection to obtain oocytes with the highest developmental potential. Through this analysis, we generated candidate gene markers for the physiological status of the follicle at the time of OPU.
Although this study was performed on SWR granulosa cells (GC), the SWR genome is incomplete and contains large gaps in annotation. Therefore, the northern white rhinoceros (NWR) genome and currently available annotations are the only resources for generating transcriptomic data on white rhinoceroses, but extensive manual curation is still necessary (15). In addition, the NWR and SWR are genetically similar at the chromosomal and mitochondrial genome levels, enabling the NWR genome to be used for this analysis (15). The original annotation for the NWR contained 14,274 different transcripts (12,252 genes) that were expressed in the GC from this dataset. After adding the transcripts from de novo sequencing, a total of 39,455 different transcripts (21,621 genes) were expressed. This was an overall improvement of 119% in transcripts and 61% in genes compared to the original annotation. In addition, it is important to note that through manual curation, multiple transcripts were reassigned to different genes, and transcripts with unknown associations were identified. This approach resulted in a substantial improvement to the available annotation, and these new annotations will be available publicly in the future. Due to their scarcity, reproductive cells are often underrepresented when generating annotations for genomes. Many new genes were annotated in the granulosa cells due to the massive transcriptional changes that occurred within the cells over time and their dynamic behavior. This additional information to the annotation will benefit other researchers in reproductive physiology in this species, and it emphasizes the importance of using reproductive cells, such as granulosa cells, as a resource for conservation research and genome improvement and accuracy. One of our objectives was to identify the genes most highly expressed across the follicle stages. Interestingly, all follicle stages had the same top three most highly expressed genes. These three genes were involved in macrophage immunometabolism regulator (MACIR), MSTRG15261.1, and collagen type I alpha 1 chain (COL1A1). The most highly expressed gene, MACIR, formally known as C5ORF30, has a well-documented role in modulating the immune response and regulating macrophage function (34, 35). In addition, blocking MACIR enhances the expression of genes involved in regulating cell migration, adhesion, and angiogenesis (36). In cattle GC, MACIR has been shown to be upregulated in follicles undergoing late-stage atresia (37). This is a potential indicator that the follicles in this study were under stress and already progressed through advanced, irreversible atresia. The second most highly expressed gene, MSTRG15261.1, was found to be unannotated in the SWR and in the domestic horse. The third most highly expressed gene was COL1A1, which is the most abundant collagen in the mammalian ovary and is regulated by transforming growth factor beta 1 (TGF-β1) (38). In addition, COL1A1 promotes the secretion of estrogen while also inhibiting the increase in progesterone, resulting in follicle growth (39) (40). Most importantly, COL1A1 plays an important role in oocyte maturation and embryo development through its association with the PI3K/Akt signaling pathway in mice and pigs (41–44). Interestingly, the most highly expressed gene was involved in follicle death, and the third highest was involved in follicle growth. This could indicate that prior to OPU, follicles received mixed signals during ovarian treatment.
Cluster analysis was performed to generate gene expression patterns across the follicle stages to identify genes that follow conventional literature for these stages and those that are not. Due to the limited opportunities to perform OPU in SWRs, the current preferred approach is to stimulate the growth of multiple follicles to increase the oocyte recovery rate. Therefore, unstimulated animals were not considered for this elective procedure in our study. Due to the lack of an unstimulated rhinoceros for comparison, the gene expression patterns were compared to published literature generated from unstimulated cows (28). The two main cluster types we focused on were patterns associated with follicle growth (clusters four and seven) and differentiation (clusters three, five, and six). For cluster three, the gene expression levels were consistent between G and D and increased from D to P; therefore, these genes are associated with follicle differentiation. Of the top 20 genes from this cluster, 13 have known expression patterns across follicle stages (28); interestingly, five of those genes should be downregulated in P, yet their expression is increased in our dataset. COL4A1 (collagen type IV alpha 1 chain) was found to be differentially expressed in this cluster compared to cattle data (28). In cattle, COL4A1 is a gene that stays in relatively low abundance throughout the follicle stages, with a decrease in P and an increase only observed in atretic follicles (28). COL4A1 is directly regulated by forkhead box L2 (FOXL2), a transcription factor specifically expressed in granulosa cells (45). High levels of FOXL2 decrease COL4A1 (45), which is required for the proliferation of GC, allowing follicle remodeling, and should remain high during folliculogenesis to allow for proliferation and differentiation (46). In our data, we observed a significant decrease in FOXL2 expression between D and P, which led to a significant increase in COL4A1 expression between D and P, a sign of late-stage atresia associated with follicular demise (31). In addition, high levels of collagen in the follicle basement membrane have been linked to an imbalance in gonadotropin stimulation (47). It is hypothesized that this delay in maturation due to increased collagen may lead to asynchrony during maturation and negative effects on oocyte developmental competence (47). Therefore, we can hypothesize that high levels of COL4A1 in our SWRs could be caused by incomplete or suboptimal stimulation prior to OPU.
Junction mediating and regulatory protein p53 cofactor (JMY) is a transcriptional cofactor and an actin nucleation factor (48) (49). Actin nucleation occurs during oocyte polarization and affects the microtubule and microfilament cytoskeleton by activating Arp 2/3 (for actin remodeling) and the assembly of filaments in the cytoskeletal spindle (48, 50, 51). If JMY is not present, the spindle fails to migrate to the cortex, and the oocyte arrests with a centrally located spindle, highlighting the crucial role of JMY during peripheral spindle migration during oocyte maturation (49). In both mouse and porcine oocytes, JMY decreases during oocyte maturation (49, 52). This same pattern of expression was observed in cluster five of our dataset, with JMY being highly significant within that cluster.
Cluster four included the genes whose expression increased from G to D to P. A gene of interest in this pattern was F-box and WD repeat domain containing 11 (FBXW11) due to its regulation by estradiol, its role in meiotic resumption, and its normal pattern of expression compared to our data. FBXW11 is regulated by many different factors, one of which is helicase for meiosis 1 (HFM1) (53). Low levels of HFM1 result in high levels of FBXW11, which corresponds to the pattern of expression we observed for both of these genes; however, according to the literature in mice, HFM1 expression should be high for oocytes to resume meiosis (53). In addition, FBXW11 is regulated by estradiol, as indicated by high levels of estradiol (as found in late stage developing follicles), which results in low levels of FBXW11 (28, 54). The primary function of FBXW11, formally known as β-TRCP2, is to degrade cell division cycle 25A (CDC25A), which is crucial for the resumption of meiosis (at the appropriate time) (55, 56). CDC25A must be downregulated for oocytes to be released from MII arrest (57), which we observed in the rhinoceros data. Therefore, from the literature, it appears that the timing of the increase in FBXW11 is important in relation to the resumption of oocyte meiosis.
NRG1 (neuregulin 1) is a member of the EGF-like factor family that is induced by gonadotropins (58, 59). There is a dramatic increase in NRG1 expression due to LH induction, which regulates proper oocyte maturation and developmental competence (60). Recent literature has shown that NRG1 supplementation during IVM can promote oocyte nuclear maturation and blastocyst development and hatching in cows (60). During development, amphiregulin (AREG) and NRG1 work synergistically (AREG is an accelerator, and NRG1 is a brake) to ensure that the required temporal changes within the oocyte occur for proper meiotic resumption and competence (59). It is important that these genes are balanced to ensure that oocytes have high developmental competence. According to our data, NRG1 is part of cluster five and remains consistent from G to D but decreases from D to P, which is the opposite of what the literature states regarding high levels of NRG1 being needed for oocyte competence acquisition (28, 60). NRG1 is a novel factor that may impact oocyte maturation; therefore, it should be considered a biomarker in GC for oocyte maturation competence and should be added to in vitro maturation medium.
Cluster six was classified as a cluster related to differentiation because there was a decrease in gene expression from G to D and a large increase from D to P. A gene that followed this cluster pattern but, based upon the literature should not be, was thymopoeitin (TMPO), also known as lamina-associated polypeptide 2 (LAP2) (61). According to both bovine and porcine data, TMPO decreases as follicle progression progresses (28, 62), which we previously reported from G to D but subsequently deviates with an increase from D to P. A specific region in TMPO [regions 32–36 known as thymopentin (TP5)] has been studied extensively due to its role in the activation of the pluripotency marker LIN28a in oocytes (63, 64); therefore, high levels of TMPO result in high levels of LIN28a. LIN28a expression decreases until the 8-cell stage of embryo development, when a dramatic increase is then noted (63). In contrast, LIN28a levels should be low until after resumption of meiosis, and high LIN28a levels in oocytes are associated with MII-arrested oocytes in mice (63, 65). High levels of LAP2 have also been associated with DNA damage in oocytes and programmed cell death in granulosa cells in mice (61). Overall, these data indicate that the P follicles in this dataset may have DNA damage, potentially leading to oocyte arrest at MII.
Based upon all the cluster analyses, the largest digression from the literature appears to be in the D to P transition, corresponding to differentiation. The five main drivers of differentiation in cattle are TGF-β1, tumor protein 53 (TP53), estradiol, hepatocyte nuclear factor 4 alpha (HNF4A) and retinoic acid (RA) (28). Our dataset included all but RA, and the expression of each gene significantly increased from D to P. RA has a positive regulatory effect on granulosa cell proliferation and oocyte maturation (66). RA levels in growing follicles are required for the follicle to reach appropriate final maturation (differentiation) in response to gonadotropins and affect gonadotropin receptor expression in granulosa cells in cows (67). It can therefore be hypothesized that with this ovarian treatment protocol, follicles are not receiving the proper stimuli to fully differentiate and produce oocytes with the highest developmental competence potential.
Based upon these results, it appears that the ovarian treatment protocol used on these animals did not result in GC or oocytes with the highest developmental competence potential. It is important to note that these comparisons are based primarily on data from unstimulated cows. We hypothesize that the gene expression changes we observed could be an outcome of the simulation protocol, but it is unknown whether some of these deviations are due to the unique expression pattern of the rhinoceros from other species. Because we performed oocyte group culture in vitro, we cannot directly correlate the in vivo GC analyzed with the corresponding oocyte maturation outcome; hence, in the future, we will match GC to the in vitro maturation outcome of oocytes by performing individual culture. Overall, for the animals treated with this stimulation protocol, we achieved an in vitro maturation rate of 33%, but none of those oocytes fertilized or developed into embryos. Although a percentage of retrieved oocytes were capable of extruding the first polar body and reaching metaphase II (MII) after in vitro maturation, we hypothesize that unpaired cytoplasmic maturation, erroneous chromosomal segregation, or morphological characteristics of atresia (68) left the oocyte incapable of fertilization and becoming a zygote. Incomplete cytoplasmic maturation can lead to failed sperm decondensation, erroneous male pronuclear formation, and ultimately unsuccessful fertilization (69). However, additional studies are needed to characterize the effects of ovarian stimulation on oocyte developmental competence acquisition/failure and cytoplasmic maturation in this species. A balance between the quantity and quality of oocytes obtained during OPU is suggested, and low, frequent doses of GnRH might be a valuable treatment.
With this dataset, we generated a list of potential biomarkers of the GC for the physiologic status of the follicle to determine which oocytes had the highest potential for developmental competence. These markers include COL1A1 (follicle growth) (41–44), JMY and FBXW11 (oocyte preparedness to resume meiosis) (49, 53, 57), NRG1 (oocyte maturation potential) (60), TMPO (DNA damage) (61), and MACIR and COL4A1 (markers for cell death and atresia) (31, 46). However, further testing will be needed to correlate these markers with the status of the follicles we collected during OPU and ultimately with oocyte maturation and cleavage.
It is important to note that one limitation in this study is grouping the follicle stages were grouped based upon follicle size at OPU. Nevertheless, for this study, we were not able to obtain follicles at all stages from the same animal due to the unknown status of the estrous cycle prior to ovarian stimulation and OPU. We attempted to overcome this limitation by obtaining technical replicates from different granulosa cell aliquots of the preovulatory follicle. In addition, the animals in this study were not ultrasounded regularly prior to ovarian treatment and OPU; therefore, the trajectory of the follicle prior to OPU was unknown (i.e., was it increasing in size or was it already regressing and undergoing atresia). We attempted to overcome this by evaluating multiple follicles per stage, and in the future, we will attempt to monitor these animals prior to OPU.
In conclusion, this was the first study to evaluate the transcriptomic differences in GC between the different stages of follicle development in the SWR, an important chapter in the reproductive knowledge of this species. In this threatened species, much remains unknown regarding ovarian physiology and follicle activity; nevertheless, it is important to continue addressing these reproductive questions to support assisted reproductive technologies for rhinoceros species. Finally, these data provide multiple biomarkers to potentially improve OPU results. This body of work contains novel information on the SWR and enriches the reproductive knowledge of the species, which will support conservation efforts for this species.