
==== Front
BMC Genomics
BMC Genomics
BMC Genomics
1471-2164
BioMed Central London

10716
10.1186/s12864-024-10716-1
Research
Whole-transcriptome sequencing analysis to identify key circRNAs, miRNAs, and mRNAs in the development of yak testes
Hu Liyan 12
Wang Xingdong 12
Guo Shaoke 12
Cao Mengli 12
Kang Yandong 12
Ding Ziqiang 12
Pei Jie 12
Ge Qianyun 12
Ma Yi tjnnyz@126.com

3
Guo Xian guoxian@caas.cn

12
1 grid.410727.7 0000 0001 0526 1937 Key Laboratory of Yak Breeding of Gansu Province, Lanzhou Institute of Husbandry and Pharmaceutical Sciences, Chinese Academy of Agricultural Sciences, Lanzhou, 730050 China
2 https://ror.org/05ckt8b96 grid.418524.e 0000 0004 0369 6250 Key Laboratory of Animal Genetics and Breeding on Tibetan Plateau, Ministry of Agriculture and Rural Affairs, Lanzhou, 730050 China
3 https://ror.org/0516wpz95 grid.464465.1 0000 0001 0103 2256 Tianjin Academy of Agriculture Sciences, Tianjin, 300192 China
2 9 2024
2 9 2024
2024
25 82414 6 2024
16 8 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ 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/.
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.

Supplementary Information

The online version contains supplementary material available at 10.1186/s12864-024-10716-1.

Keywords

Yak
Testis
circRNAs
miRNAs
mRNAs
the Xiahe County East and West Collaboration Technology Project202301 China Agriculture Research System of MOF and MARACARS-37 Innovation Project of Chinese Academy of Agricultural Sciences25-LZIHPS-01 issue-copyright-statement© BioMed Central Ltd., part of Springer Nature 2024
==== Body
pmcIntroduction

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

CircRNA/miRNA/gene	Forward (5’→3’)	Reverse (5’→3’)	
circ_002438	TCCAGACAAAGTGCGAGGTT	TTCCTGCATCTGCTGTCCTC	
circ_032419	AAGACAAAGCAGCACAAG	AAGCACATCCACCTCATC	
circ_046272	GCCCAGGACCCGGTTATT	TGCTGTGAAGCCAGGATG	
circ_046890	GTCCCATTCCCATAACTG	CTTGAGGCATATCCCACT	
circ_061335	CTTATCAGCCGCCTCTTT	CTTCCTGGTTGGTGTTCC	
circ_063280	CGCCAAGGGCAAGACCAA	ACCCGGAAGAGCAGCAAG	
circ_067598	CGGAGGTCAGTTCAATAC	TAAGAATGCTCGGTAGGT	
circ_069511	TTGGTTTATACATGGTCCTG	AGCTGTCGAAGCCACTTT	
circ_070696	CCTCAGGTAGACGGTGTT	GGGAGTTCCCTTGGTTAT	
bta-miR-146b	CCTGAGAACTGAATTCCATAGGCTGT	mRQ 3′ primer	
bta-miR-449a	TGGCAGTGTATTGTTAGCTGGT	mRQ 3′ primer	
bta-miR- 100	AACCCGTAGATCCGAACTTGTG	mRQ 3′ primer	
bta-miR-34c	CGAGGCAGTGTAGTTAGCTGATTG	mRQ 3′ primer	
bta-miR-574	TGAGTGTGTGTGTGTGAGTGTGTG	mRQ 3′ primer	
bta-miR-92b	ATATATTGCACTCGTCCCGGCCT	mRQ 3′ primer	
bta-miR-204	CGTTCCCTTTGTCATCCTATGCCT	mRQ 3′ primer	
bta-miR-92a	ATTATTGCACTTGTCCCGGCCT	mRQ 3′ primer	
bta-miR-184	CTGGACGGAGAACTGATAAGGGT	mRQ 3′ primer	
bta-miR-449a	TGGCAGTGTATTGTTAGCTGGT	mRQ 3′ primer	
bta-miR- 100	AACCCGTAGATCCGAACTTGTG	mRQ 3′ primer	
bta-miR-34c	CGAGGCAGTGTAGTTAGCTGATTG	mRQ 3′ primer	
CYP11A1	TCATCCCACTGCTGAATC	TCAAAGGCAAAGTGAAACA	
PDGFRB	TCATCTCCCTCATTATCCTCA	CATCTTCACAGCCACTTTCAT	
CTSB	CGCCAACAACGAGAAGGA	TTGCCGACCAGCCAGTAG	
KLF4	TCCCACCGCTCCATTACC	AAACTTCCACCCACAACCATC	
AURKA	GCCTTGTCTTACTGTCATTC	GCTCCAGAGGTCCACTTT	
U6	GGAACGATACAGAGAAGATTAGC	TGGAACGCTTCACGAATTTGCG	
GAPDH	ACACTGAGGACCAGGTTGTG	GACAAAGTGGTCGTTGAGGG	

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 The nucleotide length distribution of small RNA obtained from testicular tissues of yaks at different developmental stages

