Corrigendum: Construction of a circRNA-lincRNA-lncRNA-miRNA-mRNA ceRNA regulatory network identifies genes and pathways linked to goat fertility

Farzad Ghafouri ^1* Mostafa Sadeghi ^1* Abolfazl Bahrami ^1* Masoumeh Naserkheil ^1* Vahid Dehghanian Reyhan ^1* Arash Javanmard ^2* Seyed Reza Miraei Ashtiani ¹ Soheila Ghahremani ^3* Herman Wildrik Barkema ⁴ Rostam Abdollahi-Arpanahi ^1* John P Kastelic ⁴

• ¹ University of Tehran, Tehran, Tehran, Iran
• ² University of Tabriz, Tabriz, East Azerbaijan, Iran
• ³ Tarbiat Modares University, Tehran, Tehran, Iran
• ⁴ University of Calgary, Calgary, Alberta, Canada

1 INTRODUCTION Improvement of reproductive performance is a priority in the goat industry, as it is one of the most important determinants of productivity, sustainability, and profitability. The reproductive cycle involves dynamic and complex ovarian functions, characterized by progressive emergence and development of ovarian follicles (Fatet et al., 2011) under endocrine control. Both granulosa and theca cells are involved in steroidogenesis (Qiu et al., 2013). Granulosa cells (GCs) are one of the most important cell types in ovarian follicles, with crucial roles in follicular development and atresia (Matsuda et al., 2012), especially in the late stages of oocyte development and ovulation. In addition, these cells also control cytoplasmic maturation and have a key role in nuclear maturation by responding to gonadotropins (Mori et al., 2000). Thus, reproductive function in the female is inherently complex, involving various anatomical and physiological processes (Li et al., 2021a). Furthermore, litter size, an important attribute directly related to reproductive efficiency, is controlled by multiple genes and factors (Lai et al., 2016). Hence, knowledge of the genetic basis of reproductive efficiency will provide insights into components controlling ovarian folliculogenesis and fertility in goats (de Lima et al., 2020). The candidate gene approach for fertility has been extensively studied in various livestock species (Miao et al., 2016a,b; Bahrami et al., 2017; Quan et al., 2019). Furthermore, it requires well-developed tools to detect and characterize multiple genes, pathways and networks (Ahlawat et al., 2016; Zhang et al., 2017; Ghafouri et al., 2021, Naserkheil et al., 2022). RNA sequencing (RNA-Seq) has been used to characterize transcripts and differences in gene expression, identify functional genes, and analyze regulatory networks in numerous species (Chen et al., 2015; Li et al., 2021b,c). In addition, single-cell RNA sequencing (scRNA-seq) is being used for mapping and quantifying transcriptional activity at single-cell resolution for all genes in the genome (Islam et al., 2014). It is useful for analysis of cellular heterogeneity, as it can concurrently sequence thousands of cells, as well as discover novel cell types in animals (Choi and Kim, 2019). Conversely, an integrated approach is needed to manage large-scale data generated with high-throughput technologies alongside literature mining. Integrated analyses can combine multilevel views of physiology data into a total interpretation of nonlinear regulatory molecular procedures (La et al., 2019; Sadeghi et al., 2022; Reyhan et al., 2022). Currently, various bioinformatics tools, computational approaches, and algorithms are available to identify interactions and protein functions in regulatory modules in various complex biological networks (Bugrim et al., 2004). In this regard, multi-partite networks such as CircRNA–lincRNA–lncRNA–miRNA–mRNA ceRNA regulatory networks consider various RNAs and have highlighted a new regulatory mechanism of interaction between RNAs. Circular RNAs (CircRNAs) are single-stranded, covalently closed RNA molecules without free 5' and 3' ends that act as transcriptional regulators, microRNAs, and protein templates (Zhou et al., 2020). Long intergenic non-coding RNAs (lincRNAs), RNA transcripts with >200 nucleotides, have major roles in biological processes such as gene expression control, epigenetic control, and scaffold formation (Deniz and Erman, 2017). Long non-coding RNAs (lncRNAs) are non-coding RNA transcripts involved in various biological procedures, such as cell proliferation and transcriptional regulation (Wei et al., 2016). Also, microRNAs (miRNAs), regulatory molecules with 19–25 nucleotides, have vital regulatory roles in multiple biological procedures (cell differentiation and migration, oncogenesis, and apoptosis) by suppressing mRNA (Wang et al., 2015). Therefore, this approach seems well suited to understand molecular regulatory mechanisms in polygenic traits (Hallock and Thomas, 2012). Several studies have identified important candidate genes associated with hormonal regulation of the reproductive cycle and fertility traits in goats (Li et al., 2010; Lai et al., 2016; Ahlawat et al., 2016; Zhang et al., 2018a; Li et al., 2021b). Furthermore, gene ontology and systems biology enable identification of hub genes and co-expression genes with critical roles in fertility (Ahlawat et al., 2016; Zhang et al., 2018b). For instance, in a study comparing high- versus low-fertility goats, many candidate genes were identified in each group. In an analysis of the entire genome of Chinese Laoshan dairy goats, several candidate genes (CCNB2, AR, SMAD2, AMHR2, KDM6A, SOX5, and SYCP2) were associated with both high and low fertility (Lai et al., 2016). Therefore, the purpose of the present study was to examine the ovarian granulosa cell signature genomic regions of Jining Grey goats with high and low fertility, using a public repositories ScRNA-Seq dataset. In addition, we critically reviewed the literature, searching for relevant publications using keywords related to granulosa cells and fertility in goats. Genes with significant differences related to our bioinformatics analyses and candidate gene list extracted from literature mining were compiled, and the 2 gene sets were merged and used to identify protein–protein interaction networks (PPI), and gene regulatory networks (GRNs). Overall, the ceRNA regulatory network, consisting of circRNA–lincRNA–lncRNA–miRNA–mRNA, was constructed based on various interactions to explore effects of functional modules and hub differentially expressed mRNAs, miRNAs, lncRNAs, lincRNAs, and circRNAs on fertility. 2 MATERIALS AND METHODS As an overview, the general workflow for analyzing data collection and methods of identifying key genes, metabolic and signaling pathways, as well as construction of the lncRNA–miRNA–mRNA ceRNA regulatory network and Modules affecting HF and LF in domestic goats (Capra hircus) are presented in Figure 1. 2.1 Data collection Single-cell RNA-Sequencing (ScRNA-Seq) data from ovarian granulosa cells of high- and low-fertility Ji’Ning grey goats (HF and LF, respectively) were retrieved from the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) public database under accession number GSE135897 (www.ncbi.nlm.nih.gov/geo). This dataset was produced using the GPL15473 Illumina HiSeq 2000 platform (Li et al., 2021b). A case-control study approach was designed to identify differentially expressed genes (DEGs) between HF and LF goats. A total of 30 samples from Ji'ning gray goats (first group: 15 with high litter size (HF; ≥3 offspring) and second group: 15 with low litter size (LF; ≤2 offspring)) were analyzed. The Ji’Ning grey goat is a local breed and the oldest domesticated goat species with high fertility in China that has year-round estrus and mean litter size of 2.94 (Huang et al., 2012; Miao et al., 2016a,b). All ovarian follicles from each group of these Ji’Ning grey goats (from small (<3 mm) to large (˃7 mm)) were collected. Then, oocytes and GCs were mechanically separated by repeated pipetting. Oocytes from high- and low-fertility groups were labeled. This study was performed in accordance with ARRIVE guidelines. Library preparation and sequencing were performed separately for oocytes and GCs. Details regarding animal ethics approval, GCs, total RNA extraction, and sequencing and validation of RNA-Seq data have been reported (Zhao et al., 2020; Li et al., 2021b). 2.2 Quality control and detection of differentially expressed genes First, quality of the raw RNA sequences was assessed using FastQC software v0.11.5 (Andrew, 2010) and FastQ Groomer software v1.1.5 (Blankenberg et al., 2010); thereafter, these sequences were pre-processed with Trimmomatic software v0.38.0 (Bolger et al., 2014) to remove adapters, low-quality reads, and PCR primers. Alignment sequences, mapping, and identification of known and novel RNAs of reads were related to the reference genome of Capra hircus (https://ftp.ensembl.org/pub/release-108/fasta/capra_hircus/dna/) using HISAT2 software v2.2.1 with default parameters to determine number of aligned and unaligned reads (Kim et al., 2015). Regarding transcript quantification, total raw counts of mapped reads were calculated using featureCounts software (v2.0.1) (Liao et al., 2014). Subsequently, to examine whether accumulation or degradation of transcripts was related to fertility, transcripts and their expression levels were compared between ovarian granulosa cell samples of high- and low-fertility Ji’Ning grey goats. Differences in gene expression were detected from reading counts using DESeq2 software (v2.11.40.7) (Love et al., 2014). The threshold for statistical significance of the differential expression of each gene was obtained with the criteria of log2FC (fold change <5 and >5) and FDR <0.05 (false discovery rate). 2.3 Literature mining to identify candidate genes for fertility in goat We reviewed the literature up to March 12, 2023, searching for relevant publications using keywords related to granulosa cells and fertility in goats. Regulatory RNAs (i.e., circRNAs, lincRNAs, lncRNAs, miRNAs, and mRNAs) with significant differences related to candidate RNAs list extracted from literature mining were compiled as RNA sets 2-6, respectively. Finally, RNA set 1 (from our bioinformatics analyses) and RNA sets 2-6 (from literature mining) were merged and used as input files for target prediction tools, STRING website, and Cytoscape software to identify protein-protein interactions gene regulatory networks and reconstruct the circRNA–lincRNA–lncRNA–miRNA–mRNA ceRNA regulatory network, as well as modules. 2.4 Gene ontology and functional enrichment analysis Gene set annotation and enrichment analysis used DAVID (the Database for Annotation, Visualization, and Integrated Discovery; https://david.ncifcrf.gov/) v6.8 (Sherman and Lempicki, 2009), g:Profiler (https://biit.cs.ut.ee/gprofiler/gost) (Reimand et al., 2016), GeneCards (https://www.genecards.org/), and STRING database v11.0 (https://string-db.org) (Szklarczyk et al., 2019) to determine potential functions and metabolic pathways. Genes were assigned to functional categories using the Gene Ontology (GO) database under biological process (BP), molecular function (MF), and cellular component (CC). 2.5 Target prediction of differentially expressed mRNAs and other types of regulatory RNAs Functional annotation of types of regulatory RNAs i.e., circRNAs, lincRNAs, lncRNAs, and miRNAs consisted of functional annotation of their potential target mRNA genes. Predicted targeted genes and types of regulatory RNAs were predicted using: miRBase (Kozomara and Birgaoanu, 2019) (https://www.mirbase.org/); Targetscan (Grimson et al., 2007); miRanda (http://www.microrna.org/); miRWalk 3.0 (a comprehensive atlas of microRNA–target interaction tools that integrates 12 miRNA–target prediction tools; http://mirwalk.umm.uni-heidelberg.de/); NONCODE database (Volders et al., 2019) (http://www.noncode.org/); LNCipedia database (Bader et al., 2006) (https://lncipedia.org); and CircInteractome web tool (Dudekula et al., 2016) (a computational tool that enables prediction and mapping of binding sites for RBPs and miRNAs on reported circRNAs; https://circinteractome.nia.nih.gov/). Identified target genes were selected and submitted to DAVID, KEGG, Reactome pathways, and the PANTHER database for enrichment and validation of target genes for each type of RNA. 2.6 CircRNA–lincRNA–lncRNA–miRNA–mRNA ceRNA regulatory network construction Based on the ceRNA theory, global functions for all non-coding RNAs can serve as an endogenous “sponge” to regulate upregulated or downregulated circRNAs, lncRNAs, lincRNAs, miRNAs or mRNAs that are inverse relationships together in the mRNA-mRNA, mRNA-miRNA, mRNA-lncRNA, mRNA-lincRNA, mRNA-circRNA, miRNA-lncRNA, and miRNA-circRNA interaction pairs chosen to construct the circRNA–lincRNA–lncRNA–miRNA–mRNA ceRNA regulatory network (Xu et al., 2019). In this regard, protein-protein interaction (PPI) network analysis was performed using the STRING database v11.0 (Szklarczyk et al., 2019) (Search Tool for the Retrieval of Interacting Genes or Proteins; https://string-db.org), BIND (Biomolecular Interaction Network Database) (Bader et al., 2003), MIPS (Mammalian Protein-Protein Interactions Database) (Pagel et al., 2005), and BioGRID (Biological General Repository for Interaction Datasets) (Chatr-Aryamontri et al., 2012) to explore interactions between genes in Capra hircus. After identifying interactions between types of regulatory RNAs and gene expression data (Co-Expression), the circRNA–lincRNA–lncRNA–miRNA–mRNA ceRNA regulatory network was reconstructed and plotted using Cytoscape software v3.9.1. (Shannon et al., 2003; National Institute of General Medical Sciences, Bethesda Softworks, Rockville, MD, USA). Also, statistical and topological significance of the network was assessed with the Network Analyzer plugin in Cytoscape software. 2.7 Clustering of circRNA–lincRNA–lncRNA–miRNA–mRNA ceRNA regulatory network and identification of main hub regulatory RNAs Modules or subnets may have a significant role in the biologically rebuilt main ceRNA regulatory network, as they represent a set of nodes with similar functions that pursue specific biological purposes as functional modules. To evaluate topological properties and cluster nodes of the circRNA–lincRNA–lncRNA–miRNA–mRNA ceRNA regulatory network, the Cytoscape plugins ClusterONE (Shannon et al., 2003) and MCODE (Lotia et al., 2013), clustering algorithms to draw directional and directionless graphs, was used. ClusterONE is a plugin to discover densely connected sub-graphs of a network by minimizing edges between clusters and maximizing edges within a cluster. Also, the MCODE plugin can be used for directed or undirected graphs. Various parameters and the statistical and topological significance of the ceRNA regulatory modules were calculated using the Cytoscape plugin Network Analyzer v4.4.8 with the default for the directed network (Assenov et al., 2008). 3 RESULTS AND DISCUSSION 3.1 ScRNA-Seq analysis to identify differentially expressed genes (DEGs) In this study, we investigated the pattern of transcriptome profiles of ovarian granulosa cells in goats. To provide insights into the genetic basis of fertility in Ji'ning gray goats, gene expression analyses using the ScRNA-Seq dataset with access number GSE135897 (obtained from the GEO database) were used. The selected gene expression profile had 30 samples, including 15 samples of granulosa cells from high-fertility goats and 15 samples of granulosa cells from low-fertility goats. To perform this analysis, these 30 samples were divided into 2 groups (high- and low-fertility) to compare expressions of gene profiles and identify significant DEGs. A total of 3,245 significant genes were identified by processing the expression profile of granulosa cells of high- versus low-fertility goats. Finally, by considering the expression change threshold (Fold change <5 and >5, FDR <0.05), 150 genes were significantly differentially expressed in granulosa cells from goats with high versus low fertility. Of these genes, 80 genes were up-regulated and 70 were down-regulated (Supplementary Table S1). Recently, many studies in molecular genetics, bioinformatics, and biological systems have been conducted to discover candidate genes and identify molecular pathways involved in fertility (Ahlawat et al., 2020; An et al., 2021; Li et al., 2021b). Most identified genes were related to reproductive functions and strongly regulated ovarian follicular growth and secretion of hormones involved in fertility; consequently, increasing or decreasing their expression at various times can lead to complex biological processes during pregnancy (Wang et al., 2019). In this regard, in a study using RNA-Seq analysis for ovarian tissue of pregnant and non-pregnant goats, 4 genes (PGR, PRLR, STAR, and CYP19A1) were identified as having important roles in goat reproduction (Quan et al., 2019). Also, another study concluded that some lncRNAs in goats have key roles in regulating follicle development and cell growth during ovarian development (Li et al., 2021c). 3.2 Literature mining-based evidence for identified DEGs and types of regulatory RNAs Literature mining provided evidence for 81 well-known DEGs that were not included in our investigated outputs; adding this list of DEGs to our discovered DEGs created a good platform for further downstream analyses, including Gene Ontology (GO) and functional enrichment analysis. Identified DEGs based on literature mining and their role in fertility in goats are collected as RNA set 2 (Supplementary Table S2). Also, literature mining provided evidence for 58, 8, 19, and 55 well-known types of regulatory RNAs, i.e., circRNAs, lincRNAs, lncRNAs, and miRNAs with significant differences, compiled as RNA sets 3-6, respectively (Supplementary Table S3-6). The gene list in Supplementary Table S2 was merged with the gene list in Supplementary Table S1 and used to identify intergenic interactions and reconstruct the protein-protein interaction network. Then, after identifying interactions between types of regulatory RNAs in RNA sets 3-6 together and with gene expression data, the circRNA–lincRNA–lncRNA–miRNA–mRNA ceRNA regulatory network was reconstructed. 3.3 Gene ontology and pathway enrichment analysis The functional annotation of GO terms was detected based on the biological process (BP), molecular function (MF), and cellular component (CC) to identify functions and metabolic pathways as well as systematic features of the merged list of 231 DE genes (RNA sets 1 and 2), using the STRING, DAVID, PANTHER, and g:Profiler databases. Results of the GO classification of the DEGs for high- and low-fertility goats are presented in Figure 2. Identified DEGs were significantly involved in the following functions: >10 genes from DEGs had roles in cellular process, regulation of biological process, response to stimulus, cellular response to stimulus, cell communication, regulation of cellular metabolic process, cellular component organization or biogenesis, reproductive process, intracellular signal transduction, transforming growth factor beta receptor signaling pathway, BMP signaling pathway, response to steroid hormone, and SMAD protein complex assembly for biological processes (BP) (Figure 2A). The gene list, including genes SMAD2, ESR1, SOX5, BMP4, BMP15, CTNNB1, ERBB2, FGFR1, CDH26, GH, AR, and FSHB were involved in most of the molecular function terms that can be considered significant genes. Moreover, 16 molecular functions (MF) terms were identified, such as ion binding, transmembrane receptor protein kinase activity, signaling receptor binding, molecular function regulator, receptor-ligand activity, growth factor binding, and transforming growth factor beta-activated receptor activity, which were the most significant functions associated with fertility (Figure 2B). Regarding cellular components, 5 GO terms, including extracellular space, plasma membrane region, receptor complex, membrane raft, and caveola were identified (Figure 2C). In addition, KEGG-based path analysis for identified DEGs was performed using 3 online databases, DAVID, STRING, and g: Profiler. The identified DEGs involved in fertility were enriched in the signaling pathways regulating pluripotency of stem cells, cytokine-cytokine receptor interaction, ovarian steroidogenesis, oocyte meiosis, progesterone-mediated oocyte maturation, parathyroid hormone synthesis, secretion and action, growth hormone synthesis, secretion and action, cortisol synthesis and secretion, prolactin, TGF-beta, Hippo, MAPK, PI3K-Akt, and FoxO signaling pathways (Figure 3). Concerning female reproductive function, CTNNB1, BMP4, FSHR, TGFB1, BMPR1B, and ESR1 genes are jointly encoded in functions of ovulatory cycle process, ovarian follicular development, developmental process involved in reproduction, and cell differentiation. Hence, according to genes involved in the identified pathways, ESR1 and BMPR1B were considered candidate genes for reproductive functions, growth, and cell differentiation in goats. However, many studies have implicated the BMPR1B gene as one of the main genes controlling reproductive function and fertility (e.g., ovulation rate) in small ruminants, especially goats (Pramod et al., 2013; Ahlawat et al., 2014). In addition, the FSHR (Follicle-Stimulating Hormone Receptor) gene has roles in growth, differentiation, and maturation of follicles and enhances reproductive function in goats and sheep (Chen et al., 2017). Moreover, another study demonstrated that FSHR is involved in differential expression of ovarian mRNA hormone receptor genes in goat fertility (Saraiva et al., 2011). The ESR1 gene encodes estrogen receptor and ligand-activated transcription factor, and regulates the main genes involved in growth, metabolism, and pregnancy. It has also been reported that expression of this gene was highest in kidney, ovary, uterus, and testes, but lowest in brain and heart tissue (Mohammadabadi, 2020). Given the role of the BMP4 gene in reproduction, especially in growth and differentiation of ovarian follicles and ovulation, it can be considered a main candidate gene for reproductive function and fertility (Sharma et al., 2013). The CTNNB1 gene encodes a complex of proteins that constitutes adherens junctions (AJs); they are necessary for formation and maintenance of epithelial cell layers by regulating cell growth and intercellular connectivity. Also, they have an essential role in reproduction and multiplication (Zhang et al., 2018b). Therefore, our findings regarding genes FSHR, CTNNB1, BMPR1B, and ESR1 were consistent with other studies. 3.4 Construction and clustering of circRNA–lincRNA–lncRNA–miRNA–mRNA ceRNA regulatory network A decade ago, Salmena et al. presented the competitive endogenous RNA hypothesis (Salmena et al., 2011). The ceRNA regulatory networks have provided a new mechanism of interaction between RNAs and have crucial roles in multiple biological processes (Gao et al., 2021; Kfir et al., 2018; Yang et al., 2019). In this regard, many studies have been dedicated to elucidating the ceRNA roles of non-coding RNAs in some economically important traits by constructing competitive endogenous RNA networks (Salmena et al., 2011; Han et al., 2020). To detect the mechanism of how non-coding RNAs (ncRNAs) regulate mRNA through sponging miRNA, a ceRNA regulatory network was constructed with a merge of predicted mRNA-miRNA, mRNA-lncRNA, mRNA-lincRNA, mRNA-circRNA, miRNA-lncRNA, and miRNA-circRNA interaction pairs. The reconstructed ceRNA regulatory network for up- and down-regulated mRNAs/genes and types of regulatory RNAs, indicating physical connections between 2 or more protein molecules related to biochemical functions, is presented in Figure 4. Based on knowledge of interactions, this ceRNA regulatory network consisted of 310 nodes and 758 edges and the associated files with the networks are given in Supplementary Table S7 (stored in “.cys format” for further analyses). In detail, 57 circRNAs, 8 lincRNAs, 19 lncRNAs, 51 miRNAs, and 175 mRNAs were included in the network (Figure 4). As mentioned, molecular species (circRNAs, lincRNAs, lncRNAs, mRNAs, and miRNAs) in constructed networks are indicated as nodes and interactions between them as edges. Moreover, constructed networks were combined in a simple interaction format (SIF) using to Cytoscape (v3.9.1.) for topological analysis. Topological parameters of the ceRNA regulatory network and modules such as the number of nodes, number of edges, clustering coefficient, characteristic path length, and network density were evaluated to examine the state of communication and information transfer of a node with other nodes of interactive networks, as presented in Table 1. In this ceRNA regulatory network, 18 genes (SMAD1, SMAD2, SMAD3, SMAD4, TIMP1, ERBB2, BMP15, TGFB1, MAPK3, CTNNB1, BMPR2, AMHR2, TGFBR2, BMP4, ESR1, BMPR1B, AR, and TGFB2) had the most interactions with other genes in the network. Among these 18 hub genes, 4 genes (BMP4, BMPR1B, CTNNB1, and ESR1) were involved in molecular function and biological processes, in agreement with other studies (Sharma et al., 2013; Pramod et al., 2013; Zhang et al., 2018a; Mohammadabadi, 2020). Among the miRNAs, chi-miR-423-5p, chi-miR-122, chi-miR-187, and chi-miR-133b suppressed most of the identified hub genes as targets of the selected miRNAs. Also, four of the lincRNAs i.e., ENSCHIG00000000609, ENSCHIG00000000641, ENSCHIG00000000886, and ENSCHIG00000002761 interacted with TGBFR2, CTNNB1, TGFB2, and SMAD2 hub genes, respectively. Conversely, the miRNAs such as chi-miR-92a-5p, chi-miR-21-3p, chi-miR-202-3p, and chi-miR-223-3p interacted with most of the lncRNAs involved in ceRNA regulatory network and were defined as hub miRNAs that interacted with lncRNAs. In this study, clustering of the circRNA–lincRNA–lncRNA–miRNA–mRNA ceRNA regulatory network was performed using ClusterONE (Shannon et al., 2003) and MCODE plugins (Lotia et al., 2013), according to clustering algorithms utilized to determine significant sub-networks or modules by an integrated approach. There were 4 candidate modules involved in goat fertility; the node interactions of each of is described in Supplementary Table S7. Interestingly, module 1 consisted of 70 nodes and 282 edges and in detail comprised 2 lincRNAs, 12 miRNAs, and 56 mRNAs. In this module, SMAD1, SMAD2, SMAD3, SMAD4, TGFBR2, CTNNB1, and ESR1 were hub genes. Also, chi-miR-128, chi-miR-122, chi-miR-187, chi-miR-200a, chi-miR-206, and chi-miR-133b as hub miRNAs suppressed most of the involved genes such as DNMT3B, SMAD2, SMAD4, CHRD, RBCK1, BMPR1B, KDM6A, KRR1, GNA13, ERBB3, TGFB2, ERBB2, SPSB1, TGFBR2, BMPR2, CCL21, THBS1, PRLR, ESR1, and RASGRP1. Moreover, in this module, ENSCHIG00000000609 and ENSCHIG00000000886 lincRNAs interacted with TGFBR2, and TGFB2, respectively. All hub-hub genes in this module were involved in the metabolic-signaling pathways that were analyzed (Figure 5). Additionally, module 2 consisted of 42 nodes and 85 edges and in detail comprised 2 lincRNAs, 8 miRNAs, and 32 mRNAs. In this module, SMAD2, SMAD4, CTNNB1, and TIMP1 were hub genes. Also, chi-miR-182, chi-miR-200a, chi-miR-187, and chi-miR-122 as hub miRNAs suppressed most of the genes such as SMAD2, TRIP10, COPS2, PAK2, KMT2A, ERBB2, APC, GNA13, USB1, and FZD6. Among these, the SMAD2 gene was suppressed more than other genes involved in the network. Moreover, in this module, ENSCHIG00000000641 and ENSCHIG00000002761 lincRNAs interacted with CTNNB1, and SMAD2, respectively (Figure 6). Furthermore, module 3 consisted of 77 nodes and 112 edges and in detail comprised 1 circRNAs, 7 lncRNAs, 16 miRNAs, and 53 mRNAs. In this module, TUBB4B, CETN2, XAB2, and 3BHSD had the most interaction with other module genes as hub genes. Also, chi-miR-200a, chi-miR-494, and chi-miR-128 as hub miRNAs suppressed genes ANAPC7, KMT2A, PAK2, MARVELD2, SLC27A1, TADA1, FZD3, KDM6A, FSHB, and PKNOX1. Among these, chi-miR-128 miRNA suppressed most of the genes involved in the network. Moreover, in this module, XR_001297559.1 and XR_001297560.1 circRNAs interacted with most miRNAs such as chi-miR-494, chi-miR-3959-5p, chi-miR-494, chi-miR-1, chi-miR-202-5p, chi-miR-34c-5p, chi-miR-320-3p, and chi-miR-136-3p (Figure 7). Finally, module 4 consisted of 44 nodes and 59 edges and in detail comprised 5 lncRNAs, 2 miRNAs, and 37 mRNAs. In this module, MAPK3, JAK1, and ERBB3 had the most interaction with other module genes as hub genes. Also, almost all genes involved in this module were suppressed by chi-miR-423-5p and chi-miR-187 miRNAs. Moreover, chi-miR-423-5p miRNA interacted with all of the 5 lncRNAs i.e., XR_001297560.1, XR_001297559.1, LNC_000292, LNC_000417, and LNC_000492 involved in the module (Figure 8). Based on the literature, genes related to reproductive function in dairy goats, including CCNB2, AR, ADCY1, DNMT3B, SMAD2, AMHR2, ERBB2, and FGFR1 genes, were specifically selected in goats with high fertility (Lai et al., 2016; Zonaed Siddiki et al., 2020). Therefore, our results provided further evidence of the association between 8 hub genes (AR, SMAD1, SMAD2, SMAD3, SMAD4, AMHR2, TGFBR2, CTNNB1 and ERBB2) and female reproductive functions in goats. It was reported that the ERBB2 (Erb-b2 receptor tyrosine kinase 2) gene is a steroid hormone receptor involved in calcium signaling (Zwick et al., 1999). Similarly, the AR (androgen receptor) gene is a protein that involves a DNA binding transcription factor and chromatin binding; furthermore, it has a key role in reproduction by transmitting androgen signals (Zonaed Siddiki et al., 2020). TGFBR2 (Transforming Growth Factor Beta Receptor 2) is a protein-coding gene that has a major role in TGF-beta receptor signaling, activating SMADs, transferase activity, transferring phosphorus-containing groups, and protein tyrosine kinase activity (Yao et al., 2022). CTNNB1 (catenin beta 1) and TGF-β2 genes were associated with transforming growth factor-beta (TGF-β), signaling pathways regulating pluripotency of stem cells, HTLV-I infection, neuroactive ligand-receptor interaction, Wnt, and Hippo signaling pathways (Li et al., 2021b). The SMAD2, SMAD3, and SMAD4 genes have critical roles in growth and differentiation of ovarian cells, consistent with some aspects of ovulation (Li et al., 2008; Fortin et al., 2014). Notably, the AMHR2 (Anti-Mullerian Hormone Receptor Type 2) gene is associated with transferase activity, transport of phosphorus-containing groups, and protein tyrosine kinase activity. This gene is also involved in growth and development of ovarian follicles in cattle and goats (Monniaux et al., 2011). Furthermore, genes BMP15 and GDF9 as candidate genes, are members of the beta-growth factor (TGF-β) family, directly related to twinning, increasing ovulation, and growth and development of ovarian follicles in sheep and goats (Pramod et al., 2013). The FOXL2 gene is involved in ovarian growth and function, as well as early stages of mammalian ovarian growth (Baron et al., 2005). FSHB (Follicle Stimulating Hormone Subunit Beta) is a critical gene in follicle-stimulating hormone activity and peptide hormone metabolism. Moreover, variations in this gene may affect signaling of follicle differentiation and ovulation (Zi et al., 2020). The GH (Growth Hormone) gene is directly involved in nutrition-induced changes in the control of reproductive functions, e.g., ovarian follicular growth and development, cell division, and ovulation (Zhang et al., 2011). The PRLR (prolactin receptor) gene has been detected in various tissues, including brain, ovary, placenta, and uterus in various mammals, especially small ruminants. This hormone is involved in many endocrine activities and is essential for reproductive function (Ozmen et al., 2011). Therefore, according to the functions of hub genes identified in the circRNA–lincRNA–lncRNA–miRNA–mRNA ceRNA regulatory network and modules, we concluded that these genes have important roles in reproductive performance and fertility of goats. In that regard, they are involved in endocrine glands, growth, cell differentiation, as well as follicle maturation, and increased ovulation and could be selected in breeding programs to increase economic benefits. In addition, most genes involved in the ceRNA regulatory network encode signaling pathways regulating pluripotency of stem cells, cytokine-cytokine receptor interaction, ovarian steroidogenesis, and neuroactive ligand-receptor interaction, thereby confirming functions of the involved genes, especially hub genes. Signaling and metabolic pathways encoded by the genes involved in the circRNA–lincRNA–lncRNA–miRNA–mRNA ceRNA regulatory network and modules are presented in Figures 2 and 3. The signaling pathways of TGF-beta, regulating pluripotency of stem cells, Hippo, MAPK, PI3K-Akt, and FoxO are encoded by ceRNA regulatory network and modules. The TGF-beta signaling pathway encodes a large group of related structural proteins, including bone morphogenetic proteins, growth factor, and differentiation (Liu et al., 2018). Signaling pathways regulating pluripotency of stem cells are encoded by pluripotent stem cells (PSCs), which have potential to produce all 3 germ cell layers. It is noteworthy that embryonic stem cells (ESCs) are derived from the inner cell mass (ICM) of blastocyst-stage embryos (Mossahebi-Mohammadi et al., 2020). The Hippo signaling pathway is involved in inhibiting cell proliferation and enhancing apoptosis (Saucedo and Edgar, 2007). MAPK3 (Mitogen-Activated Protein Kinase 3) is a gene in a MAP kinase family. This gene has a major role in the signaling cascade that regulates various cellular processes such as proliferation, differentiation, and cell cycle progression in response to a variety of extracellular signals (Miao et al., 2016a). Functional metabolic pathways such as cytokine-cytokine receptor interaction, ovarian steroidogenesis, neuroactive ligand-receptor interaction, and growth hormone synthesis secretion and action are also encoded in these subnetworks, with important roles in reproduction and fertility. In this study, a computational approach with a circRNA–lincRNA–lncRNA–miRNA–mRNA ceRNA regulatory network was performed using predicted and validated expression profiles of RNAs. Spatiotemporal differential expression in various tissues, especially ovarian follicles, covered the potential roles of types of RNAs in the transcriptional and post-transcriptional regulation of genes involved in fertility. A common explanation for inconsistencies in our results were differences in applied molecular techniques (GWAS, Microarray, RNA-Seq, and simple relative gene expression), differences in ovarian tissue, time of sampling, and bioinformatics algorithms. Limitations of the present studies include the lack of a single, comprehensive dataset with similar environmental conditions and ovarian tissue from similar goat breeds. In summary, we concluded that identified transcriptomic signatures are potentially important biomarkers to better understand functional pathways involved in fertility in female goats. Further efforts are needed to elucidate the specific biological functional types of RNAs in reproduction and fertility. Moreover, our findings integrated circRNAs, lincRNAs, lncRNAs, miRNAs, and mRNAs based on an integrated approach from bioinformatics analyses and literature mining to construct ceRNA regulatory networks. Although this can be considered a robust approach to detect greater insights into biological processes, further research will be needed to confirm our results. 4 Conclusions This study used a novel approach to combining various types of RNA as an integrated network in goat fertility. Analyses of single-cell RNA sequencing data resulted in identification of 150 DEGs in goats with high versus low fertility. Among them, 80 genes were up-regulated and 70 were down-regulated. Moreover, 81 mRNAs/genes, 58 circRNAs, 8 lincRNAs, 19 lncRNAs, and 55 miRNAs, all well-known types of regulatory RNAs, were obtained from literature mining. Using circRNA–lincRNA–lncRNA–miRNA–mRNA ceRNA regulatory network was constructed and these identified RNAs were mainly associated with transcriptional regulatory activities and signalling receptor binding activities in terms of molecular functions, as well as reproductive functions such as ovulation cycle, ovarian follicle development, growth, and differentiation cells based on biological processes. Furthermore, our results are a valuable resource to elucidate molecular networks and the functions of DEGs underlying ovarian follicular development and they increase understanding of the genetic basis of high- versus low-fertility goats.