Skip to main content

Whole-transcriptome sequencing analysis to identify key circRNAs, miRNAs, and mRNAs in the development of yak testes

Abstract

Background

The Testis is an important reproductive organ in male mammals and the site for spermatogenesis, androgen synthesis, and secretion. Non-coding RNAs (ncRNAs) play an important regulatory role in various biological processes. However, the regulatory role of ncRNAs in the development of yak testes and spermatogenesis remains largely unclear.

Result

In this study, we compared the expression profiles of circular RNAs (circRNAs), microRNAs (miRNAs), and messenger RNAs (mRNAs) in yak testicular tissue samples collected at 6 months (Y6M), 18 months (Y18M), and 4 years (Y4Y). Using RNA sequencing (RNA-Seq), we observed a significant difference in the expression patterns of ncRNAs in the samples collected at different testicular development stages. Twenty-two differentially expressed (DE) circRNAs, 69 DE miRNAs, and 64 DE mRNAs were detected in Y6M, Y18M, and Y4Y testicular samples, respectively. The results of gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses showed that the source genes of DE circRNAs, predicted target genes of DE miRNAs, and DE mRNAs were specifically associated with signaling pathways and GO terms that were related to sperm synthesis, sperm vitality, and testicular development, such as cell cycle, Wnt signaling pathway, MAPK signaling pathway, GnRH signaling pathway, and spermatogenesis. The analysis of the circRNA-miRNA-mRNA network revealed that some DE ncRNAs, including miR-574, miR-449a, CDC42, and CYP11A1, among others, may be involved in testicular spermatogenesis. Concurrently, various circRNA-miRNA interaction pairs were observed.

Conclusion

Our findings provide a database of circRNAs, miRNAs, and mRNAs expression profiles in testicular tissue of yaks at different developmental stages and a detailed understanding of the regulatory network of ncRNAs in yak testicular development and provide data that can help elucidate the molecular mechanisms underlying yak testicular development.

Peer Review reports

Introduction

The yak (Bos grunniens) is a unique breed of livestock in the Qinghai Tibet Plateau (QTP) and its adjacent areas. They have the compact body structures, are resistant to cold, and exhibit strong stress resistance. They serve as “versatile livestock” for people in high-altitude pastoral areas and produce a multitude of products, including meat, milk, plush, leather, and fuel, among others [1, 2].The testis is one of the most important reproductive organs in male animals and the foundation of male fertility. The level of male fertility directly influences the economic benefits from a specific species in animal husbandry [3]. Testicular development can be divided into prenatal and postnatal stages. The prenatal stage primarily involves the formation of testicular tissue. The postnatal stage involves the growth and differentiation of testicular tissue. When male animals are born, the testes are mainly composed of Sertoli cells (SCs) and Spermatogonia (SPG), and the diameter of the seminiferous tubules is relatively small. The continuous development of the testicles leads to the gradual enlargement of the diameter of the seminiferous tubules. Until adolescence, spermatogenesis gradually performe within the seminiferous tubules [4, 5]. Spermatogenesis is a highly regulated process in germ cell differentiation and comprises three stages: spermatogonium differentiation, spermatocyte meiosis, and sperm cell formation [6]. In recent years, several functional genes and ncRNAs have been found to regulate testicular development and spermatogenesis [7]. However, the regulatory mechanism remains unclear.

Circular RNAs (circRNAs) and microRNAs (miRNAs) are the two types of endogenous non-coding RNAs (ncRNAs) that have garnered significant attention in recent years. Most of these circRNAs have been shown to function as miRNA sponges that prevent miRNAs from binding to targets. miRNAs play a post-transcriptional regulatory role by specifically binding to protein-coding mRNAs [8, 9]. Several studies have shown that circRNAs and miRNAs are widely present in the testicular tissues of various animals. These play a crucial role in testicular development and spermatogenesis. For example, a total of 8,140 differentially expressed (DE) circRNAs were found in goat testicular tissues collected at four developmental stages [10]. Similarly, 223 DE miRNAs were identified in the testicular tissues of neonatal and mature cattle, with significantly enriched target genes in the calcium signaling, oxytocin signaling, PI3K Akt signaling, and MAPK signaling pathways, which are associated with male reproduction [11]. Each species has differences in the composition of testicular tissues, and the sperm production pathways may also differ.

In this study, we analyzed the expression of circRNAs, miRNAs, and mRNAs in yak testicular tissues collected at three developmental stages using whole-transcriptome RNA sequencing to reveal the molecular regulatory mechanisms of protein coding RNA (mRNA) and non-coding levels (miRNAs and circRNAs) on testicular development. At the same time, we conducted the functional enrichment analysis to screen for DE ncRNAs with a potential role in testicular development and spermatogenesis. Our results indicate that the potential function of ncRNAs in reproduction regulation in male yaks.

Materials and methods

Ethics statement

All experimental animals used in this study were treated strictly in accordance with the Animal Ethics Procedures and Guidelines of the People’s Republic of China. The experimental procedures were reviewed and approved by the Animal Administration and Ethics Committee of Lanzhou Institute of Husbandry and Pharmaceutical Sciences, Chinese Academy of Agricultural Sciences (approval number: SYXK-2014-0002).

Collection of testicular tissue samples

Nine healthy 6-month-old (Y6M), 18-month-old (Y18M), and 4-year-old (Y4Y) male yaks from the same cooperative in Xiahe County (N34°51’50’, E102°26’9’) in Gannan Tibetan Autonomous Prefecture were selected. Three male yaks were taken at each period and testicular tissue was collected via surgical biopsy. Briefly, the yaks were anesthetized with Jingsongling (Zhengfeng, Lanzhou, China) and wiped with iodophor. After local anesthetization, the yaks were castrated using sterile scalpels. The wounds of these yaks were immediately sutured with surgical needles, and penicillin and streptomycin were administered to prevent infection. The collected testicular tissue was washed with PBS (SolarBio, Beijing, China), immediately placed into a cryotube, and rapidly frozen in liquid nitrogen.

