Ibrar Muhammad Khan1†
Hongyu Liu1*†
Jingyi Zhuang1Nazir Muhammad Khan2Dandan Zhang1Jingmeng Chen1Tengteng Xu1Lourdes Felicidad Córdova Avalos1Xinqi Zhou1
Yunhai Zhang1*- 1Anhui Provincial Laboratory of Local Livestock and Poultry Genetical Resource Conservation and Breeding, College of Animal Science and Technology, Anhui Agricultural University, Hefei, China
- 2Department of Zoology, University of Science and Technology, Bannu, Pakistan
Wandong cattle are an autochthonous Chinese breed used extensively for beef production. The breed tolerates extreme weather conditions and raw feed and is resistant to tick-borne diseases. However, the genetic basis of testis development and sperm production as well as breeding management is not well established in local cattle. Therefore, improving the reproductive efficiency of bulls via genetic selection is crucial as a single bull can breed thousands of cows through artificial insemination (AI). Testis development and spermatogenesis are regulated by hundreds of genes and transcriptomes. However, circular RNAs (circRNAs) are the key players in many biological developmental processes that have not been methodically described and compared between immature and mature stages in Bovine testes. In this study, we performed total RNA-seq and comprehensively analyzed the circRNA expression profiling of the testis samples of six bulls at 3 years and 3 months of developmental age. In total, 17,013 circRNAs were identified, of which 681 circRNAs (p-adjust < 0.05) were differentially expressed (DE). Among these DE circRNAs, 579 were upregulated and 103 were downregulated in calf and bull testes. The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses revealed that the identified target genes were classified into three broad functional categories, including biological process, cellular component, and molecular function, and were enriched in the lysine degradation, cell cycle, and cell adhesion molecule pathways. The binding interactions between DE circRNAs and microRNAs (miRNAs) were subsequently constructed using bioinformatics approaches. The source genes ATM, CCNA1, GSK3B, KMT2C, KMT2E, NSD2, SUCLG2, QKI, HOMER1, and SNAP91 were found to be actively associated with bull sexual maturity and spermatogenesis. In addition, a real-time quantitative polymerase chain reaction (RT-qPCR) analysis showed a strong correlation with the sequencing data. Moreover, the developed model of Bovine testes in the current study provides a suitable framework for understanding the mechanism of circRNAs in the development of testes and spermatogenesis.
Introduction
The local Wandong cattle are viewed as one of the primary breeds in Anhui, China, and moreover, they are considered as one of the most economically important breeds in the beef industry. However, their development and utilization are still in the initial stages (Zhou et al., 2001). The reproductive efficiency of breeding bulls is an important consideration in livestock production systems, and the genetic mainstream is obtained through the selective application of germ cells from the desirable sire (Waqas et al., 2019). Hence, molecular strategies are required to improve the reproduction traits and to explore the genetic potential for economically important traits in local cattle breeds.
Mammalian spermatogenesis can be divided into three distinct stages: i) self-renewal and mitosis in spermatogonia to form spermatocytes, ii) meiosis in spermatocytes to generate elongated haploid spermatids, and iii) postmeiotic modifications in haploid spermatids to form spermatozoa (Hecht, 1998; Griswold, 2016; Ibtisham et al., 2017). The crucial stages of spermatogenesis are regulated at the transcriptional, post-transcriptional, and epigenetic levels by a well-coordinated genomics network (Eddy and O'Brien, 1997; Yu et al., 2003). Male gonads consist of complicated transcriptomic elements, with more than 15,000 differentially expressed genes (DEGs) involved in spermatogonial maturation (Ramsköld et al., 2009; Soumillon et al., 2013). Recent studies have demonstrated that noncoding RNAs (nc-RNAs) play an important role at the post-transcriptional level during spermatogenesis (Mukherjee et al., 2014; Robles et al., 2017; Gao et al., 2019).
Circular RNAs (circRNAs) are nc-RNAs with a circular closed loop and play a vital regulatory role in many biological processes. The circRNAs are resistant to RNase degradation because of the lack of 5′- and 3′-ends (Sanger et al., 1976; Lasda and Parker, 2014). In the modern era of bioinformatics and sequencing technology, a significantly high number of circRNAs have been revealed in the mouse brain and testis. It has been shown that the testis consists of the second highest number of tissue-specific circRNA-level genes (You et al., 2015). A total of 15,996 circRNAs have been identified in humans, of which 10,792 (67%) have been observed in the testis (Guo et al., 2014; Dong et al., 2016). Remarkably, a primary report has shown that 1,017 circRNA-linked genes are mostly related to spermatogenesis, sperm motility, and fertilization (Dong et al., 2016). However, little is known about the global profile and characteristics of circRNAs in mammalian testis development and spermatogenesis.
In the present study, we provided a comprehensive catalog of circRNAs in the testes of bulls at two different developmental stages and detailed new insights into the functional activities of circRNAs in testicular development and spermatogenesis. The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomics (KEGG) analyses were performed to ascertain the molecular mechanisms of reproduction that were modulated during growth and development in bulls.
Materials and Methods
Animals and Sample Collection
Three physically healthy individuals (n = 3) of two different age groups, andrologically mature bulls (3 ± 0.014 years) and immature calves (3 ± 0.24 months), were selected for the current study (Tomlinson et al., 2017). The native cattle are also known as Wandong in the Anhui province of China. All testicular samples were obtained immediately after the slaughtering process under the veterinary surgical protocol in Fengyang County, Anhui. All pairs of testes were processed by incising the scrotum medially and uncovering the right and left testicles within the tunica vaginalis (Lunstra and Echternkamp, 1982). After removing the fascia tissues, tunica vaginalis, and tunica albuginea from the testis, three slices of testicular samples were cross-cut from the middle aspect of the testis through fine-scale dissection. One cross-cut section was stored in 4% formaldehyde solution for histological assessment, and the others were immediately immersed in liquid nitrogen and stored until total RNA extraction (Wu et al., 2020).
Histological Assessment
Testicular tissue samples that had been preserved in 4% formaldehyde (Wuhan Servicebio Technology Co., Ltd.) for 72 h were used for histological sections (Fischer et al., 2008). The tissues were sectioned into 6-μm-thick sections and stained with Masson’s trichrome stain. The histomorphology of testicular tissues was assessed using different magnification powers.
RNA Isolation, Quantification, and Qualification
The testicular samples were sent to a commercial sequencing facility (Beijing Novogene Corporation, China) for Seq-analysis. Total RNA was isolated from the testis samples using the Trizol reagent (Takara, Beijing, China) under sterile conditions and treated with the DNase I enzyme to remove endogenous DNA contamination according to the manufacturer’s protocol. RNA quality and contamination were assessed by 1% agarose gel electrophoresis. The purity was checked using a NanoPhotometer® spectrophotometer (IMPLEN, CA, United States). The total RNA concentration was assessed using a Qubit® RNA Assay Kit in a Qubit® 2.0 fluorometer (Life Technologies, CA, United States). RNA integrity was measured using an RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, United States).
cDNA Library Construction for circRNA Sequencing Analysis
A total of 3 μg of RNA from each testis sample was used to construct the cDNA libraries using the NEB-Next® Ultra-TM RNA Library Prep Kit for Illumina® (NEB, United States). Novogene constructed a chain-specific library by removing ribosomal RNA (Parkhomchuk et al., 2009). First, the ribosomal RNA was removed from the total RNA, and then, the RNA was divided into many short slices of 250–300 bp. The first cDNA strand manufactured from the fragmented RNA was used as the template strand, and random oligonucleotides were used as primers. Subsequently, the RNA strand was degraded with RNase H, and then, the second cDNA strand was manufactured along with dNTPs (dATP, dGTP, dUTP, and dCTP) using the DNA polymerase I system. After purification, the double-stranded cDNA was fixed, followed by the addition of a poly A-tail and sequencing adaptors. A cDNA of approximately 200 bp was separated using the AMPure XP beads. Finally, the USER enzyme was used to degrade the second (uracil-containing) strand of the cDNA, and PCR amplification was performed to obtain the library.
Alignment of RNA-Sequencing Reads and circRNA Identification
To ensure the quality and reliability of the data analysis, it is necessary to filter the original raw data. The low-quality reads (exceeding 50% low-quality bases, e.g., where Q-phred score ≤ 20), reads exceeding 10% possessing (N) nucleotides (where N, the proportion of reads that cannot be determined, is greater than 0.002), and reads with seq-adaptors were removed to obtain the clean reads, which were aligned to the ribosomal RNA database using the Bowtie2 tool (Langmead and Salzberg, 2012). The ribosomal RNA-free reads of each calf and bull sample were charted to the reference genome (Bos taurus-UMD 3.1.1) by TopHat2 (v. 2.0.3.12) as described previously (Kim et al., 2013). The charted reads of each sample were assembled using String-Tie v1.3.1 (Pertea et al., 2016). Subsequently, 20-mer sequences were removed from both ends of the unmapped reads and aligned according to the reference genomes to identify the exclusive anchor positions within the splicing position. The anchor reads were aligned in the reverse head-to-tail direction that exhibited back splicing (Hansen et al., 2013) and led to key identification of circRNAs (Figure 1). The identified circRNAs were blast for the specific functions against circBase, a database for circRNAs (http://www.circbase.org/) (Glažar et al., 2014), and assigned to different features, such as genomic classification, length distribution, and chromosomal distribution. The circRNAs that did not find their annotations were considered novel.
FIGURE 1. Scheme showing the difference between linear and back-splicing and also illustrating the RNA reads’ accommodation to each other that makes the stable circular RNAs.
Putative Identification and Coding Potential of circRNAs
Putative circRNAs were identified by filtration of anonymous transcripts. To reduce the rates of false positives, the assembled transcripts were filtrated to gain putative circRNAs by the following steps: i) the transcripts having more than one exon were selected, ii) transcripts having more than 200 bp in length and free from exon region overlapping were desired, and iii) transcripts with high expression levels where the FPKM value was more than 2 log2 (fold change) were selected. All known transcripts of the database were analyzed using the Cuffcompare software. The coding potential of circRNAs was determined using the software programs Coding–Non-Coding Index (CNCI) (Sun et al., 2013), Pfam-scan (Bateman et al., 2004; Finn et al., 2014), and coding potential calculator (CPC) (Kong et al., 2007). All transcripts identified with coding potential were refined, and those without coding potential were recognized as putative circRNAs.
Expression Pattern of Differentially Expressed circRNAs and miRNA Endogenous Sponge
The circRNA gene expression level was assessed using the fragments per kilobase of transcript sequence per million base pairs sequenced (FPKM) value. The significance of DE circRNAs at the gene or transcript level was analyzed, and functional genes relevant to the developmental groups were identified. Cuffdiff (v2.1.1) software was used to project the TPM levels of circRNAs (Trapnell et al., 2010). An edgeR kit (http://www.rproject.org/) was used to classify differentially expressed circRNAs through samples. CircRNAs with a fold change greater than or equal to 2 and a p-value < 0.05 in comparison between the samples are differentially expressed. We identified the endogenous sponge interaction between the differentially expressed circRNA and miRNA by means of three software programs: miRanda (v. 3.3a), MIREAP, and TargetScan (v. 7.0).
Bioinformatics Analysis of GO and KEGG
The GO term analysis and classification provides significant source genes that are involved in different biological functions, and DE circRNA genes were employed by the GO-seq R package (Young et al., 2010). All source genes in the pathways and the background genes were charted to the GO terms in the database (http://www.geneontology.org/). The GO terms with the corrected p-value (p < 0.05) were recognized as relatively enriched by DE genes according to the definition of the hyper-geometric test. The significance analysis of the term enrichment was corrected by FDR, and the corrected p-value (Q-value) was obtained (Ashburner et al., 2000). The KEGG pathway-based study further elaborated the explanation of source genes and their biological functions (http://www.genome.jp). KOBAS software was used to test the enriched DE source genes in the KEGG pathways, as suggested by (Mao et al., 2005). The enrichment analysis significantly identified the signaling and metabolic pathways, to which circRNA source genes and total background genes were charted.
Validation of the circRNAs Gene via RT-qPCR
Whole RNA was extracted from immature and mature groups using the Trizol reagent (Life Technologies, 182805, United States), and its concentration was measured using a Nanodrop spectrophotometer. The good-quality RNAs were then converted into cDNAs using a QuantiTect Reverse Transcription Kit (Qiagen, 205311, Germany) in accordance with the manufacturer’s protocol. The software programs Primer 3 web version 4.0.0 and basic local alignment search tool (BLAST; https://blast.ncbi.nlm.nih) were used to design gene-specific primers. The obtained PCR products were detected by 3% agarose gel and Sanger sequencing (Sangon Biotech, Shanghai, China). The divergent primers were used in this study and are presented in Supplementary File 1. The expression of circRNAs was detected by using a StepOnePlus Real-Time PCR System (Applied Biosystems, United States) using the FastStart SYBR Green Master Mix (Roche, Germany). Each PCR comprised 7.5 μl 2× SYBR Green PCR Master Mix, 1.5 μl cDNA, 0.5 μl each of reverse and forward primers, 4.7 μl nuclease-free water, and 10.3 μl ROX dye. The following q-PCR thermal cycles were carried out: i) predenaturation at 95°C for 10 min, ii) followed by 45 denaturation amplification cycles at 95°C for 15 s, iii) annealing at 60°C for 10 s, and iv) an extension cycle at 72°C for 20 s. The quantification cycle (Cq) of each target gene was normalized to that of the reference GAPDH gene. The q-PCR experiment was performed in triplicate to minimize the risk of error in the experiments. The Cq values were obtained and transferred to a Microsoft Excel sheet for further relative quantification analysis using the 2−ΔΔCt method (Sherman and Lempicki, 2009).
Statistical Analysis
Data on circRNA transcripts of interest from immature and mature testicular tissues were analyzed using the student t-test (SPSS 17.0) and presented as mean ± SEM. Before conducting the t-test, the data distribution and variances between the two groups were evaluated and found to be normal and homogeneous, respectively. Mean values were believed to be significantly different with *p < 0.05 and **p < 0.01.
Results
Histomorphology of Testes
A histological study of bull (B. taurus) testes showed a noteworthy difference between immature and mature stages. Under 100X microscopic examination, the diameter of seminiferous tubules was much smaller in immature than in mature testes. At the same time, similar conditions were observed in the interstitial connective tissue of immature and mature testes, as highlighted in Figures 2A,B. Under 400X microscopic examination, we obtained further details and observed more developmental stages of biological spermatozoa in mature than in immature testes. Sertoli cells were significantly increased in number and were found very near the basement membrane of seminiferous tubules in mature testes, as shown in Figures 2E,F.
FIGURE 2. Histomorphological analysis of bull testicular tissues at ages of 3 months and 3 years was performed under a microscope at ×100, × 200, and × 400 magnifications. The segments (A,C,E) represent the morphology of 3-month-old testes, whereas (B,D,F) represents the morphology of 3-year-old bulls. The red arrows indicate different cell types. 1: spermatogonia, 2: Leydig cells, 3: Sertoli cells, 4: spermatocytes, 5: round spermatids, 6: elongated spermatids, and 7: sperms.
Transcriptome Sequencing of circRNAs in Immature and Mature Bull Testes
A total of six experimental animals, three immature (3 ± 0.014 months) and three mature (3 ± 0.24 years), were selected from a local cattle station, and their testis samples were collected. To investigate the circRNA profiles of the immature and mature testis tissues, we prepared the RNA-seq libraries of all tissue samples. The raw reads of sequencing data were acquired using an Illumina HiSeq 4000 platform (Illumina, Inc, San Diego, CA, United States). Raw reads were filtered and classified into different segments, including clean reads (70,333,218, 98.95%), N-containing reads (38,769, 0.05%), low-quality reads (126,760, 0.18%), and adopter-related reads (584,096, 0.82%) (Figure 3). The raw reads for each sample ranged between 118,785,372 and 142,165,686, while the clean reads per sample ranged between 116,004,150 and 140,666,436. Thus, the raw and clean reads together yielded 90 GB of data, and the GC content ranged from 47.09 to 53.57% (Supplementary File 2). All Q20 values of the subjected reads in the six samples exceeded 97%, as shown in Supplementary File 3. More than 95.5% of the clean reads were coordinated to B. taurus UMD3.1.
FIGURE 3. Pie charts show the raw reads of sequencing data from immature and mature bull testis samples that were subjected to RNA sequencing analysis.
Correlation and Differential Analysis Between the Testis Samples
The testis samples were collected from three calves (Calf1, Calf2, and Calf3) and three bulls (Bull1, Bull2, and Bull3). The Pearson’s correlation heatmap was constructed according to the expression profiles for samples’ quality control. All samples in the calf and bull groups showed the linear correlation (Figure 4A). The differences of transcripts of two groups demonstrated in the sample are shown in the PCA 3D map (Figure 4B). We also discovered that bull and calf groups have noticeable differences, whereas the testis samples of calves and bulls clustered separately and showed the difference.
FIGURE 4. Pearson correlation and principle component analysis were constructed between groups and samples. (A) Correlation among the sequencing samples of calf and bull testes. (B) PCA 3D plot among the sequencing samples.
Identification and Characterization of circRNAs
A total of 17,013 putative circRNAs were recognized with at least one read spanning and head-to-tail splicing in immature and mature testes. Based on their genomic location and features, circRNAs were divided into three subclasses, denoted as circ-exon (79.0%), circ-intergenic (14.0%), and circ-intron (7.0%) (Figure 5A). The total length plot quantile of circRNAs (approximately 75%) was not more than 1,000 nt, and the average length was 400 nt, as shown in Figure 5B. According to their host gene site, all the circRNAs were widely distributed on all sets of chromosomes, while none of the circRNA was mapped to the mitochondrial genome. The total available set of chromosomes generated 100 circRNAs in both developmental stages of the Bovine testis (Figure 5C). The coding potential of circRNAs was analyzed using the software programs CPC2, CNCI, and PFAM, and 6,011 coding circRNAs were identified in all samples (Figure 5D). Detailed information on the top 40 up- and downcoding circRNAs is listed in Table 1.
FIGURE 5. General features of circRNAs in immature and mature bull testes. (A) Distribution of circRNAs in the genomic region and high magnitude identified from exon, intergenic, and intron parts of the genome. (B) Length distribution of circRNAs, while colors represent the different lengths of circRNAs in calf and bull testis samples. (C) Distribution of testis circRNAs on the cattle chromosomal sets. (D) Venn diagram shows the coding potential results of circRNAs and intersection of coding tools such as coding–non-coding index (CNCI), PFAM, and coding potential calculator (CPC2).
TABLE 1. Top 40 regulated coding circular RNAs in immature and mature testicular tissues of Wandong cattle.
Differentially Expressed Gene and circRNAs in Testicular Tissues
After the quantification process, the expression pattern of circRNAs was identified using Cuffdiff and Ballgown tools. A total of 681 DE circRNAs at the gene level were detected in immature and mature stages. The expression levels of circRNAs in the testis samples were calculated, taking into account the parameter of significance, whereas the log2 fold change was considered higher than or equal to 2 and p-adjusted < 0.05. According to the significance criteria, we detected 578 upregulated and 102 downregulated circRNAs between bull and calf testis samples. These up and down highlights of DE genes are displayed in volcanic plots (Figures 6A,B), and the list of total and DE circRNAs is shown in Supplementary Files 4–6. We also analyzed the DE circRNA by the hierarchical clustering method, which is another way to display DE genes and to cluster all the genes with similar expression patterns that may have common functions in metabolic and signaling pathways. The gene clusters on the left side were formed because of similar expression patterns (fold change >2, p < 0.05); the columns presented calves and bulls, and the expression from blue to red was gradually upregulated (Figure 6C).
FIGURE 6. Expression patterns of circular RNAs (circRNAs) and protein-coding genes. (A) Volcanic plots represent log10 (p-value) vs log2 fold difference in circRNAs in the fragments per kilobase of transcript per million mapped reads (FPKM) values between immature and mature testes. (B) Total abundance of circRNAs and differentially expressed (DE) highlights between calf vs bull. (C) Heatmap of DE circRNAs among calves and bulls. Red indicates upregulated coding genes, and blue indicates downregulated gene products.
GO and Enriched KEGG Pathways Analysis
CircRNAs are closed-loop structures that are most commonly generated by back-splicing events between the exons of coding genes. In a limited way, circRNAs regulate the host gene and thus perform the desired function at the genetic level. Hence, the analysis of GO terms and circRNA source gene functions might provide new insights. Many GO classification terms were significantly enriched in the source genes of circRNAs. A list of 26,815 source genes significantly associated with DE circRNAs was submitted to the database for functional annotation and GO classification analysis, resulting in the identification of 50 significant GO terms, as shown in Figure 7. The functional annotation study revealed that several source genes actively participate in spermatogenesis, sperm morphology, developmental stages, and reproduction. The functional classification of source genes showed that their products were assigned to the three GO categories of biological processes, cellular components, and molecular functions (Supplementary File 7). To further understand the role of circRNA host genes in the developmental stages of the Bovine, a pathway analysis was applied using the KEGG pathway database and a total of 147 signaling pathways linked to circRNA gene products were identified (Supplementary File 8). The top 20 statistically enriched pathways (p < 0.05) were subjected to further analysis and were found to be significantly associated with signaling pathways such as lysine degradation, cell cycle, propanoate metabolism, adherens junction, and cell adhesion molecules (Figure 8). The host genes such as KMT2E, EHMT1, and NSD2 were mediated by the lysine degradation signaling pathways but actively participated in the reproductive traits of the Bovine. Similarly, ATM, GSK3B, and CCNA1 are host genes of the cell cycle signaling pathway and influence the biological descriptions relevant to spermatogenesis and testis development. The most significantly enriched pathways and DEGs are listed in Table 2.
FIGURE 7. GO enrichment analysis of differentially expressed (DE) circular RNA (circRNA) genes between immature and mature testicles was performed. The DE circRNA genes are divided into the following three functional categories: biological processes (BP), cellular components (CC), and molecular functions (MF), while the left and right y-axes show the percentage and numbers of circRNAs host genes, respectively.
FIGURE 8. KEGG enrichment pathway analysis was performed. (A) Top 20 enriched pathways of differentially expressed (DE) circRNA host genes in immature and mature stages of testicular growth. The color dots represent different signaling pathways that are regulated by the circRNA host genes. The size and different colors of dots show the numbers of genes and significance levels in enriched pathways. (B) Top 10 most relevant reproduction pathways representing the total background genes and significant source genes of circRNAs.
TABLE 2. Most enriched biological pathway terms and regulated DE genes in immature and mature testes.
Prediction of Deferentially Expressed circRNA and miRNA Endogenous Sponges
Some studies in recent years have suggested that miRNAs play a role in all aspects of spermatogenesis and testis formation; thus, miRNAs could be used as biomarkers for reproductive traits (Yadav and Kotaja, 2014). Reports also suggested that circRNAs act as miRNA sponges and play an important role in regulating gene expression in posttranscriptional processes (Kulcheski et al., 2016). We analyzed further to determine whether circRNAs found in immature and mature testes act as endogenous miRNA sponges. We noticed that a total of 4,300 circRNAs have potential binding miRNA sites, but other circRNAs were predicted to have no possible binding targets for miRNAs. Each circRNA may attach to one or more targeted miRNAs. As a result, the percentage of circRNAs consisting of different numbers of miRNA targets was further tested. The majority of circRNAs have at least two miRNA-binding sites; here, the proportion of circRNAs containing 6–10 miRNA targets is the largest, as illustrated in Figure 9A. The binding interactions between selected DE circRNAs in testicular tissues and their predicted miRNA targets are shown in Figure 9B. Recent studies have proposed that circRNAs act as miRNA sponges and play a critical role in regulating gene expression in pathways through complex (circRNAs–miRNAs genes) chains (Hansen et al., 2013; Memczak et al., 2013). Approximately 758 miRNAs have been predicted, and many of them, such as bta-miR-204, bta-miR-532, and bta-miR-34 groups, have been correlated with spermatogenesis (Zhang et al., 2015).
FIGURE 9. Differentially expressed (DE) circRNA and miRNA networks in the bull testis were identified. (A) During binding interaction, different numbers of miRNA targets were processed by circRNAs. (B) DE circRNAs and the predicted targeted miRNAs network were identified for calf and bull testes. Green nodes represent the targeted miRNAs, and yellow nodes show circRNAs.
Validation of DE circRNA-Seq Results by RT-qPCR
To justify the sequencing analysis of DE circRNAs through RT-qPCR, we selected the truly significant and exonic DE circRNAs for validation. We randomly selected nine circRNAs of different nucleotide richness values and lengths, and divergent primers were prepared to verify the expression levels of circRNAs in the testis. Ten host genes were proposed for circRNA detection, of which seven, including ATM serine/threonine kinase (ATM), cyclin A1 (CCNA1), glycogen synthase kinase 3 beta (GSK3B), lysine methyltransferase 2C (KMT2C), lysine methyltransferase 2E (KMT2E), nuclear receptor binding SET domain protein 2 (NSD2), and succinate-CoA ligase GDP-forming subunit beta (SUCLG2), were upregulated. The remaining three genes, including QKI, KH domain-containing RNA binding (QKI), homer scaffold protein 1 (HOMER1), and synaptosome-associated protein 91 (SNAP91), were downregulated. The RT-qPCR analysis validated the expression of all nine circRNAs, which were found to regulate the functions of the aforementioned genes, as shown in Figures 10A,B. The gel electrophoresis findings showed that the PCR products of circRNA were single bands, which designated the presence of specific circRNAs, Figure 10C. Additionally, the back-splice junction of circRNAs was highlighted by further Sanger sequencing of the PCR products, Figure 10D. These genes work in specific reproductive-related pathways, such as lysine degradation, cell adhesion molecules (CAMs), focal adhesion, cell cycles, and adherens junctions, thus supporting reproductive functions in animals. Therefore, the obtained results confirmed the accuracy of the ribo-depleted RNA-seq results.
FIGURE 10. Confirmation of the differentially expressed (DE) circRNAs was done by RT-qPCR at two different developmental stages of the Bovine testis. The figures (A,B) show the expression patterns of ten DE circRNAs and their host genes (ATM, CCNA1, GSK3B, KMT2C, KMT2E, NSD2, SUCLG2, QKI, HOMER1, and SNAP91) in the immature and mature groups. The data were measured by the 2−ΔΔCt method, and GAPDH was used as a housekeeping gene. The data are presented as the mean ± SEM, *p < 0.05, **p < 0.01. (C) PCR products of circRNA were confirmed via gel electrophoresis, where M represents marker I. (D) Back-splice junction of circRNA was confirmed by the Sanger sequencing analysis.
Discussion
Mammalian testicular growth and spermatogenesis are complex biological transformation processes that are regulated by a strong combination of coding and noncoding genes and transcriptomes (Card et al., 2013; Djureinovic et al., 2014). Testicular growth includes the growth of testicular tissues during the embryonic and postnatal phases (Svingen and Koopman, 2013). The bull testes first gradually expand until approximately 25 weeks of growth and then rapidly expand until adolescence at 37–50 weeks of age (Rawlings et al., 2008). To analyze the molecular biology of testis development and reproductive transformation in Bovine species, we selected n = 3 animals of two different age groups: 3-month-old (immature) calves and 3-year-old (adult) bulls of a local Wandong cattle breed. Gao et al. (2018) tested the similar experimental design in the Qinchuan cattle testis and selected n = 1 animals of two different age groups: neonatal (1 week) and mature bulls (4 years).
Previously, nc-RNAs were regarded as useless stocks or junk RNAs. However, various types of ncRNAs, such as miRNAs, PIWI-interacting RNAs (pi-RNAs), long noncoding RNAs (lnc-RNAs), and circRNAs (Lee et al., 1993; Girard et al., 2006; Memczak et al., 2013; Ulitsky and Bartel, 2013), have been identified and reported to play important roles in cellular and tissue development. In recent years, research has focused on the functional features of nc-RNAs, and much work has been reported specifically in animal reproduction (Svingen and Koopman, 2013). In addition, lnc-RNAs play a regulatory role in the testis and are expressed more strongly in mammalian testes than in other organs (Soumillon et al., 2013). While circRNAs are essential members of the nc-RNA family, they are also known to be involved in animal reproduction and disease control (Wang et al., 2016). Limited studies have characterized circRNAs in reproductive organs, such as placentae, embryos, and testes, and in cell types including oocytes, granulosa cells, immature sperm cells, seminal plasma cells, and mature sperm cells (Dong et al., 2016; Quan and Li, 2018; Chioccarelli et al., 2019). circRNAs were identified between the good and poor sperm subpopulations in terms of kinetic parameters and morphological characteristics. These comparative data highlighted 148 DE circRNAs between the two semen quality classes (Dong et al., 2016). Previous studies have reported that circRNAs are the primary regulators of multiple biological processes, but they have not been methodically explored in the Bovine testis during reproductive development. In this study, we identified and marked circRNAs in Bovine testes and discovered possible genes that determine the development of testes.
A total of 17,013 and 681 coding circRNAs were identified in the Bovine testes at two different immature and mature stages. You et al. (2015) reported that the second highest number of circRNAs was found in the testes compared to other body organs. The magnitude of circRNAs in testicular samples was significantly greater than that of circRNAs found in other organs, for example, chicken theca cells contained 14,502 circRNAs (Shen et al., 2020), pig pituitaries had 5,148 circRNAs (Ulitsky and Bartel, 2013; Chen et al., 2020), and buffalo fat adipose tissues had 5,141 circRNAs (Huang et al., 2019). We identified 681 coding circRNAs, of which 579 circRNAs were upregulated and 103 were downregulated in immature and mature bull testes. The study composed of Qinchuan cattle identified a total of 21,753 candidate circRNAs, where 2023 circRNAs were downregulated and 2,225 circRNAs were upregulated (Gao et al., 2018). A total of 2,326 differentially expressed circRNAs (DECs) were found in testicular development, of which 1,003 circRNAs were upregulated in adult boar testes and 1,323 circRNAs were downregulated; an analysis of transcriptomic changes in boar testicular tissues revealed an agreement with our results (Zhang et al., 2021). A study of transcriptomic changes in cattle-yak testes showed a tendency toward upregulation rather than downregulation and identified 679 upregulated and 2,281 downregulated genes (Cai et al., 2017). An RNA-seq analysis showed that a total of 10 ,095 genes were substantially differentially expressed in porcine testes, of which 5,199 were upregulated and 4,896 were downregulated in immature and mature porcine testes (Song et al., 2015). Approximately 242 genes were associated with porcine spermatogenesis. The key reason for testis transcriptome profiling is to discover the transcriptional factors that play crucial roles in testis development and spermatogenesis and thus clarify the genetic structure of the testis.
We performed GO annotation and KEGG pathway enrichment analysis to better describe the circRNAs involved in the biological processes of testis development and spermatogenesis. These processes may indicate the roles of circRNAs in the testis or reflect the cellular differences between immature and mature bull testes. In GO annotation, 60 host genes were assigned to spermatogenesis processes, 12 genes to testis development, and 23 host genes were found solely relevant to reproduction phenomena, such as behavior and primary sexual characteristics. Gao et al. (2018) executed the GO annotation, and 44 host genes were assigned to reproduction and the reproductive process, and 25 host genes were associated with spermatogenesis, including PIWIL1 and the spermatogenesis-associated protein 6 (SPATA6). We found different candidate genes by analyzing these GO terms, and some of these genes, such as ATM serine/threonine kinase (ATM), glycogen synthase kinase 3 beta (GSK3B), and cyclin A1 (CCNA1), were found to be significantly enriched in the cell cycle pathway and associated with spermatozoa differentiation. The ATM protein is associated with low-fertile buffalo bull spermatozoa (M. K et al., 2019); in mouse spermatocytes, the meiotic recombination is initiated by SPO11-induced double-strand breaks (DSBs), and ATM controls this process (Huang et al., 2019; Paiano et al., 2020). GSK3B showed significant expression changes across these pubertal stages in the hypothalamus of piglets and is also involved in crucial biological processes that regulate economically relevant traits associated with beef cattle fertility (Fonseca et al., 2018). The GO and KEGG pathway enrichment analyses showed that the DE circRNA genes were associated with pathways, such as the cell cycle, lysine degradation, and tight junction, whereas the top 20 most significant enriched pathways were evaluated in this study. A total of 47 significant signaling pathways in the Qinchuan cattle testis were found, and key host genes including PIWIL1, DPY19L2, SLC26A8, IFT81, SMC1B, IQCG, TTLL5, and ACVR2A belong to the tight junction, the adherens junction, the TGFβ signaling pathway, progesterone-mediated oocyte maturation, and oocyte meiosis (Gao et al., 2018).
The probable reasons of variation between the findings of the present study and the results of Gao et al. regarding candidate genes and signaling pathways are due to the breed variations. In the current study, we collected samples from Wandong, while Gao et al. used Qinchuan breed testicular tissues for RNA-seq analysis. A study found that many candidate genes were differentially expressed (DE) and also differential alternative spliced in three different cattle breeds, including the Holstein, Jersey, and Cholistani. It suggests that breed-specific gene expression occurs at the transcriptional level and posttranscriptional modifications such as splicing, 5’ capping, and poly A-tail addition. The differentially expressed genes are enriched in KEGG pathways including translation, electron transport, and immune responses (Huang et al., 2012). In the Bovine breeds, in Merino vs Poll Dorset, SLC35A5 and ITM2C were identified as key DEGs; both have been reported for their association with conception and embryonic development (Hodge et al., 2021). We found different candidate genes and signaling pathways from the findings of Gao et al. due to variation in age, where they analyzed the testicular tissue in days 7 and in 4 years of age. However, we analyzed the same tissue in the age of 3 months and 3 years. As gene expression is greatly influenced by age, especially genes related to developmental phases, a comparative study was conducted on candidate genes having functional roles in regulation of muscle development at various growth stages in goat. In this study, i) 2 vs 9 months groups were compared and intramuscular fat-linked genes, e.g., MYH13, IFIT1, METTL21C, SERPINE1, EGR1, and ESRRG, were found; ii) 2 vs 24 months groups were matched and significant DEGs RPS25, MYH13, COMP, LOC102186300, and NR4A2 were found; and iii) 9 vs 24 months goat muscles were transcriptionally matched and PARM1, LOC102177715, MRF1-like, ARID5B, and PFKFB3 were found (Lin et al., 2017). Therefore, the candidate genes and signaling pathways are different in these two different studies of similar nature.
In mammalian testes, the cell cycle, lysine degradation, and tight junction and adherens junction pathways and host genes function together to control testis growth, spermatogenesis, and primary sexual activity (Young et al., 2015). The gene expression profiling of testes showed that some lysin degradation pathway genes, including lysine methyltransferase 2E (KMT2E) and nuclear receptor-binding SET domain protein 2 (NSD2), are active during spermatogenesis. Zhang et al. (2015) reported that the KMT2E gene plays key roles in various biological processes, including cell cycle progression, adult hematopoiesis, and spermatogenesis. A vital question raised in developmental biology is that how pluripotent stem cells can give rise to the cells derived from the three germ layers. The NSD2 gene has a dual role in pluripotency exit and germ layer specification of embryonic stem cells (Tian et al., 2019).
Exonic circRNAs have been shown to act preferentially in post-transcriptional regulation in mice and humans (Bose and Ain, 2018). Similarly, we found that the majority of DE circRNAs identified in this study have putative target miRNA-binding sites, suggesting that most circRNAs are likely to function as miRNA sponges. circRNA–miRNA network analysis further revealed that circRNAs can interact with miRNAs in multiple directions at different rates in the testes, which is consistent with observations in previous studies (Liang et al., 2017). Therefore, when translated, some circRNAs can inhibit or relieve miRNA repression. In this transcriptomic network, we selected 11 DE circRNAs that were significantly matched with 160 miRNAs, wherein the bta-miR-34, bta-miR-532, and bta-miR-204 families were involved in cattle spermatogenesis (Tscherner et al., 2014; Zhang et al., 2015). The related circRNAs may function as miRNA sponges to control the growth of testes and spermatogenesis. The circRNAs act as miRNA sponges, and research on Qinchuan cattle showed that 758 miRNAs matched with significantly differentially expressed circRNAs which play an important role in regulating gene expression in bull testes (Gao et al., 2018). The RT-qPCR validation results indicated that the dynamic expression of the circRNA had a strong interaction with the RNA-seq data. Previous studies have shown that circRNAs exhibit significant developmental stage-precision activity in human testes, suggesting their potential contribution to spermatogenesis (Dong et al., 2016). In rodents, circRNA patterns have been analyzed in three distinct stages, mitotic, meiotic, and postmeiotic germ cells (Lin et al., 2016), and correlated with testis development (Zhou et al., 2018).
Conclusion
In conclusion, our study presented an extensive analysis of circRNA expression profiles in immature and mature Wandong cattle bull testes. The landscapes of circRNA expression in testes were detected and used for interpreting their regulatory structure in testis development and spermatogenesis in cattle bulls. To date, there is a lack of the genetical data to establish baseline factors related to reproductive performance of local cattle bulls in China. Nevertheless, the insights from this study could lead to a breakthrough in the reproductive efficiency of local cattle resources.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE183347, https://www.ncbi.nlm.nih.gov/geo/info/linking.html.
Ethics Statement
The animal study was reviewed and approved by the ethical committee of the faculty of veterinary, permit number SYXK 2016–007. Every effort was made to reduce the number of animals used and minimize the sufferings of the animals.
Author Contributions
IK and HL contributed to conceptualization. JZ contibuted to the methodology. IK contributed to software. DZ, TX, and XZ contributed to validation. JC contributed to formal analysis. HL contributed to investigation. YZ contributed to resources. HL contributed to data curation. IK contributed to the original draft. YZ contributed to final draft and editing. NK and LA contributed in visualization. HL contributed in supervision. YZ contributed to project administration. YZ contributed in funding acquisition.
Funding
This research was funded by the National Research and Development Program of China (No. 2018YFD0501700), the National Natural Science Foundation of China (No. 31101696), and Anhui Province Modern Agriculture Industry (Cattle, Goat, and Sheep) Technology System (No. AHCYJSTX-07).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.
Acknowledgments
We highly acknowledge Fengyang County Daming Agriculture and Animal Husbandry Technology Development Co., Ltd., who provided the animals and physical support.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.685541/full#supplementary-material
Abbreviations
GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; DE, differentially expressed; RT-qPCR, real-time quantitative polymerase chain reaction; FPKM, fragments per kilobase of transcript sequence per millions base pairs sequenced; CNCI, coding–non-coding index; CPC, coding potential calculator; nc-RNAs, noncoding RNAs; cds, coding sequence; BP, biological process; CC, cellular components; MF, molecular function; dNTP, deoxyribonucleotide triphosphate; dATP, deoxyadenosine triphosphate; dGTP, deoxyguanosine triphosphate; dUTP, deoxyuridine triphosphate; dCTP, deoxycytidine triphosphate.
References
Ashburner, M., Ball, C. A., Blake, J. A., Botstein, D., Butler, H., Cherry, J. M., et al. (2000). Gene Ontology: Tool for the Unification of Biology. Nat. Genet. 25, 25–29. doi:10.1038/75556
Bateman, A., Coin, L., Durbin, R., Finn, R. D., Hollich, V., Griffiths-Jones, S., et al. (2004). The Pfam Protein Families Database. Nucleic Acids Res. 32, D138–D141. doi:10.1093/nar/gkh121
Bose, R., and Ain, R. (2018). Regulation of Transcription by Circular RNAs. Circular RNAs. 1087, 81–94. doi:10.1007/978-981-13-1426-1_7
Cai, X., Yu, S., Mipam, T., Yang, F., Zhao, W., Liu, W., et al. (2017). Comparative Analysis of Testis Transcriptomes Associated with Male Infertility in Cattleyak. Theriogenology 88, 28–42. doi:10.1016/j.theriogenology.2016.09.047
Card, C. J., Anderson, E. J., Zamberlan, S., Krieger, K. E., Kaproth, M., and Sartini, B. L. (2013). Cryopreserved Bovine Spermatozoal Transcript Profile as Revealed by High-Throughput Ribonucleic Acid Sequencing. Biol. Reprod. 88, 49. doi:10.1095/biolreprod.112.103788
Chen, Z., Pan, X., Kong, Y., Jiang, Y., Zhong, Y., Zhang, H., et al. (2020). Pituitary-Derived Circular RNAs Expression and Regulatory Network Prediction during the Onset of Puberty in Landrace × Yorkshire Crossbred Pigs. Front. Genet. 11, 135. doi:10.3389/fgene.2020.00135
Chioccarelli, T., Manfrevola, F., Ferraro, B., Sellitto, C., Cobellis, G., Migliaccio, M., et al. (2019). Expression Patterns of Circular RNAs in High Quality and Poor Quality Human Spermatozoa. Front. Endocrinol. 10, 435. doi:10.3389/fendo.2019.00435
Djureinovic, D., Fagerberg, L., Hallström, B., Danielsson, A., Lindskog, C., Uhlén, M., et al. (2014). The Human Testis-specific Proteome Defined by Transcriptomics and Antibody-Based Profiling. Mol. Hum. Reprod. 20, 476–488. doi:10.1093/molehr/gau018
Dong, W. W., Li, H. M., Qing, X. R., Huang, D. H., and Li, H. G. (2016). Identification and Characterization of Human Testis Derived Circular RNAs and Their Existence in Seminal Plasma. Sci. Rep. 6, 39080–39111. doi:10.1038/srep39080
Eddy, E., and Obrien, D. (1997). 5 Gene Expression during Mammalian Meiosis. Curr. Top. Dev. Biol. 37, 141–200. doi:10.1016/s0070-2153(08)60174-x
Finn, R. D., Bateman, A., Clements, J., Coggill, P., Eberhardt, R. Y., Eddy, S. R., et al. (2014). Pfam: the Protein Families Database. Nucl. Acids Res. 42, D222–D230. doi:10.1093/nar/gkt1223
Fischer, A. H., Jacobson, K. A., Rose, J., and Zeller, R. (2008). Hematoxylin and Eosin Staining of Tissue and Cell Sections. Cold Spring harbor Protoc. 2008, pdb.prot4986. doi:10.1101/pdb.prot4986
Fonseca, P. A. d. S., Id-Lahoucine, S., Reverter, A., Medrano, J. F., Fortes, M. S., Casellas, J., et al. (2018). Combining Multi-OMICs Information to Identify Key-Regulator Genes for Pleiotropic Effect on Fertility and Production Traits in Beef Cattle. PLoS One 13, e0205295. doi:10.1371/journal.pone.0205295
Gao, Y., Li, S., Lai, Z., Zhou, Z., Wu, F., Huang, Y., et al. (2019). Analysis of Long Non-coding RNA and mRNA Expression Profiling in Immature and Mature Bovine (Bos taurus) Testes. Front. Genet. 10, 646. doi:10.3389/fgene.2019.00646
Gao, Y., Wu, M., Fan, Y., Li, S., Lai, Z., Huang, Y., et al. (2018). Identification and Characterization of Circular RNAs in Qinchuan Cattle Testis. Royal society open science. 5, 180413. doi:10.1098/rsos.180413
Girard, A., Sachidanandam, R., Hannon, G. J., and Carmell, M. A. (2006). A Germline-specific Class of Small RNAs Binds Mammalian Piwi Proteins. Nature 442, 199–202. doi:10.1038/nature04917
Glažar, P., Papavasileiou, P., and Rajewsky, N. (2014). circBase: a Database for Circular RNAs. Rna 20, 1666–1670. doi:10.1261/rna.043687.113
Griswold, M. D. (2016). Spermatogenesis: The Commitment to Meiosis. Physiol. Rev. 96, 1–17. doi:10.1152/physrev.00013.2015
Guo, J. U., Agarwal, V., Guo, H., and Bartel, D. P. (2014). Expanded Identification and Characterization of Mammalian Circular RNAs. Genome Biol. 15, 409. doi:10.1186/s13059-014-0409-z
Hansen, T. B., Jensen, T. I., Clausen, B. H., Bramsen, J. B., Finsen, B., Damgaard, C. K., et al. (2013). Natural RNA Circles Function as Efficient microRNA Sponges. Nature 495, 384–388. doi:10.1038/nature11993
Hecht, N. B. (1998). Molecular Mechanisms of Male Germ Cell Differentiation. Bioessays 20, 555–561. doi:10.1002/(sici)1521-1878(199807)20:7<555:aid-bies6>3.0.co;2-j
Hodge, M. J., de Las Heras-Saldana, S., Rindfleish, S. J., Stephen, C. P., and Pant, S. D. (2021). Characterization of Breed Specific Differences in Spermatozoal Transcriptomes of Sheep in Australia. Genes 12 (2), 203. doi:10.3390/genes12020203
Huang, J., Zhao, J., Zheng, Q., Wang, S., Wei, X., Li, F., et al. (2019). Characterization of Circular RNAs in Chinese Buffalo (Bubalus Bubalis) Adipose Tissue: a Focus on Circular RNAs Involved in Fat Deposition. Animals 9, 403. doi:10.3390/ani9070403
Huang, W., Nadeem, A., Zhang, B., Babar, M., Soller, M., and Khatib, H. (2012). Characterization and Comparison of the Leukocyte Transcriptomes of Three Cattle Breeds. PLoS ONE 7 (1), e30244. doi:10.1371/journal.pone.0030244
Ibtisham, F., Wu, J., Xiao, M., An, L., Banker, Z., Nawab, A., et al. (2017). Progress and Future prospect of In Vitro Spermatogenesis. Oncotarget 8, 66709–66727. doi:10.18632/oncotarget.19640
Kim, D., Pertea, G., Trapnell, C., Pimentel, H., Kelley, R., and Salzberg, S. L. (2013). TopHat2: Accurate Alignment of Transcriptomes in the Presence of Insertions, Deletions and Gene Fusions. Genome Biol. 14, R36–R13. doi:10.1186/gb-2013-14-4-r36
Kong, L., Zhang, Y., Ye, Z.-Q., Liu, X.-Q., Zhao, S.-Q., Wei, L., et al. (2007). CPC: Assess the Protein-Coding Potential of Transcripts Using Sequence Features and Support Vector Machine. Nucleic Acids Res. 35, W345–W349. doi:10.1093/nar/gkm391
Kulcheski, F. R., Christoff, A. P., and Margis, R. (2016). Circular RNAs Are miRNA Sponges and Can Be Used as a New Class of Biomarker. J. Biotechnol. 238, 42–51. doi:10.1016/j.jbiotec.2016.09.011
Langmead, B., and Salzberg, S. L. (2012). Fast Gapped-Read Alignment with Bowtie 2. Nat. Methods 9, 357–359. doi:10.1038/nmeth.1923
Lasda, E., and Parker, R. (2014). Circular RNAs: Diversity of Form and Function. Rna. 20, 1829–1842. doi:10.1261/rna.047126.114
Lee, R. C., Feinbaum, R. L., and Ambros, V. (1993). The C. elegans Heterochronic Gene Lin-4 Encodes Small RNAs with Antisense Complementarity to Lin-14. cell 75, 843–854. doi:10.1016/0092-8674(93)90529-y
Liang, G., Yang, Y., Niu, G., Tang, Z., and Li, K. (2017). Genome-wide Profiling of Sus scrofa Circular RNAs across Nine Organs and Three Developmental Stages. DNA Res. 24, 523–535. doi:10.1093/dnares/dsx022
Lin, X., Han, M., Cheng, L., Chen, J., Zhang, Z., Shen, T., et al. (2016). Expression Dynamics, Relationships, and Transcriptional Regulations of Diverse Transcripts in Mouse Spermatogenic Cells. RNA Biol. 13, 1011–1024. doi:10.1080/15476286.2016.1218588
Lin, Y., Zhu, J., Wang, Y., Li, Q., and Lin, S. (2017). Identification of Differentially Expressed Genes through RNA Sequencing in Goats (Capra hircus) at Different Postnatal Stages. PloS one 12 (8), e0182602. doi:10.1371/journal.pone.0182602
Lunstra, D.-D., and Echternkamp, S. E. (1982). Puberty in Beef Bulls: Acrosome Morphology and Semen Quality in Bulls of Different Breeds1. J. Anim. Sci. 55, 638–648. doi:10.2527/jas1982.553638x
Mao, X., Cai, T., Olyarchuk, J. G., and Wei, L. (2005). Automated Genome Annotation and Pathway Identification Using the KEGG Orthology (KO) as a Controlled Vocabulary. Bioinformatics 21, 3787–3793. doi:10.1093/bioinformatics/bti430
Memczak, S., Jens, M., Elefsinioti, A., Torti, F., Krueger, J., Rybak, A., et al. (2013). Circular RNAs Are a Large Class of Animal RNAs with Regulatory Potency. Nature 495, 333–338. doi:10.1038/nature11928
M. K, M. A., Kumaresan, A., Yadav, S., Mohanty, T. K., and Datta, T. K. (2019). Comparative Proteomic Analysis of High- and Low-fertile buffalo Bull Spermatozoa for Identification of Fertility-Associated Proteins. Reprod. Domest. Anim. 54, 786–794. doi:10.1111/rda.13426
Mukherjee, A., Koli, S., and Reddy, K. V. R. (2014). Regulatory Non-coding Transcripts in Spermatogenesis: Shedding Light on 'dark Matter'. Andrology 2, 360–369. doi:10.1111/j.2047-2927.2014.00183.x
Paiano, J., Wu, W., Yamada, S., Sciascia, N., Callen, E., Paola Cotrim, A., et al. (2020). ATM and PRDM9 Regulate SPO11-Bound Recombination Intermediates during Meiosis. Nat. Commun. 11, 857–915. doi:10.1038/s41467-020-14654-w
Parkhomchuk, D., Borodina, T., Amstislavskiy, V., Banaru, M., Hallen, L., Krobitsch, S., et al. (2009). Transcriptome Analysis by Strand-specific Sequencing of Complementary DNA. Nucleic Acids Res. 37, e123. doi:10.1093/nar/gkp596
Pertea, M., Kim, D., Pertea, G. M., Leek, J. T., and Salzberg, S. L. (2016). Transcript-level Expression Analysis of RNA-Seq Experiments with HISAT, StringTie and Ballgown. Nat. Protoc. 11, 1650–1667. doi:10.1038/nprot.2016.095
Quan, G., and Li, J. (2018). Circular RNAs: Biogenesis, Expression and Their Potential Roles in Reproduction. J. Ovarian Res. 11, 9–12. doi:10.1186/s13048-018-0381-4
Ramsköld, D., Wang, E. T., Burge, C. B., and Sandberg, R. (2009). An Abundance of Ubiquitously Expressed Genes Revealed by Tissue Transcriptome Sequence Data. Plos Comput. Biol. 5, e1000598. doi:10.1371/journal.pcbi.1000598
Rawlings, N., Evans, A., Chandolia, R., and Bagu, E. (2008). Sexual Maturation in the Bull. Reprod. Domest. Anim. 43, 295–301. doi:10.1111/j.1439-0531.2008.01177.x
Robles, V., Herráez, P., Labbé, C., Cabrita, E., Pšenička, M., Valcarce, D. G., et al. (2017). Molecular Basis of Spermatogenesis and Sperm Quality. Gen. Comp. Endocrinol. 245, 5–9. doi:10.1016/j.ygcen.2016.04.026
Sanger, H. L., Klotz, G., Riesner, D., Gross, H. J., and Kleinschmidt, A. K. (1976). Viroids Are Single-Stranded Covalently Closed Circular RNA Molecules Existing as Highly Base-Paired Rod-like Structures. Proc. Natl. Acad. Sci. 73, 3852–3856. doi:10.1073/pnas.73.11.3852
Shen, M., Wu, P., Li, T., Wu, P., Chen, F., Chen, L., et al. (2020). Transcriptome Analysis of circRNA and mRNA in Theca Cells during Follicular Development in Chickens. Genes 11, 489. doi:10.3390/genes11050489
Sherman, B. T., and Lempicki, R. A. (2009). Systematic and Integrative Analysis of Large Gene Lists Using DAVID Bioinformatics Resources. Nat. Protoc. 4, 44. doi:10.1038/nprot.2008.211
Song, H., Zhu, L., Li, Y., Ma, C., Guan, K., Xia, X., et al. (2015). Exploiting RNA-Sequencing Data from the Porcine Testes to Identify the Key Genes Involved in Spermatogenesis in Large White Pigs. Gene 573, 303–309. doi:10.1016/j.gene.2015.07.057
Soumillon, M., Necsulea, A., Weier, M., Brawand, D., Zhang, X., Gu, H., et al. (2013). Cellular Source and Mechanisms of High Transcriptome Complexity in the Mammalian Testis. Cel Rep. 3, 2179–2190. doi:10.1016/j.celrep.2013.05.031
Sun, L., Luo, H., Bu, D., Zhao, G., Yu, K., Zhang, C., et al. (2013). Utilizing Sequence Intrinsic Composition to Classify Protein-Coding and Long Non-coding Transcripts. Nucleic Acids Res. 41, e166. doi:10.1093/nar/gkt646
Svingen, T., and Koopman, P. (2013). Building the Mammalian Testis: Origins, Differentiation, and Assembly of the Component Cell Populations. Genes Develop. 27, 2409–2426. doi:10.1101/gad.228080.113
Tian, T. V., Di Stefano, B., Stik, G., Vila-Casadesús, M., Sardina, J. L., Vidal, E., et al. (2019). Whsc1 Links Pluripotency Exit with Mesendoderm Specification. Nat. Cel Biol 21, 824–834. doi:10.1038/s41556-019-0342-1
Tomlinson, M., Jennings, A., Macrae, A., and Truyers, I. (2017). The Value of Trans-scrotal Ultrasonography at Bull Breeding Soundness Evaluation (BBSE): The Relationship between Testicular Parenchymal Pixel Intensity and Semen Quality. Theriogenology 89, 169–177. doi:10.1016/j.theriogenology.2016.10.020
Trapnell, C., Williams, B. A., Pertea, G., Mortazavi, A., Kwan, G., Van Baren, M. J., et al. (2010). Transcript Assembly and Quantification by RNA-Seq Reveals Unannotated Transcripts and Isoform Switching during Cell Differentiation. Nat. Biotechnol. 28, 511–515. doi:10.1038/nbt.1621
Tscherner, A., Gilchrist, G., Smith, N., Blondin, P., Gillis, D., and Lamarre, J. (2014). MicroRNA-34 Family Expression in Bovine Gametes and Preimplantation Embryos. Reprod. Biol. Endocrinol. 12, 85–89. doi:10.1186/1477-7827-12-85
Ulitsky, I., and Bartel, D. P. (2013). lincRNAs: Genomics, Evolution, and Mechanisms. Cell 154, 26–46. doi:10.1016/j.cell.2013.06.020
Wang, K., Long, B., Liu, F., Wang, J.-X., Liu, C.-Y., Zhao, B., et al. (2016). A Circular RNA Protects the Heart from Pathological Hypertrophy and Heart Failure by Targeting miR-223. Eur. Heart J. 37, 2602–2611. doi:10.1093/eurheartj/ehv713
Waqas, M. S., Ciccarelli, M., Oatley, M. J., Kaucher, A. V., Tibary, A., and Oatley, J. M. (2019). Enhanced Sperm Production in Bulls Following Transient Induction of Hypothyroidism during Prepubertal Development. J. Anim. Sci. 97, 1468–1477. doi:10.1093/jas/sky480
Wu, S., Mipam, T., Xu, C., Zhao, W., Shah, M. A., Yi, C., et al. (2020). Testis Transcriptome Profiling Identified Genes Involved in Spermatogenic Arrest of Cattleyak. PloS one 15, e0229503. doi:10.1371/journal.pone.0229503
Yadav, R. P., and Kotaja, N. (2014). Small RNAs in Spermatogenesis. Mol. Cell. Endocrinol. 382, 498–508. doi:10.1016/j.mce.2013.04.015
You, X., Vlatkovic, I., Babic, A., Will, T., Epstein, I., Tushev, G., et al. (2015). Neural Circular RNAs Are Derived from Synaptic Genes and Regulated by Development and Plasticity. Nat. Neurosci. 18, 603–610. doi:10.1038/nn.3975
Young, J. C., Wakitani, S., and Loveland, K. L. (2015). TGF-β Superfamily Signaling in Testis Formation and Early Male Germline Development. Semin. Cel. Dev. Biol. 45, 94–103. doi:10.1016/j.semcdb.2015.10.029
Young, M. D., Wakefield, M. J., Smyth, G. K., and Oshlack, A. (2010). Gene Ontology Analysis for RNA-Seq: Accounting for Selection Bias. Genome Biol. 11, R14–R12. doi:10.1186/gb-2010-11-2-r14
Yu, Z., Guo, R., Ge, Y., Ma, J., Guan, J., Li, S., et al. (2003). Gene Expression Profiles in Different Stages of Mouse Spermatogenic Cells during Spermatogenesis1. Biol. Reprod. 69, 37–47. doi:10.1095/biolreprod.102.012609
Zhang, F., Zhang, X., Ning, W., Zhang, X., Ru, Z., Wang, S., et al. (2021). Expression Analysis of Circular RNAs in Young and Sexually Mature Boar Testes. Animals 11, 1430. doi:10.3390/ani11051430
Zhang, S., Zhang, Y., Yang, C., Zhang, W., Ju, Z., Wang, X., et al. (2015). TNP1 Functional SNPs in Bta-miR-532 and Bta-miR-204 Target Sites Are Associated with Semen Quality Traits in Chinese Holstein Bulls. Biol. Reprod. 92 (139), 139–110. doi:10.1095/biolreprod.114.126672
Zhou, G. H., Liu, L., Xiu, X. L., Jian, H. M., Wang, L. Z., Sun, B. Z., et al. (2001). Productivity and Carcass Characteristics of Pure and Crossbred Chinese Yellow Cattle. Meat Sci. 58, 359–362. doi:10.1016/s0309-1740(00)00160-1
Keywords: wandong cattle, testicular growth, spermatogenesis, circRNAs, total RNA sequencing
Citation: Khan IM, Liu H, Zhuang J, Khan NM, Zhang D, Chen J, Xu T, Avalos LFC, Zhou X and Zhang Y (2021) Circular RNA Expression and Regulation Profiling in Testicular Tissues of Immature and Mature Wandong Cattle (Bos taurus). Front. Genet. 12:685541. doi: 10.3389/fgene.2021.685541
Received: 25 March 2021; Accepted: 13 October 2021;
Published: 22 November 2021.
Edited by:
Mudasir Ahmad Syed, Sher-e-Kashmir University of Agricultural Sciences and Technology, IndiaReviewed by:
Olanrewaju B. Morenikeji, University of Pittsburgh at Bradford, United StatesRuihua Dang, Northwest A&F University, China
Copyright © 2021 Khan, Liu, Zhuang, Khan, Zhang, Chen, Xu, Avalos, Zhou and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Hongyu Liu, liuhongyu@ahau.edu.cn; Yunhai Zhang, yunhaizhang@ahau.edu.cn
†These authors have contributed equally to this work