Fig. 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 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 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 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 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 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 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 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 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–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.

Electronic supplementary material

Below is the link to the electronic supplementary material.

Supplementary Material 1

Supplementary Material 2

Supplementary Material 3

Supplementary Material 4

Supplementary Material 5

Acknowledgements

Not applicable.

Author 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.

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).

Data availability

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

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.

Publisher’s note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
==== Refs
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 30606130
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.30606130
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 10.1007/s00484-018-1549-8 29728761
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.29728761 10.1007/s00484-018-1549-8
3. Mäkelä J-A Koskenniemi JJ Virtanen HE Toppari J Testis development Endocr Rev 2019 40 4 857 905 10.1210/er.2018-00140 30590466
Mäkelä J-A, Koskenniemi JJ, Virtanen HE, Toppari J. Testis development. Endocr Rev. 2019;40(4):857–905.30590466 10.1210/er.2018-00140
4. Mأ¤kelأ¤ J-A Koskenniemi JJ Virtanen HE Toppari J Testis development Endocr Rev 2019 40 4 857 905 10.1210/er.2018-00140 30590466
Mأ¤kelأ¤ J-A, Koskenniemi JJ, Virtanen HE, Toppari J. Testis development. Endocr Rev. 2019;40(4):857–905.30590466 10.1210/er.2018-00140
5. de Kretser DM Loveland KL Meinhardt A Simorangkir D Wreford N Spermatogenesis Hum Reprod 1998 13 suppl1 1 8 10.1093/humrep/13.suppl_1.1 9663765
de Kretser DM, Loveland KL, Meinhardt A, Simorangkir D, Wreford N. Spermatogenesis. Hum Reprod. 1998;13(suppl1):1–8.9663765 10.1093/humrep/13.suppl_1.1
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 10.1186/1471-213X-11-64 21194500
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.21194500 10.1186/1471-213X-11-64
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 10.1007/s12041-020-01248-0 32089520
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.32089520 10.1007/s12041-020-01248-0
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 10.1038/nature11993 23446346
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.23446346 10.1038/nature11993
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 10.3389/fvets.2023.1167758
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.10.3389/fvets.2023.1167758
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 10.3390/ijms21113846 32481702
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.32481702 10.3390/ijms21113846
12. Kim D Langmead B Salzberg SL HISAT: a fast spliced aligner with low memory requirements Nat Methods 2015 12 4 357 60 10.1038/nmeth.3317 25751142
Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60.25751142 10.1038/nmeth.3317
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 10.1038/nature11928 23446348
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.23446348 10.1038/nature11928
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 10.1371/journal.pone.0015224 21253009
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.21253009 10.1371/journal.pone.0015224
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 10.1038/srep14024 26357995
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.26357995 10.1038/srep14024
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 10.1101/gr.1239303 14597658
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.14597658 10.1101/gr.1239303
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 10.3390/ijms232012320 36293174
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.36293174 10.3390/ijms232012320
20. Kimmins S Sassone-Corsi P Chromatin remodelling and epigenetic features of germ cells Nature 2005 434 7033 583 9 10.1038/nature03368 15800613
Kimmins S, Sassone-Corsi P. Chromatin remodelling and epigenetic features of germ cells. Nature. 2005;434(7033):583–9.15800613 10.1038/nature03368
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 10.1038/nrd.2017.243
Pardi N, Hogan MJ, Porter FW, Weissman D. mRNA vaccinesقa new era in vaccinology. Nat Rev Drug Discovery. 2018;17(4):261–79.10.1038/nrd.2017.243
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 10.3390/genes14071334 37510239
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.37510239 10.3390/genes14071334
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 10.1089/jir.2018.0155 30998423
Mahesh G, Biswas R. MicroRNA-155: a master regulator of inflammation. J Interferon Cytokine Res. 2019;39(6):321–30.30998423 10.1089/jir.2018.0155
25. Wilusz JE A 360 view of circular RNAs: from biogenesis to functions Wiley Interdisciplinary Reviews: RNA 2018 9 4 e1478 10.1002/wrna.1478 29655315
Wilusz JE. A 360 view of circular RNAs: from biogenesis to functions. Wiley Interdisciplinary Reviews: RNA. 2018;9(4):e1478.29655315 10.1002/wrna.1478
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 10.1007/s10815-014-0406-x 25563581
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.25563581 10.1007/s10815-014-0406-x
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 10.3390/ani10020338 32098036
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.32098036 10.3390/ani10020338
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 10.3389/fcell.2022.837542
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.10.3389/fcell.2022.837542
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 10.1095/biolreprod.115.128082 26134865
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.26134865 10.1095/biolreprod.115.128082
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 10.1095/biolreprod.109.079335 19794154
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.19794154 10.1095/biolreprod.109.079335
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 10.1038/s41419-019-1782-z 31316051
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.31316051 10.1038/s41419-019-1782-z
33. Dibble CC Cantley LC Regulation of mTORC1 by PI3K signaling Trends Cell Biol 2015 25 9 545 55 10.1016/j.tcb.2015.06.002 26159692
Dibble CC, Cantley LC. Regulation of mTORC1 by PI3K signaling. Trends Cell Biol. 2015;25(9):545–55.26159692 10.1016/j.tcb.2015.06.002
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 10.1016/j.rbmo.2015.11.007 26706460
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.26706460 10.1016/j.rbmo.2015.11.007
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 10.1021/acs.jafc.9b04059 31532202
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.31532202 10.1021/acs.jafc.9b04059
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 10.1038/nature09144 20577206
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.20577206 10.1038/nature09144
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 10.1016/j.cell.2011.07.014 21802130
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.21802130 10.1016/j.cell.2011.07.014
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 10.1186/s12864-020-07335-x 33388042
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.33388042 10.1186/s12864-020-07335-x
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 10.1093/jas/skad115 37083698
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.37083698 10.1093/jas/skad115
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 10.1074/jbc.M111.328054 22570483
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.22570483 10.1074/jbc.M111.328054
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 10.1007/s00018-022-04174-9 35181820
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.35181820 10.1007/s00018-022-04174-9
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 10.1080/10495398.2022.2125401
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.10.1080/10495398.2022.2125401
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 10.1095/biolreprod.114.117473 24740597
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.24740597 10.1095/biolreprod.114.117473
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 10.1038/sj.emboj.7600912
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.10.1038/sj.emboj.7600912
48. Peter M Targeting of mRNAs by multiple miRNAs: the next step Oncogene 2010 29 15 2161 4 10.1038/onc.2010.59 20190803
Peter M. Targeting of mRNAs by multiple miRNAs: the next step. Oncogene. 2010;29(15):2161–4.20190803 10.1038/onc.2010.59