RNA isolation, library preparation, and sequencing

Total RNA was isolated from yak testicular tissues and purified using a TRIzol reagent kit (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. A 1% agarose gel was used to detect RNA degradation and contamination. The purity and concentration of total RNA were measured using a Nanodrop 2000 (Thermo Scientific, MA, USA) instrument and a Qubit® 2.0 fluorometer (Life Technologies, Carlsbad, CA, USA). The RNA integrity number (RIN) was determined using an Agilent 2100 Bioanalyzer (Agilent, Palo Alto, CA, USA). Only high-quality RNA samples with RNA concentration > 80 ng/µL, a purity ranging from 1.80 to 2.10, and RIN > 7 were used for the subsequent experiments. Ribosomal RNA (rRNA) was removed from high-quality RNA samples using a Ribo-Zero Gold rRNA Removal Kit (Illumina, San Diego, CA, USA). The residual RNA was used to generate 9 complementary DNA (cDNA) libraries using a NEBNext Ultra RNA Library Prep Kit (New England Biolabs, Ipswich, MA, USA). The libraries constructed were subjected to paired-end sequencing using the Illumina HiSeq2500 sequencer (Illumina, San Diego, CA, USA) from Gene Denovo Biotechnology Co., Ltd. (Guangzhou, China).

Processing and analysis of RNA-seq data from testicular tissues

High-quality clean reads were obtained using fastp v0.18.0 to eliminate unqualified reads, including adapters-containing reads, low-quality reads with quality score < Q20, and reads with more than 10% unknown nucleotides. First, the clean reads were aligned to the Bos Taurus reference genome (version: Ensembl_release106) using HISAT2 ver.2.1.0 with default parameters [12]. Second, the mapping reads for each sample were assembled using StringTie (version 1.3.4).

Library preparation and construction of small RNA-Seq

After total RNA extraction, RNA with 18–30 nt was concentrated by polyacrylamide gel electrophoresis. Next, the 3’ and 5’ adapters were connected separately, following which reverse transcription and PCR amplification of small RNAs were performed to connect nucleotides to both sides of the adapters. Finally, PCR products with a size of approximately 140 bp were recovered and purified for constructing cDNA libraries. The constructed libraries were assessed for quality and yield using an Agilent 2100 instrument and RT-quantitative PCR (RT-qPCR). Finally, small RNA-Seq was performed.

Identification and annotation of circRNAs

For reads that had not been aligned to the B. taurus reference genome, 20 bp at both ends were intercepted as anchor reads. Anchored reads were also realigned with the reference genome using bowtie2 v2.2.8. To avoid high false positives from circRNAs, the alignment results were identified using Find_circ [13]. To correct the expression levels of annotated circRNA, the Reads Per Million mapped reads (RPM) were used to standardize the raw counts.

Bioinformatics analysis of miRNA

Clean reads were obtained according to the method in 2.4, and then these were mapped to GenBank database v209.0 and Rfam database v11.0 for annotation, excluding other non-coding RNAs such as ribosomal RNA (rRNA), transfer RNA (tRNA), small nucleolar RNA (snoRNA), small nuclear RNA (snRNA), and small cytoplasmic RNA (scRNA). To eliminate exons, introns, and repeat sequences, these removed clean reads were mapped to the reference genome, and the remaining clean reads were annotated with miRbase v22.0 to obtain known bovine miRNAs and known miRNAs from other species (named known species-conserved miRNA). Finally, miReap v.0.2 was used to predict novel miRNAs using the signature hairpin structure of miRNA precursors. miRNA abundance was normalized using the transcripts per million (TPM) value according to established criteria [14].

Analysis of DE genes

The use of FPKM standardizes the mRNA expression level. CircRNAs, miRNAs, and mRNAs in the yak testicular tissue samples collected at three developmental stages were subjected to differential expression analysis using the DESeq v2.0 software. Only genes with |Fold Change| > 2.0 and p-value < 0.05 between any comparison groups were considered to show differential expression.

Enrichment analysis of GO and KEGG pathways for DE genes

The miReap v0.2 [15], Miranda v3.3a [16], and TargetScan v7.0 [17] were used to predict the target DE genes of miRNAs (DE miRNAs) and determine the intersection of the predicted results. The GO database was used to map DE genes, host genes of DE circRNAs, and target genes of DE miRNAs to GO terms. In addition, KEGG analysis was conducted to analyze whether the functional pathways of these DE ncRNAs are better for exploring their potential roles.

Construction of co-expression networks

A circRNA-miRNA-mRNA interaction network was developed using Cytoscape v3.5.1 [18] to investigate the interactions among circRNA, miRNA, and mRNA during spermatogenesis and testicular development in yaks at different developmental stages.

Assessing the reliability of RNA-Seq results using reverse transcriptase-PCR (RT-PCR)

Fifteen DE ncRNAs (five each for DE circRNAs, DE miRNAs, and DE mRNAs) were randomly selected for RT-qPCR to confirm the reliability of the sequencing results. First, primers were designed using Primer 5, followed by the reverse transcription of cDNA using the Transcriptor First Strand cDNA synthesis kit and Mir-X miRNA First-Strand Synthesis kit, respectively. The authenticity of circRNAs was assessed using RT-PCR and Sanger sequencing. The sequencing results were compared using MAG 5.0. RT-qPCR in an Applied Biosystems QuantStudio 6 Flex (Thermo Lifetech, MA, USA) platform 2 × ChamQ SYBR qPCR Master system (Vazyme, Nanjing, China), with three replicates used for each experiment. GAPDH expression was used as a reference to normalize the circRNA and mRNA expression levels, and U6 were used as the internal references to normalize the expression levels of miRNA. The threshold cycles were calculated using the 2− ∆∆ Ct method. The primer information is shown in Table 1.

Table 1 Primer information for RT-PCR and PT-qPCR

Results

Overview of RNA-Seq data

The cDNA libraries of testicular tissues collected at the three developmental stages of yaks were constructed and sequenced for each tissue sample. Transcriptome sequencing data were submitted to the GEO database with GSE268999. In all, 850,634,688 raw reads were collected in the three libraries. After low-quality and unqualified reads were eliminated, an average of 96,121,452, 94,463,369 and 92,257,787 clean reads were obtained from the testicular tissues of Y6M, Y18M, and Y4Y yaks, respectively. More than 93% of the clean reads were successfully mapped to the reference genome, of which 87.45% (Y6M), 88.82% (Y18M), and 89.79% (Y4Y) were uniquely aligned to the Bos Taurus reference genome. These findings indicated that the obtained data had high reliability and helped ensure the accuracy of subsequent analyses.

Overview of small RNA-Seq data

An average of 10,712,387, 9,575,322, and 14,836,913 raw reads were obtained from the testicular tissue samples of Y6M, Y18M, and Y4Y yaks, respectively. All raw reads generated in this investigation were submitted to GEO with GSE268999. After low-quality and unqualified reads were removed, an average of 10,419,535, 9,223,855, and 14,123,021 clean reads were obtained from the Y6M, Y18M, and Y4Y libraries, of which 80.50%, 73.62%, and 68.63% clean reads were mapped to the yak reference genome, respectively. Among the reads obtained from the small RNA library, most small RNAs had a length of 20–23 nt, and the most common length for small RNAs was 22 nt, accounting for at least 43.11%, 32.25%, and 26.65% of the total reads from Y6M, Y18M, and Y4Y samples, respectively (Fig. 1). We also found that among the small RNAs annotated in this study, the known miRNAs (including mature yak miRNAs and species-conserved miRNAs) constituted the largest proportion of the small RNA reads from Y6M, Y18M, and Y4Y samples, at 82.39%, 63.08%, and 51.74%, respectively (Fig. 2).

Fig. 1
figure 1

The nucleotide length distribution of small RNA obtained from testicular tissues of yaks at different developmental stages

Fig. 2
figure 2

The percentage of small RNA types obtained from testicular tissues of yaks at different developmental stages

Analysis of circRNA, miRNA, and mRNA characteristics

In this study, a total of 73,679 circRNAs, 1,536 miRNAs, and 22,183 mRNA were identified in testicular tissue samples collected in the three developmental stages. Among them, 31,386, 44,215, and 39,776 circRNAs were identified in the testicular tissues of Y6M, Y18M, and Y4Y yaks, respectively, and 12,825 circRNAs were co-expressed in three groups (Fig. 3A, Supplementary File S1). We found that among the different types of circRNAs, annot-exon was the most common sequence, accounting for 67% of the total number of circRNAs. This was followed by exon intron (19%) and intergenic (9%) circRNAs (Fig. 3B). CircRNAs with no more than 1000 bp accounted for the majority (Fig. 3C). In addition, based on their host gene location, the circRNAs identified in the testicular tissues were widely distributed in the yak chromosomes, with the greatest proportion observed in chromosomes 1–11 and least in MT chromosomes (Fig. 3D). The miRNAs identified in this study included 558 known yak miRNAs, 717 known species-conserved miRNAs, and 261 novel miRNAs (Supplementary File S2). Among the identified miRNAs, 1,210 were co-expressed in yak testicular tissues at the three developmental stages, whereas only 10, 1, and 34 miRNAs were expressed in the Y6M, Y18M, and Y4Y testicular tissues, respectively (Fig. 4).

Fig. 3
figure 3

The general characteristics of circRNA in testicular tissues of yaks at different developmental stages (A) The number of expressed circRNAs in testicular tissues of yaks at Y6M, Y18M, and Y4Y. (B) The type of circRNA. (C) The length distribution of circRNA. (D) Chromosomal distribution of circRNAs

Fig. 4
figure 4

The expression characteristics of miRNAs

Differential expression of circRNA, miRNA, and mRNA in yak testicular tissues collected at the three developmental stages

We compared the yak testicular tissue samples collected at the three developmental stages pairwise to evaluate the different circRNA expression levels. Of the 73,679 circRNAs identified, 1,841 DE circRNAs were identified in Y6M and Y18M samples, of which 672 were up-regulated and 1,169 were down-regulated. Next, 3,962 DE circRNAs were identified in Y4Y and Y6M samples, of which 2,255 were up-regulated and 1,707 were downregulated. Lastly, 908 DE circRNAs were identified in Y4Y and Y18M samples, of which 398 were up-regulated and 510 were down-regulated (Fig. 5A).

Among all miRNAs identified in the yak testicular tissues collected at the three different developmental stages, 462 miRNAs showed significant differential expression between the Y6M and Y18M stages (154 up-regulated miRNAs and 308 downregulated miRNAs). A total of 676 DE miRNAs were identified in Y4Y and Y6M samples (440 up-regulated miRNAs and 236 downregulated miRNAs). Similarly, 176 DE miRNAs were identified in Y4Y and Y18M samples (102 up-regulated miRNAs and 74 downregulated miRNAs) (Fig. 5B).

Subsequently, DE mRNAs were screened using DESeq with P < 0.05 and |log2 (fold change) | ≥ 1 as the screening criteria. Through screening, 7,668 DE mRNAs (2,315 up-regulated mRNAs and 5,353 down-regulated mRNAs) and 10,565 DE mRNAs (including 6,391 up-regulated mRNAs and 4,174 down-regulated mRNAs) were identified in the Y6M vs. Y18M and Y4Y vs. Y6M comparisons, respectively. Meanwhile, 247 DE mRNAs were identified in the Y4Y vs. Y18M comparison (64 up-regulated and 183 down-regulated), which was relatively lesser than those identified in other pairs (Fig. 5C).

Fig. 5
figure 5

The differential expression of DE ncRNAs in yak testicular tissues of Y6M, Y18M, and Y4Y

Prediction of target genes for DE miRNAs

To further investigate the potential functions and mechanisms of action of DE miRNAs, miRanda and Targetscan were used to predict the target genes of the DE miRNAs. We identified 462, 676, and 176 DE miRNAs from Y6M vs. Y18M, Y4Y vs. Y6M, and Y4Y vs. Y18M comparisons and predicted 20,132, 20,375, and 18,924 target genes, respectively. The vast majority of miRNAs obtained had more than one target gene and could function by cleaving multiple targets. A single target gene can simultaneously target multiple miRNAs. For instance, the newly discovered novel-m0007-3p was predicted to have 1,153 targets, whereas ITGA6 was predicted to target 124 miRNAs simultaneously.

Functional enrichment analysis of DE ncRNAs in Y6M vs. Y18M comparison

We conducted GO and KEGG enrichment analyses of DE ncRNAs (including DE circRNAs, miRNAs, and genes) to further explore the alterations in yak testicular tissues at different developmental stages. Results of the GO enrichment analysis showed that the source genes of DE circRNAs, predicted target genes of DE miRNAs, and DE genes between Y6M and Y18M samples were significantly concentrated in 1,111, 3,849, and 1,828 GO terms, respectively (Supplementary File S3), primarily involving intracellular, metal ion binding, single-organism organelle organization, cellular component organization, and cell migration terms (Fig. 6A and C). Interestingly, some terms were closely associated with testicular development and reproduction traits, including spermatogenesis, sexual reproduction, transformation growth factor beta binding, and spermatid development.

The metabolic processes related to DE ncRNAs were analyzed by using the KEGG pathway database. The source genes of DE circRNAs, predicted target genes of DE miRNAs, and DE genes were significantly enriched in 64, 128, and 75 KEGG pathways, respectively. These were primarily involved in cell cycle, growth hormone synthesis, secretion and action, Wnt signaling pathway, GnRH signaling pathway, and MAPK signaling pathway (Fig. 6D and F).

Fig. 6
figure 6

GO and KEGG enrichment analysis of DE ncRNAs (including DE circRNAs, DE miRNAs, and DEGs) in Y6M vs. Y18M. (A) GO analysis of the source genes of DE circRNAs in the Y6M vs. Y18M group. (B) GO analysis of the predicted target genes of DE miRNAs in the Y6M vs. Y18M group. (C) GO analysis of DEGs in the Y6M vs. Y18M group. (D) KEGG analysis of the source genes of DE circRNAs in the Y6M vs. Y18M group. (E) KEGG analysis of the predicted target genes of DE miRNAs in the Y6M vs. Y18M group. (F) KEGG analysis of DEGs in the Y6M vs. Y18M group

Functional enrichment analysis of DE ncRNAs in Y4Y vs. Y6M comparison

Between Y4Y and Y6M samples, the source genes of DE circRNAs, predicted target genes of DE miRNAs, and DE genes were significantly enriched in 1,590, 3,837, and 2,631 GO terms, respectively (Supplementary File S4). Among these, the most important cellular component functions involved the intracellular, intracellular membrane-bounded organelle, organelle, and cytoplasm. The molecular functions comprised binding, protein binding, catalytic activity, MAP kinase activity, and transforming growth factor beta receptor binding. With respect to biological processes, the most important GO terms primarily included spermatogenesis, phosphatidylinositol 3-kinase signaling, sperm mobility, DNA methylation, and Wnt signaling pathway (Fig. 7A and C). KEGG analysis of DE ncRNAs (including DE circRNAs, miRNAs, and genes) indicated that the parental genes were involved in the FoxO signaling pathway, mTOR signaling pathway, PI3K Akt signaling pathway, MAPK signaling pathway, and GnRH signaling pathway (Fig. 7D and F).

Fig. 7
figure 7

GO and KEGG enrichment analysis of DE ncRNAs (including DE circRNAs, DE miRNAs, and DEGs) in Y4Y vs. Y6M. (A) GO analysis of the source genes of DE circRNAs in the Y4Y vs. Y6M group. (B) GO analysis of the predicted target genes of DE miRNAs in the Y4Y vs. Y6M group. (C) GO analysis of DEGs in the Y4Y vs. Y6M group. (D) KEGG analysis of the source genes of DE circRNAs in the Y4Y vs. Y6M group. (E) KEGG analysis of the predicted target genes of DE miRNAs in the Y4Y vs. Y6M group. (F) KEGG analysis of DEGs in the Y4Y vs. Y6M group

Functional enrichment analysis of DE ncRNA in Y4Y vs. Y18M

Similarly, between Y4Y and Y18M samples, DE ncRNAs were significantly enriched in 5,815 GO terms, with the source genes of DE circRNAs, the predicted target genes of DE miRNAs, and DEGs significantly enriched in 1,200, 3,973, and 766 GO terms (Supplementary File S5), respectively, among which spermatid development, cell migration, sex differentiation, transforming growth factor beta receptor signaling pathway, and catalytic activity may play a role in testicular development (Fig. 8A and C). The results of KEGG pathway annotation revealed that the source genes of DE circRNAs in the Y4Y and Y18M samples were significantly enriched in 56 pathways, among which GnRH signaling pathway, cell cycle, Wnt signaling pathway, MAPK signaling pathway, and PI3K-Akt signaling pathway may contribute to testicular development (Fig. 8D and F).

Fig. 8
figure 8

GO and KEGG enrichment analysis of DE ncRNAs (including DE circRNAs, DE miRNAs, and DEGs) in Y4Y vs. Y18M. (A) GO analysis of the source genes of DE circRNAs in the Y4Y vs. Y18M group. (B) GO analysis of the predicted target genes of DE miRNAs in the Y4Y vs. Y6M group. (C) GO analysis of DEGs in the Y4Y vs. Y18M group. (D) KEGG analysis of the source genes of DE circRNAs in the Y4Y vs. Y18M group. (E) KEGG analysis of the predicted target genes of DE miRNAs in the Y4Y vs. Y18M group. (F) KEGG analysis of DEGs in the Y4Y vs. Y18M group

Construction and investigation of the competitive endogenous RNA (ceRNA) regulatory network

To further determine the potential regulatory effects of DE ncRNAs on testicular development and male reproduction in yaks at different developmental stages, we randomly selected 12 DE circRNAs of type annot exons and their targeted DE miRNAs and mRNAs to construct a circRNA -miRNA-mRNA interaction network based on the ceRNA theory, the results of which are shown in Fig. 9.

Fig. 9
figure 9

The circRNA -miRNA-mRNA interaction network constructed with DE circRNAs and its target miRNA and mRNA in leaves. The gray circles represent DE circRNA, the green inverted triangles represent the predicted target miRNA of circRNA, and the red triangles represent the target gene of miRNA

Validating the authenticity and reliability of RNA-Seq data using RT-qPCR

The results of RT-PCR and Sanger sequencing showed that the randomly selected circRNAs in the RNA-seq data had the head-to-tail splice junction sequences, further demonstrating the authenticity of circRNAs in the RNA-seq data (Fig. 10A). RT-qPCR is used to detect the expression level of DE ncRNAs in yak testicular tissue. Although there is a certain quantitative difference between RT-qPCR and RNA-seq, the relative expression levels (log2 FC) between the two methods are similar (Fig. 10B). Therefore, the RNA-seq data are reliable and can be used for subsequent analysis.

Fig. 10
figure 10

Verification of the authenticity and reliability of RNA Seq data. (A) RT-PCR and Sanger sequencing were used to verify the authenticity of circRNAs. The red arrow indicates the junction site on the DNA sequence. (B) RT-qPCR was used to verify the expression level of DE ncRNA in yak testicular tissue

Discussion

In male animals, the production, synthesis, and secretion of sperm occur in the testicles. Normal testicular development is crucial for sperm production and fertility maintenance in male mammals [19]. In this case, it is crucial to study the molecular mechanisms underlying testicular development and spermatogenesis for improving livestock fertility. We investigated the expression profiles of circRNA, miRNA, and mRNA in the testicular tissues of yaks during infancy (Y6M), adolescence (Y18M), and adulthood (Y4Y), and constructed a ceRNA regulatory network indicative of yak testicular growth and development and spermatogenesis. We also explored the impact of expression changes on testicular development and spermatogenesis.

Previous studies have shown that testicular development and spermatogenesis were not only regulated by hormones but also rely on the stringent regulation of the spatiotemporal expression of protein-coding RNA and ncRNA [20]. It is known that mRNA serves as a bridge between DNA and proteins, responsible for transmitting genetic information stored in DNA and guiding protein synthesis in cells [21]. Up to now, an increasing number of mRNAs have been shown to be involved in testicular development and spermatogenesis processes, and ncRNAs also play an important role in protein translation [22, 23]. MiRNAs induce gene degradation or suppresses gene expression by targeting mRNAs. CircRNAs regulate gene expression by competitively binding to miRNAs [24, 25]. In this study, only 64 DE mRNAs, 69 DE miRNAs, and 22 DE circRNAs were identified to have differential expression in all three stages, further confirming the temporal expression pattern of ncRNAs. In addition, we compared Y4Y vs. Y18 samples with Y6M vs. Y18M samples and Y4Y vs. Y6M samples and found that there were a greater number of DE circRNAs, DEmiRNAs, and DE mRNAs in the Y6M vs. Y18M and Y4Y vs. Y6M groups, with significant differences in expression levels. It can be seen that the testicular development and spermatogenesis of yaks at different stages may be achieved by regulating the quantity and expression level of various RNAs. This is consistent with the research results of La et al. (2023).

Here, we conducted functional enrichment analysis to explore the potential regulatory effects of ncRNA on testicular development and spermatogenesis. As expected, DE ncRNAs in the yak testicular tissues samples collected at 6 months, 18 months, and 4 years were involved in various biological processes, including cell cycle, protein binding, spermatogenesis, metal ion binding, and several metabolic processes. An expanding body of evidence supports the fact that metal ion binding is involved in testicular development and spermatogenesis [26, 27]. Besides, for GO enrichment analysis of the three control groups, the parental genes of DE circRNAs in Y4Y vs. Y18M samples were enriched in cell morphogenesis and cell shape regulation. Similarly, many target mRNAs of DE miRNAs in the Y4Y vs. Y6M group were also enriched in cell shape regulation. This may be attributed to the complex deformation of sperm cells during spermatogenesis after the second meiotic division post puberty [28, 29].

Moreover, the KEGG analysis results of the three comparison groups showed that among the significantly enriched KEGG pathways, the GnRH signaling pathway, Wnt signaling pathway, MAPK signaling pathway, PI3K-Akt signaling pathway, p53 signaling pathway, and cell cycle were the most frequently identified. Most of these pathways were associated with spermatogenesis and testicular development. GnRH was reported to regulate testicular growth and development by regulating hormone synthesis and secretion [30]. The Wnt signaling pathway is involved in the regulation of various biological processes, including spermatogenesis. Tanwar et al. (2010) [31] showed that the Wnt signaling pathway influences spermatogenesis by regulating the differentiation of SCs. The MAPK signaling is involved in the regulation of SC self-renewal and meiosis and is closely associated with spermatogenesis [32]. PI3K-Akt regulates spermatogonia proliferation and spermatogenesis by regulating mTOR signaling and downstream target protein p70S6K/4EBP1 expression [33, 34]. Interestingly, DE genes in Y4Y vs. Y6M and Y6M vs. Y18M comparisons were enriched in the Hippo signaling pathway. Reportedly, the Hippo signaling pathway is involved in regulating steroid production and plays a crucial role in adolescent development and male sexual maturity [35]. Therefore, the above functional enrichment results indicate that the identified DE ncRNAs are involved in yak testes development and spermatogenesis, and the degrees of development in the yak testes at the three stages differ significantly.

The ceRNA hypothesis reveals a novel mechanism of RNA-to-RNA interaction, which focuses on miRNAs and connects mRNAs with lncRNAs/circRNAs through microRNA response elements (MREs) to jointly regulate gene transcription and expression in the body [36, 37]. CeRNAs can act as a bait to competitively bind to miRNA, further relieving the inhibitory effect of miRNAs on target genes [38]. The involvement of the ceRNA regulatory network in the regulation of testicular growth, development, and spermatogenesis in cattle [19], goats [10], camels [39], and pigs [40] has been reported in recent years. In this study, we constructed a ceRNA regulatory network related to testicular development and spermatogenesis to better elucidate the molecular mechanisms underlying testicular development in yaks. In the ceRNA regulatory network, most miRNAs were shown to function by specifically binding to targets. For example, miRNA-449 and miRNA-34b/c regulate the development of male germ cells in mice by targeting genes in the E2F Transcription Factor-Retinoblastoma Protein (E2F-pRb) pathway [41]. MiR-184 regulates mouse spermatogenesis by downregulating nuclear receptor corepressor 2 (Ncor2) expression [6]. KLF4 downregulation mediated by miR-92a-3p may play a key role in limiting the maturation of mouse SCs [42]. Here, we found that miR-92a, miR-92b, and miR-574 target KLF4. As widely known, KLF4 is a key transcription factor that regulates apoptosis and can induce the pluripotency of spermatogonial stem cells, thus playing an important role in germ cells [43]. CYP11A1, PDGFRB, LIMK2, and AXL in the ceRNA network are also reportedly involved in testicular development, spermatogenesis, and androgen synthesis [44,45,46,47]. Therefore, in the ceRNA network, the interactions between novel_circ_002426-bta-miR-92a-KLF4, novel_circ_016117-bta-miR-34c-CYP11A1, and novel_circ_000745-bta-miR-449a-CYP11A1 may be involved in testicular development and spermatogenesis. Additionally, we observed a significant overlap between the potential mRNA and circRNA targets of some miRNAs. CircRNAs are widely believed to act as a “sponge” for multiple miRNAs simultaneously, and an miRNA can also target multiple mRNAs simultaneously [48].

To summarize, we constructed a ceRNA regulatory network involved in testicular development and spermatogenesis. During the regulatory process, many circRNAs act as ceRNAs and compete with mRNAs to bind to miRNAs, indirectly altering the expression level of mRNAs related to testicular growth and development and eventually regulating testicular development and spermatogenesis. In future research, we will focus more on the circRNA-miRNA-mRNA interactions potentially involved in testicular development and spermatogenesis, further exploring their biological functions and providing comprehensive recommendations for investigations on the regulatory mechanisms underlying testicular development and spermatogenesis.

Conclusion

In this study, we comprehensively evaluated the expression patterns of circRNAs, miRNAs, and mRNAs during the growth and development of yaks using testicular tissues. ncRNAs, such as circRNAs and miRNAs, were found to be actively expressed in yak testicles. The functions of DE circRNAs, DE miRNAs, and DE mRNAs in the testicular tissues collected during the three developmental stages of yaks were associated with testicular development and spermatogenesis. In addition, we constructed a ceRNA network related to testicular development and spermatogenesis and screened the interactions between novel_circ_002426-bta-miR-92a-KLF4, novel_circ_016117-bta-miR-34c-CYP11A1, and novel_circ_000745-bta-miR-449a-CYP11A1, all of which potentially regulate testicular development and spermatogenesis. However, the specific roles of these candidate ncRNAs require further validation. Our results provide novel insights into the regulatory mechanisms underlying yak testicular development and spermatogenesis.

Data availability

Sequence data that support the findings of this study have been deposited in the GEO repository with the primary accession code GSE268999.

References

  1. Wang H, Chai Z, Hu D, Ji Q, Xin J, Zhang C, Zhong J. A global analysis of CNVs in diverse yak populations using whole-genome resequencing. BMC Genomics. 2019;20:1–12.

    Google Scholar 

  2. Krishnan G, Paul V, Biswas T, Chouhan V, Das P, Sejian V. Adaptation strategies of yak to seasonally driven environmental temperatures in its natural habitat. Int J Biometeorol. 2018;62:1497–506.

    Article  CAS  PubMed  Google Scholar 

  3. Mäkelä J-A, Koskenniemi JJ, Virtanen HE, Toppari J. Testis development. Endocr Rev. 2019;40(4):857–905.

    Article  PubMed  Google Scholar 

  4. Mأ¤kelأ¤ J-A, Koskenniemi JJ, Virtanen HE, Toppari J. Testis development. Endocr Rev. 2019;40(4):857–905.

    Article  Google Scholar 

  5. de Kretser DM, Loveland KL, Meinhardt A, Simorangkir D, Wreford N. Spermatogenesis. Hum Reprod. 1998;13(suppl1):1–8.

    Article  PubMed  Google Scholar 

  6. Wu J, Bao J, Wang L, Hu Y, Xu C. MicroRNA-184 downregulates nuclear receptor corepressor 2 in mouse spermatogenesis. BMC Dev Biol. 2011;11:1–10.

    Article  Google Scholar 

  7. Das PP, Begum SS, Choudhury M, Medhi D, Paul V, Das PJ. Characterizing miRNA and mse-tsRNA in fertile and subfertile yak bull spermatozoa from Arunachal Pradesh. J Genet. 2020;99:1–9.

    Article  Google Scholar 

  8. Hansen TB, Jensen TI, Clausen BH, Bramsen JB, Finsen B, Damgaard CK, Kjems J. Natural RNA circles function as efficient microRNA sponges. Nature. 2013;495(7441):384–8.

    Article  CAS  PubMed  Google Scholar 

  9. Bartel DP. MicroRNAs: target recognition and regulatory functions. cell 2009, 136(2):215–233.

  10. Tang W, Xu QH, Chen X, Guo W, Ao Z, Fu K, Ji T, Zou Y, Chen JJ, Zhang Y. Transcriptome sequencing reveals the effects of circRNA on testicular development and spermatogenesis in Qianbei Ma Goats. Front Veterinary Sci. 2023;10:1167758.

    Article  Google Scholar 

  11. Gao Y, Wu F, Ren Y, Zhou Z, Chen N, Huang Y, Lei C, Chen H, Dang R. MiRNAs expression profiling of bovine (Bos taurus) testes and effect of Bta-miR-146b on proliferation and apoptosis in bovine male germline stem cells. Int J Mol Sci. 2020;21(11):3846.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Memczak S, Jens M, Elefsinioti A, Torti F, Krueger J, Rybak A, Maier L, Mackowiak SD, Gregersen LH, Munschauer M. Circular RNAs are a large class of animal RNAs with regulatory potency. Nature. 2013;495(7441):333–8.

    Article  CAS  PubMed  Google Scholar 

  14. Zhou L, Chen J, Li Z, Li X, Hu X, Huang Y, Zhao X, Liang C, Wang Y, Sun L. Integrated profiling of microRNAs and mRNAs: microRNAs located on Xq27. 3 associate with clear cell renal cell carcinoma. PLoS ONE. 2010;5(12):e15224.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Liu W, Xu L, Wang Y, Shen H, Zhu X, Zhang K, Chen Y, Yu R, Limera C, Liu L. Transcriptome-wide analysis of chromium-stress responsive microRNAs to explore miRNA-mediated regulatory networks in radish (Raphanus sativus L). Sci Rep. 2015;5(1):14024.

    Article  PubMed  PubMed Central  Google Scholar 

  16. Turner DA. Miranda: A non-strict functional language with polymorphic types. In: Conference on Functional Programming Languages and Computer Architecture: 1985: Springer; 1985: 1–16.

  17. Lewis BP, Burge CB, Bartel DP. Conserved seed pairing, often flanked by adenosines, indicates that thousands of human genes are microRNA targets. cell 2005, 120(1):15–20.

  18. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. La Y, Ma X, Bao P, Chu M, Yan P, Guo X, Liang C. Identification and characterization of piwi-interacting RNAs for early testicular development in yak. Int J Mol Sci. 2022;23(20):12320.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Kimmins S, Sassone-Corsi P. Chromatin remodelling and epigenetic features of germ cells. Nature. 2005;434(7033):583–9.

    Article  CAS  PubMed  Google Scholar 

  21. Pardi N, Hogan MJ, Porter FW, Weissman D. mRNA vaccinesقa new era in vaccinology. Nat Rev Drug Discovery. 2018;17(4):261–79.

    Article  CAS  PubMed  Google Scholar 

  22. Zou Y, Chen X, Tian X, Guo W, Ruan Y, Tang W, Fu K, Ji T. Transcriptomic analysis of the developing testis and spermatogenesis in Qianbei Ma Goats. Genes. 2023;14(7):1334.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Hombach S, Kretz M. Non-coding RNAs: classification, biology and functioning. Non-coding RNAs Colorectal cancer 2016:3–17.

  24. Mahesh G, Biswas R. MicroRNA-155: a master regulator of inflammation. J Interferon Cytokine Res. 2019;39(6):321–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Wilusz JE. A 360 view of circular RNAs: from biogenesis to functions. Wiley Interdisciplinary Reviews: RNA. 2018;9(4):e1478.

    Article  PubMed  Google Scholar 

  26. Luo Z, Liu Y, Chen L, Ellis M, Li M, Wang J, Zhang Y, Fu P, Wang K, Li X. microRNA profiling in three main stages during porcine spermatogenesis. J Assist Reprod Genet. 2015;32:451–60.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Han H, Chen Q, Gao Y, Li J, Li W, Dang R, Lei C. Comparative transcriptomics analysis of testicular miRNA from cryptorchid and normal horses. Animals. 2020;10(2):338.

    Article  PubMed  PubMed Central  Google Scholar 

  28. Hess RA, De Franca LR. Spermatogenesis and cycle of the seminiferous epithelium. Mol Mech Spermatogenesis 2009:1–15.

  29. Yao M, Qu H, Han Y, Cheng CY, Xiao X. Kinesins in mammalian spermatogenesis and germ cell transport. Front Cell Dev Biology. 2022;10:837542.

    Article  Google Scholar 

  30. Desaulniers AT, Cederberg RA, Mills GA, Ford JJ, Lents CA, White BR. LH-independent testosterone secretion is mediated by the interaction between GnRH2 and its receptor within porcine testes. Biol Reprod. 2015;93(2):45.

    Article  PubMed  Google Scholar 

  31. Tanwar PS, Kaneko-Tarui T, Zhang L, Rani P, Taketo MM, Teixeira J. Constitutive WNT/beta-catenin signaling in murine sertoli cells disrupts their differentiation and ability to support spermatogenesis. Biol Reprod. 2010;82(2):422–32.

    Article  CAS  PubMed  Google Scholar 

  32. Ni F-D, Hao S-L, Yang W-X. Multiple signaling pathways in sertoli cells: recent findings in spermatogenesis. Cell Death Dis. 2019;10(8):541.

    Article  PubMed  PubMed Central  Google Scholar 

  33. Dibble CC, Cantley LC. Regulation of mTORC1 by PI3K signaling. Trends Cell Biol. 2015;25(9):545–55.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Xu H, Shen L, Chen X, Ding Y, He J, Zhu J, Wang Y, Liu X. mTOR/P70S6K promotes spermatogonia proliferation and spermatogenesis in Sprague Dawley rats. Reprod Biomed Online. 2016;32(2):207–17.

    Article  CAS  PubMed  Google Scholar 

  35. Zhang X, Cheng Z, Wang L, Jiao B, Yang H, Wang X. MiR-21-3p centric regulatory network in dairy cow mammary epithelial cell proliferation. J Agric Food Chem. 2019;67(40):11137–47.

    Article  CAS  PubMed  Google Scholar 

  36. Poliseno L, Salmena L, Zhang J, Carver B, Haveman WJ, Pandolfi PP. A coding-independent function of gene and pseudogene mRNAs regulates tumour biology. Nature. 2010;465(7301):1033–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi PP. A ceRNA hypothesis: the Rosetta Stone of a hidden RNA language? Cell. 2011;146(3):353–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Quan J, Kang Y, Luo Z, Zhao G, Li L, Liu Z. Integrated analysis of the responses of a circRNA-miRNA-mRNA ceRNA network to heat stress in rainbow trout (Oncorhynchus mykiss) liver. BMC Genomics. 2021;22:1–10.

    Article  Google Scholar 

  39. Hasi G, Sodnompil T, Na H, Liu H, Ji M, Xie W, Nasenochir N. Whole transcriptome sequencing reveals core genes related to spermatogenesis in bactrian camels. J Anim Sci. 2023;101:skad115.

    Article  PubMed  PubMed Central  Google Scholar 

  40. Wang P, Liu Z, Zhang X, Huo H, Wang L, Dai H, Yang F, Zhao G, Huo J. Integrated analysis of lncRNA, miRNA and mRNA expression profiles reveals regulatory pathways associated with pig testis function. Genomics 2024:110819.

  41. Bao J, Li D, Wang L, Wu J, Hu Y, Wang Z, Chen Y, Cao X, Jiang C, Yan W. MicroRNA-449 and microRNA-34b/c function redundantly in murine testes by targeting E2F transcription factor-retinoblastoma protein (E2F-pRb) pathway. J Biol Chem. 2012;287(26):21686–98.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Gupta A, Vats A, Ghosal A, Mandal K, Sarkar R, Bhattacharya I, Das S, Pal R, Majumdar SS. Follicle-stimulating hormone-mediated decline in miR-92a-3p expression in pubertal mice sertoli cells is crucial for germ cell differentiation and fertility. Cell Mol Life Sci. 2022;79(3):136.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Corbineau Sb, Lassalle B, Givelet M, Souissi-Sarahoui Is, Firlej V, Romeo PH, Allemand I, Riou L, Fouchet P. Spermatogonial stem cells and progenitors are refractory to reprogramming to pluripotency by the transcription factors Oct3/4, c-Myc, Sox2 and Klf4. Oncotarget 2017, 8(6):10050.

  44. Wang H, Wang X, Li T, An X, Chen N, Shi H, Su M, Ma K, Hao Z, Duan X. Differential tissue expression of sex steroid-synthesizing enzyme CYP11A1 in male tibetan sheep (Ovis aries). Animal Biotechnol. 2023;34(7):2900–9.

    Article  CAS  Google Scholar 

  45. Odeh HM, Kleinguetl C, Ge R, Zirkin BR, Chen H. Regulation of the proliferation and differentiation of Leydig stem cells in the adult testis. Biol Reprod. 2014;90(6):123.

    Article  PubMed  PubMed Central  Google Scholar 

  46. Takahashi H, Funakoshi H, Nakamura T. LIM-kinase as a regulator of actin dynamics in spermatogenesis. Cytogenetic Genome Res 2003, 103.

  47. Sasaki T, Knyazev PG, Clout NJ, Cheburkin Y, Gأhring W, Ullrich A, Timpl R, Hohenester E. Structural basis for Gas6قôAxl signalling. EMBO J. 2006;25(1):80–7.

    Article  CAS  PubMed  Google Scholar 

  48. Peter M. Targeting of mRNAs by multiple miRNAs: the next step. Oncogene. 2010;29(15):2161–4.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

This study was financially supported by the Xiahe County East and West Collaboration Technology Project (202301), China Agriculture Research System of MOF and MARA (CARS-37), and Innovation Project of Chinese Academy of Agricultural Sciences (25-LZIHPS-01).

Author information

Authors and Affiliations

Authors

Contributions

L.H. conceptualized this study and helped in data curation, formal analysis, software, validation, visualization and writing - original draft. X.W. helped in writing – review. S.G. helped in writing - review. M.C. helped in validation. Y.K. helped in sample collection. Z.D. helped in sample collection. J.P. helped in Supervision. Q.G. helped in Methodology. Y.M. helped in software. X.G. helped in funding acquisition, project administration and writing - review & editing.

Corresponding authors

Correspondence to Yi Ma or Xian Guo.

Ethics declarations

Ethics approval and consent to participate

All experimental animals used in this study were treated strictly in accordance with the Animal Ethics Procedures and Guidelines of the People’s Republic of China. The experimental procedures were reviewed and approved by the Animal Administration and Ethics Committee of Lanzhou Institute of Husbandry and Pharmaceutical Sciences, Chinese Academy of Agricultural Sciences (approval number: SYXK-2014-0002).

Consent for publication

Not applicable.

Competing interest

The authors declare no conflict of interest.

Additional information

Publisher’s note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Electronic supplementary material

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Hu, L., Wang, X., Guo, S. et al. Whole-transcriptome sequencing analysis to identify key circRNAs, miRNAs, and mRNAs in the development of yak testes. BMC Genomics 25, 824 (2024). https://doi.org/10.1186/s12864-024-10716-1

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-024-10716-1

Keywords