Unveiling the Transcriptomic Landscape of Holstein Bull Testes Post-Sexual Maturation: Insights into Spermatogenesis-Associated MicroRNAs and mRNA Profiles

1. Introduction

Male reproductive efficiency is an essential component of cattle’s production, biological breeding or reproduction, and economic efficiency (Ansari-Lari et al., 2010; Han & Peñagaricano, 2016; Liu et al., 2022). The elite genetic gain is made through enhanced production of sperm and selective gametes utilization from desirable sires to improve production and reproduction characteristics in dairy and beef industry on large scales. Bull fertility is an important factor in dairy farm businesses with strong influence on breeding and conception of cows (Han & Peñagaricano, 2016; Waqas et al., 2019). Selection of good breeding soundness bulls may have increase 20-40% of fertility in getting cows pregnant. Improved bulls fertility increases conception rate and calf weaning weight, reduce culling females, and reduce calving season. Artificial insemination and single-sire mating groups can also increase fertility of bulls in dairy farms and ranches. Therefore, standard breeding soundness evaluation of bulls and their testis must be considered to solve the substantial deficit of fertility with limitation to detect sub fertile bulls in a consistent manner. Variations in fertility can be suggested with combinations of laboratory tests at cellular, molecular and animal levels (Kastelic & Thundathil, 2008).The testis is the site of the spermatogenesis process performing two major functions namely, exocrine function and endocrine function. The endocrine function of testis is production of differentiated gametocytes facilitated by Sertoli cells, which then released in to the lumen of seminiferous tubules and stored in the epididymis, while the endocrine function is supported by leydig cells to produce steroids (Ge et al., 2021; Staub & Johnson, 2018).

Spermatogenesis is a long, orderly sophisticated developmental process of gametes multiplication and differentiation into spermatozoa in the seminiferous tubules. Mammalian spermatogenesis is an androgen dependent process in which male germ line cells divide and differentiate (Grande et al., 2022; Sahlu et al., 2020; Smith & Walker, 2014; Zhao, Wang, et al., 2022). The undifferentiated spermatogonial germ cells starts to split and divided in to spermatocytes, spermatids and finally in to the production of mature spermatozoa, released in seminiferous tubules, at the apical pole of Sertoli cells and stored in the epididymis (Sahlu et al., 2020; Staub & Johnson, 2018). This process involves mitotic division, differentiation meiotic divisions and post meiotic divisions. Spermatogonial stem cells to produce type A spermatogonia developing in to type B spermatogonia that differentiate in to diploid primary spermatocytes (Zhao et al., 2021). The primary spermatocytes then divided into secondary spermatocytes during the meiosis I. The secondary spermatocytes also divided to produce four haploid spermatids. During spermiogenesis process, the haploid spermatids undergo morphological changes, including Chromatid reorganization, elimination of cytoplasm, flagellum development, nuclear formation and acrosomal formation to produce mature spermatozoa (Sahlu et al., 2020; Zhao et al., 2021). Testosterone (T), a well-known and prominent sexual steroid hormone that functions prepubertal development for the expression of specific transcripts, in males that plays a crucial role for Sertoli cells, function in pubertal period prior to testicular maturation, largely produced by the Leydig cells in the testes (Johnston et al., 2004). At different stages during the development of the testes, Leydig, Sertoli, and Germ cells are also capable of producing estrogen. Estrogen has an important function in early testicular development and is primarily produced by Sertoli cells, which are found in the immature testis regulating Leydig cells development and number (Baker et al., 2003; Carreau et al., 2012). As their testes continue to generate estrogen and control spermatogenesis throughout this time, prepubertal animals, have a significant quantity of estradiol (E2) in their serum (Schulster et al., 2016).

MicroRNAs(miRNAs) are endogenous transcriptomic elements of small noncoding RNAs that are identified in animals, plants and viruses with 20-25 nucleotide length controlling most of the genes expression by blocking translation or inducing degradation of target messenger RNAs at posttranscriptional and transcriptional levels (Baohong et al., 2007; Liu et al., 2022; Sun et al., 2019; Vigneault et al., 2012). The miRNAs help to mediate coding and noncoding RNA molecules during fine-tuning crosstalk by sharing miRNA response elements (MRE) (Ala, 2020). These are responsible for RNA silencing, gene expression through targeting mRNA for post-transcriptional endonucleolytic cleavage (Wu et al., 2014), and translational repression or mRNA decay (Ala, 2020; Baohong et al., 2007).They also play regulatory and mediation functions in the cross-talks between the coding and noncoding transcriptome components by recognizing the target sites in various biological functions and physiological processes (Ala, 2020; Sun et al., 2019).The regulation of miRNAs is modulated by noncoding RNAs that are significantly expressed in many biological and cellular processes.

The high-tech sequencing technologies such as second generation sequencing technologies, micro RNA deep sequencing and Illumina-Solexa high throughput sequencing technologies speed up the identification and discovery of new miRNA species (Li et al., 2012; Wu et al., 2014).Multiple studies shown that miRNAs play an important role in the regulation of various biological processes and tissue systems that function as regulation of cellular differentiation, self-renewal and apoptosis (Niu et al., 2011), including mammalian spermatogenesis and reproductive system development (La et al., 2023; Liu et al., 2022; Wang et al., 2018; Wu et al., 2014; Yin et al., 2021; Zhang et al., 2022). Lots of studies indicate the involvement of miRNAs modulating gene expression in multiple cellular and biological processes, including proliferation (Kim et al., 2017), organ development, stem cell differentiation (Luo et al., 2020), hormone secretion, metabolism and reproductive regulation (Gao et al., 2020; La et al., 2023; Liu et al., 2022; Wu et al., 2014). Additionally, miRNAs can serve as epigenetic modulators, influencing target mRNAs’ protein levels without altering gene sequences. The expression of miRNAs is influenced by the regulation of epigenetic machineries, such as DNA methylation and RNA modification (Yao et al., 2019). Su et al. reported the expression profiles of miRNAs in different tissues (bovine sera, exosomes, mammary gland, rumen and gut) (Sun et al., 2019). This study identified 639 miRNAs providing a clue to the understanding of miRNAs and their regulatory roles in tissue specific and tissue development functions that can be applied in animal production, animal health and human health research. Puberty is closely related with sexual and testicular maturity, spermatogenesis and capability of reproduction. There is increased evidence of miRNAs involvement in regulating testicular development and spermatogenesis. Zhang et al. screened out 322 novel miRNAs, 767 known miRNAs and 359 exist miRNAs regulating precocious sexual maturity traits providing new light on miRNAs role and regulation mechanisms in Land race and Hezuo boars for precocious puberty (Zhang et al., 2022). A study by Zhu et al. on Saanen Dairy Goat Testis found conserved miRNAs through Illumina high-throughput sequencing technology and miRNA sequences blasted to miRNAs of cow and sheep in miRBase 19.0 identifying 373 conserved miRNAs and 91 novel paired-miRNAs while those blasted to miRNAs of cow and mouse resulted in 128 conserved miRNA, which may reveal their involvement in meiosis and testis development in the dairy goat testis (Wu et al., 2014). An investigation by Liu et al. conducted a comprehensive sequencing analysis study using testicular tissue samples from three immature calves and three sexually mature bulls of the Wandong cattle breed identified 151 differentially expressed (DE) miRNAs and miRNA-linked genes. The miRNA-linked genes found involved in PI3KAkt signaling, TNF signaling and T cell receptor of biological pathways and functions affecting spermatogenesis and testicular development. The DE expressed miRNAs were also associated with sexual maturation and sperm production in Wandong Cattle (Bos taurus) (Liu et al., 2022).

The identification of single stranded noncoding RNA is important to determine the functional mechanisms during sexual maturation and testicular growth at transcriptional and post transcriptional levels (Liu et al., 2022).The improvement of reproductive efficiency in breeding Holstein bulls with emphasis on potential genetic resources regulating spermatogenesis and testis development at molecular level remains unexplored. Even though, the comprehensive expression profiles of miRNA in different age groups of cattle breeds and different species have been studied, the spermatogenesis related miRNAs expression in Holstein bulls has not been reported. Therefore, there is a gap on the expression of miRNAs, regulatory mechanisms and expression profiles in the developing Holstein bull testis and remains uncertain. The final goal of the study was to determine potential biomarkers related to Holstein bulls testes and find miRNA associated candidate mRNAs or target genes so as to improve Holstein bull reproductive performance for the future.

I study the insights of miRNA and its related mRNA accompanying spermatogenesis and testicular growth to further broaden the understanding of miRNA regulation of spermatogenesis and testicular development in immature and mature Holstein bull testes was briefly described.

2. Materials and Methods

2.1. Ethical Approval

All animal protocols, experimental designs and procedures required for study were approved by the Animal Care and Use Committee of the Institute of Animal Sciences, Chinese Academy of Agricultural Sciences, Beijing, China. All the required standards with necessary measures were made to minimize and avoid suffering of the experimental animal from pain.

2.2. Experimental Animals and Testis Tissue Sample Collection

Before the start of the experiment, the experimental animals were checked for their health status and six animals selected from those healthy and free of disease. Animals were based on similar genetic background of the same breed and body condition according to their female parent’s better calving interval and categorized in to two according to their age. Three 80 weeks old mature Holstein bull (n=3; designated as D5, D6 and D7) and three 8 weeks young Holstein bulls (n=3; designated as X1, X2 and X3) were filtered and surgically castrated to obtain intact testis after being given general anesthesia (Zoletil 50, Virbac Co., France). The experimental animals were properly restrained in a clean environment before and after sedation. Then castration started following the excision of the testis with removing the lower third of the scrotum and testis cutting the spermatic cord by surgical blade in a manner of reduced or with no bleeding. The connective tissue, tunica albuginea and fat were removed using forceps under the stereomicroscope cutting the middle segments of the testis parenchymal tissue, quarter of each pair testes near the epididymis tail, incised into1-2cm³ size tissues, collected with RNAlater and snap frozen in (-196 ^oC) liquid nitrogen immediately for isolation. The frozen sample then transported to Institute of Animal Science, Key Laboratory of Animal Genetics, Breeding and Reproduction and stored at -80 ^oC until RNA extraction.

2.3. Determination of Hormonal Concentration

A total of six Holstein bulls were used, comprising three young and three mature Holstein bulls sourced from the same dairy farm in Beijing, China. These bulls were divided into two groups, with three bulls assigned to each group, namely 8 week old and 80 week old. Blood samples were collected from all six bulls through the tail vein venipuncture. Within one hour of collection, the samples were centrifuged at 3,000 rpm for 5 minutes to collect serum and then stored at -20°C for further analysis for T (testosterone), E2 (estradiol), and P4 (progesterone) hormone concentrations.

To determine the concentrations of T, E2, and P4 in each bull’s serum, a radioimmunoassay system (XH6080, China) was used. The measurements were performed in triplicate as per the manufacturer’s instructions. Subsequently, the serum levels of T, E2, and P4 were analyzed using one-way analysis of variance (ANOVA) in the SPSS software.

2.4. Histological Assessment

Following the removal of the fascial tissue, tunica vaginalis, and tunica albuginea, a testicular sample was obtained by cross-cutting a slice from the middle region of the testes using precise dissection techniques. The slice was preserved in a 4% formaldehyde solution to enable subsequent histological evaluation. The testicular tissues were further processed, sectioned into sections measuring 6μm thickness, and stained with Masson’s trichrome stain. The stained sections were examined under different magnification powers to assess the histomorphology of the testicular tissues. A comprehensive analysis of the structural characteristics and composition of the testis in the two groups was achieved, allowing for a detailed understanding of its intricate histological features in between the group.

2.5. RNA Extraction and Qualification

Initially, the frozen testicular tissues were pulverized into a fine powder. Subsequently, the extraction of tissue RNA was carried out using a trizol reagent kit (Thermo Fisher, 15596-026), and the quantification of total RNA was accurately performed using the qubit RNA assay br kit (Invitrogen, Q10211). The ribo-zero magnetic kit (Lucigen, MRZPL116) was employed to eliminate ribosomal RNA (rRNA) from the total RNA. The agilent RNA 6000 pico chip (Agilent, 5067-1511) was utilized for quality control of the rRNA that had been removed from the total RNA.

2.6. Quality Control of Sequence Reads

The raw data was filtered using third party software with the filtering criteria of discarding reads with 5’-end adapter contaminants, trim 3’-end adapter; reads with 3’-end adapter only, reads with length <18nt, reads with N base, reads with poly-A in 3’-end, reads with length >30nt and low quality reads (more than 70% bases with quality score <20). The FastQC was used for the quality checks on the raw data of the six testicular tissue samples resulting in high quality sequences. The adapter sequences (Adapter, Oligonucleotide sequences for TrueSeq^TM small RNA sample prep kits of RNA 5’ Adapter (RA5), part# 15013205 GUUCAGAGUUCUACAGUCCGACG AUC and RNA 3’ Adapter (RA3), part # 15013207 AGATCGGAAGAGCACACGTCT were used to filter out the sequence reads. The Cutadapt-1.7.1 and fastx-toolkit-0.0.14 were used to remove the adaptors or the connector sequences and low quality reads, respectively. The clean reads statistical results were used for comparisons and sequence analysis.

2.7. Reads Mapping and Sequence Analysis

The cDNA libraries were constructed from testicular tissues of three young Holstein bulls (X1, X2, X3; group YB) and three mature Holstein bulls (D5, D6, D7; group AB). The sequence reads for each sample was from the Illumina high-throughput sequencing. They compared to reference genome with HISAT22.1.0 software (Institute Pasteur, Paris, France) (Kim et al., 2015) according to the criteria. Cutadapt-1.7.1 was used to remove the adapters and fastx-toolkit-0.0.14 to eliminate the low quality reads. The reads were aligned to the Rfam database and mapped using the Blastin. The number and percentage of clean reads of pre-miRNA, rRNA, tRNA, snRNA, other ncRNAs and repeats were compared to sample species miRNA precursor sequence database (miRBase v21), Rfam database (release11.0, rRNAsequences) and the Silva database (rRNAdatabase), Rfam database (release 11.0, tRNA sequences), Rfam database (release 11.0, snRNA sequences), Rfam database (release11.0, except miRNA, tRNA, snRNA, rRNA sequences) and repeat region with the RepBase (version 20.04) database, respectively, as a base for classification statistics of high quality clean data. The filtered sequences that do not matched to those databases are not classified. Bovine rRNA, sRNA, snRNA, and snoRNA mapping reads were eliminated. The sRNAs were categorized, and readings were mapped to the bovine reference genome (ARS-UCD1.3), both in exon and intron sequences, using Bowtie v1.1.1; Boao Crytal Code Program. The retained reads that were mapped to the bovine genome and intron region were further processed for novel miRNAs prediction and quantification of known miRNAs via the miRDeep2 v2.0.0.6 module and quantifier module of the miRDeep2 software (Friedländer et al., 2012).

2.8. Predicting miRNAs

In order to predict miRNAs, the first step involved was extracting candidate precursor sequences from the reference genome by extending the sequences of Clean Reads at 20nt upstream and 70nt downstream, as well as 70nt upstream and 20nt downstream, with a total length of less than 120nt. The extracted sequences were comprehensively scored based on the thermodynamic stability of the hairpin structure and secondary structure of the candidate pre-miRNA sequences. The miRDeep2 software predicted precursor pre-miRNAs that satisfy specific conditions, including more than 60% of the bases in the mature miRNA sequence being paired with star miRNAs, a length difference between mature miRNA and star miRNA of less than 6nt, and no branching in the precursor structure. The predicted pre-miRNAs were filtered based on criteria such as the total number of Clean Reads compared to pre-miRNA being greater than or equal to 2 and the total number of Clean Reads compared to mature miRNAs and star miRNAs as a percentage of the total number of precursor Clean Reads being greater than or equal to 20%. Additionally, the starting position of Clean Reads was within 2nt overhang upstream of the 5’-end of the known miRNA sequence and the ending position of clean reads were within 2nt overhang downstream of the 3’end of the known miRNA sequence. MiRDeep2 analysis software was used and resulted in novel miRNA prediction, secondary structure, length and base distribution, and comparison display map.

2.9. Expression Normalization and Correlation Analysis of Samples

After localizing clean reads specifically to miRNAs, the numbers of localized reads in the miRNAs were utilized to determine the expression level of miRNAs. Since the number of reads is also dependent on the length and sequencing depth, normalization using TPM (Transcripts Reads Number per Million) was done before expression analysis. The denominator in the equation represents the total number of reads in a sample as compared to the reference genome, while the numerator represents the total number of reads in a specific miRNA. The normalized TPM values (TPM=1,000,000*(The miRNA Mapped Reads/Total Mapped Reads) represent the expression of miRNA in a given sample, and these values were applied in subsequent analysis. The Boao Jingdian program is the expression normalization analysis software utilized in the process to determine the expression of miRNA in the sample groups more accurately and creatively. To conduct analysis of the sample correlation utilizing the R language; the Pearson Correlation Coefficient was calculated. This value, which falls between the range of -1 and 1, was determined using the TPM values of both the samples of known miRNA and novel miRNA. The Pearson Correlation Coefficient serves as an essential indicator of biological replicate correlation, with a value closer to 1 suggest a stronger correlation between the two replicates.

2.10. Differential Expression Analysis

The expression level of miRNAs was compared by the R package edgeR/limma package for the two grouped testicular tissue samples. Differential miRNA expression analysis was performed on grouped samples to find up- and down-regulated miRNAs. Specifically, the lowly expressed miRNAs were filtered with the TPM normalization method. To screen the differentially expressed miRNAs between the two group samples, the significant P-value and fold-change (FC) were used as the criteria for differential expression: (1) Judged Up-regulated miRNAs when P-value < 0.05 and FC> 2 (log2 (FC) > 1); (2) Down-regulated miRNAs was judged when p-value < 0.05 and FC < 0.5 (log2(FC) < -1). Larger difference multiples of FC values indicated that miRNA expression was more different in the two grouped samples. The differentially expressed miRNAs in known miRNAs and novel miRNAs were systematically clustered separately to construct differential expression profiles, and miRNAs with similar expression patterns in multiple samples would be clustered together. Volcano plot was carried out to merge common DE miRNAs between the two groups for further analysis. The differentially expressed miRNAs in known miRNAs and novel miRNAs were systematically clustered separately to construct differential expression profiles, and miRNAs with similar expression patterns in multiple samples would be clustered together.

2.11. Target Prediction and Function Annotation

With the miRNAs of the two age groups (X vs D) being filtered, miRanda was employed to predict their target genes. The threshold of the prediction software miRanda was set to Score > 140 and Energy < -20 kcal/mol. Finally, the target gene prediction results were taken from the top 10 of the miRanda software prediction results. The miRanda Prediction principle was used for miRanda Target Gene Prediction Software by Memorial Sloan-Kettering Cancer Research Center (https://www.mskcc.org/), published in 2003. miRanda predicts target genes based on the complementarity of miRNA sequences and candidate 3’UTR or mRNA sequences, as well as the thermal stability of miRNA-mRNA duplexes. The miRNAs interact with mRNAs and can regulate mRNA expression at the translational level through complementary pairing. Each miRNA can have multiple target genes, and several miRNAs can regulate the same gene. This complex regulatory network can either regulate the expression of multiple genes through a single miRNA or finely regulate the expression of a gene through a combination of several miRNAs. Thus, the 3’UTR region of mRNAs as the target gene prediction database was used for selection of the miRNAs. Therefore, the gene set enrichment analysis (GSEA) was performed to gain a deeper understanding on the DEGs expression patterns of miRNAs and their function extrapolating the data for profound processes and pathways.

2.12. Network Analysis of Differentially Expressed miRNA Target Genes

In order to comprehensively understand the complex relationships between functional processes and genes with their respective pathways the powerful tool of igraph R package (the network analysis package) was utilized to illustrate net plots and provide a clear depiction of the correlations among pathways. In particular, the study focused on creating diagrammatic representations of the GO to GO, gene to pathways, and pathway to pathway plots.

2.13. RT-qPCR Validation and Statistical Analysis

To validate the sequencing results, an experiment was conducted on miRNAs and target genes chosen randomly for qRT-PCR validation and U6 was used as a house keeping gene. Further experimentation was done on the selected up/down regulated miRNAs and mRNAs using specific primers and conditions. The 2^-∆∆CT method was utilized for quantifying the relative expression of miRNAs and target genes. A student’s t-test was conducted using Graph Pad Prism 8 software for comparing the differences among the young and adult Holstein bull testis tissue samples. The p-value < 0.05 was considered statistically significant while p-value < 0.01 was considered extremely significant.

Table 1. The information on the primers for mRNA genes used in RT-qPCR validation.

Gene name	Sequence of primers (5’- 3’)	Product length (bp)	Tm (◦C)
THAP3	F: TCTGGGGAGAGAAAGGTCCTC R: AGGGCATAGCTGTGATCGGA	227	60.27
FOSB	F: GAGAAGAGAAGGGTTCGCCG R: CTAGCTGATCTGTCTCCGCC	106	60.46
SLC2A5	F: CTAGCTGATCTGTCTCCGCC R: GTCGACGGTGGAAACTCCTT	203	60.03
CCNT2	F: ACCCAGTTAGTAAGAGCAAGCA R: CAGGGATCTCCCAATTGGACC	160	59.36
CMC1	F: ACCCCTCAGAGCAGCAT R: TGGGAAGCTTCTGTAGCCTTT	300	57.33
GAPDH	F: CTTCGGCATTGTGGAGGG R: GGAGGCAGGGATGATGTTCT	130	61.3

Table 1. The information on the primers for mRNA genes used in RT-qPCR validation.

Gene name	Sequence of primers (5’- 3’)	Product length (bp)	Tm (◦C)
THAP3	F: TCTGGGGAGAGAAAGGTCCTC R: AGGGCATAGCTGTGATCGGA	227	60.27
FOSB	F: GAGAAGAGAAGGGTTCGCCG R: CTAGCTGATCTGTCTCCGCC	106	60.46
SLC2A5	F: CTAGCTGATCTGTCTCCGCC R: GTCGACGGTGGAAACTCCTT	203	60.03
CCNT2	F: ACCCAGTTAGTAAGAGCAAGCA R: CAGGGATCTCCCAATTGGACC	160	59.36
CMC1	F: ACCCCTCAGAGCAGCAT R: TGGGAAGCTTCTGTAGCCTTT	300	57.33
GAPDH	F: CTTCGGCATTGTGGAGGG R: GGAGGCAGGGATGATGTTCT	130	61.3

The U6 for RT primer: 5’--CTCAACTGGTGTCGTGGAGTCGGCAATTCAGTTGAGAAAAATATG--3’; Forward primer: 5’--CTC GCTTCGGCAGCACA--3’; Reverse qPCR Primer: 5’--CTCGCTTCACG AATTTGCGTGTCAT--3’. Thus, U6 was used as housekeeping gene for the validation of novel and known miRNAs.

Table 2. Stem-loop RT-PCR primers for real-time quantification of selected differentially expressed miRNAs in testis between Young bull and Adult bull.

MiRNA ID (5′-3′)	Stem-loop RT primer (5′-3′)	MiRNA specific forward primer
bta-miR-11971	5’--GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACTTCTCA --3’	5’--CGCGTGAGGGGCAGAGAG --3’
bta-miR-2284c	5’-- GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACAAAACC --3’	5’--GCGCGAAAAAGTTCGTTTT --3’
bta-miR-3956	5’--GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACCCTCTG --3’	5’--CGACGTGGATGCTGAAGGT --3’
bta-miR-135b	5’--GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACTCACAT --3’	5’--CGCGTATGGCTTTTCATTCCT --3’
bta-miR-299-2	5’--GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACATGGTT --3’	5’--CGCGGTATGTGGGACGGTA --3’
bta-miR-433	5’--GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACACACCG --3’	5’--GCGATCATGATGGGCTCCT --3’
bta-miR-191	5’--GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACCAGCTG --3’	5’--CGCAACGGAATCCCAAAAG --3’
bta-miR-20a	5’--GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACCTACCT --3’	5’--GCGCGTAAAGTGCTTATAGTGC --3’
4_37054-5p(cli-miR-9632-3p)	5’--GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACAGTCCG --3’	5’--GCGAACTTTTGCCCCTAGTAA --3’
7_42582-5p	5’--GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACTGTTAG --3’	5’--GCGCGTGAGAGGTCTGTAATAA --3’
6_40791-3p	5’--GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACGGTCCC --3’	5’--GCGTGAGTGTGCATCCACG --3’
21_23590-5p(mmu-miR-665-3p)	5’--GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACGCTCAT --3’	5’-- CGCCCAGGAGAATGAGGG --3’

Table 2. Stem-loop RT-PCR primers for real-time quantification of selected differentially expressed miRNAs in testis between Young bull and Adult bull.

MiRNA ID (5′-3′)	Stem-loop RT primer (5′-3′)	MiRNA specific forward primer
bta-miR-11971	5’--GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACTTCTCA --3’	5’--CGCGTGAGGGGCAGAGAG --3’
bta-miR-2284c	5’-- GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACAAAACC --3’	5’--GCGCGAAAAAGTTCGTTTT --3’
bta-miR-3956	5’--GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACCCTCTG --3’	5’--CGACGTGGATGCTGAAGGT --3’
bta-miR-135b	5’--GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACTCACAT --3’	5’--CGCGTATGGCTTTTCATTCCT --3’
bta-miR-299-2	5’--GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACATGGTT --3’	5’--CGCGGTATGTGGGACGGTA --3’
bta-miR-433	5’--GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACACACCG --3’	5’--GCGATCATGATGGGCTCCT --3’
bta-miR-191	5’--GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACCAGCTG --3’	5’--CGCAACGGAATCCCAAAAG --3’
bta-miR-20a	5’--GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACCTACCT --3’	5’--GCGCGTAAAGTGCTTATAGTGC --3’
4_37054-5p(cli-miR-9632-3p)	5’--GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACAGTCCG --3’	5’--GCGAACTTTTGCCCCTAGTAA --3’
7_42582-5p	5’--GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACTGTTAG --3’	5’--GCGCGTGAGAGGTCTGTAATAA --3’
6_40791-3p	5’--GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACGGTCCC --3’	5’--GCGTGAGTGTGCATCCACG --3’
21_23590-5p(mmu-miR-665-3p)	5’--GTCGTATCCAGTGCAGGGTCCGAGGTATTCGCACTGGATACGACGCTCAT --3’	5’-- CGCCCAGGAGAATGAGGG --3’

3.2. Results

3.2.1. Hormonal Concentration

Hormonal profiles of the Holstein young bull and mature bull testis as the testes developed, the testosterone, estrogen and P4 (progesterone) levels underwent significant changes, with a noteworthy surge from 0.003 ± 0.001 ng/mL to 0.167 ± 0.0078 ng/mL, 62.12 ± 7.07 ng/mL to 57.22 ± 21.41 ng/mL and 0.12 ± 0.01 ng/mL to 0.16± 0.04 ng/mL, respectively, as revealed by hormone testing, indicating a crucial role of hormonal profiles in the developmental stages (Table 3).

3.2.2. Testicular Tissue Morphologic Characteristics of Six Samples

The testicular tissue sections of adult Holstein bull were found to possess seminiferous tubules that were highly vacuolated, with thinner walls as compared to young Holstein bull. The morphological structure of the testes indicated significant difference between the growing young and mature bulls with considerably larger diameter of the seminiferous tubules in sexually matured bulls when observed in 10X and 40 X microscopic magnifications. The interstitial connective tissues of the two groups had remarkably very comparable condition. It can be observed from Figure 1 that the area of the longitudinal section of the testis was expanded in the case of adult Holstein bull. At 40X magnification, more clear structures and more developed stages of spermatogenesis were observed than the young Holstein bulls. It can be seen that there are very few germ cells or spermatocytes present in the seminiferous tubules of young Holstein bull, which is in contrast to Adult Holstein bulls. Moreover, significantly well developed and larger numbers of Sertoli cells were notably seen close to the seminiferous tubules basement (Figure 1).

3.3.2. Quality Control

After high-throughput sequencing filtering the base quality values with FastQC software to remove the 3’ adaptors, the clean reads were then changed to FASTA format indicating the sequencing accuracy. A large number of clean reads, ranging from 83,607,760 to 88,281,454 were expertly generated through the meticulous process of filtering low-quality and covering repeated reads. These reads exhibited a remarkable clean reads rate percentage ranging from 98.17% to 98.33%, and all boasted Q20 and Q30 values that exceeded an impressive 97.98% and 93.77%, respectively (Table 4). After filtering the low quality reads with fastx-toolkit-0.0.14 and with Cutadapt-1.7.1 to remove the 3’ adaptors, the clean reads were displayed with quality control statistics (Table 5).

3.3.3. Reads Mapping and Filtration

Blastn was utilized to precisely align the reads to the Rfam database, which includes Hairpin, repeats, rRNA, snRNA, snoRNA, tRNA and other ncRNAs (Figure 2b). The mapped reads were then filtered using stringent criteria such as gapopen =0, value <0.01, and mismatch ≤1. Subsequently, Bowtie was used to map the exon and intron sequences, where they filtered the reads that could be mapped to the exon sequence but not to the intron sequence, keeping those with mismatch ≤1. Finally, the clean reads were mapped to the bovine reference genome (ARS-UCD1.3) with rigorous mismatch ≤1 criteria (Table 6). The BeoCrystal Program; R language was also employed for the statistical analysis of the clean reads length sequence distribution of each sample, small RNA (sRNA) classification and length distribution (Figure 2a). The clean reads obtained within certain length range from the two groups (X and D group) were further screened for sRNA and analyzed. The range length of the clean reads of sRNA in animals was from 18 to 35, and similarly, in this study also falls between 18 and 30nt. Thus, difference in peaks length distribution helps to the identification of different classes of sRNAs mainly the miRNA, which have predominantly concentrated within the range of 21~22 (Figure 2e-f). A total of 1114 miRNAs were identified from the 8 week and 80 week old Holstein bull testes. The known and novel miRNAs base distribution map with the horizontal and vertical coordinate indicating the number of base and base positions, respectively, for the nucleotides was indicated (Figure 2c-d). In addition, the principal component of these six samples was done (Figure 2k-l), and determined that they could be separated into two distinct groups: X, and D.

3.3.5. The miRNAs Family and Expression Level Analysis

After removing the reads that were not miRNAs, mirDeep2 was utilized to make predictions about novel miRNAs and the secondary structure of miRNAs. Furthermore, quantification of miRNAs and an analysis of miRNA families were conducted in order to gain a deeper understanding of their function in relation to target genes. The counts obtained were then normalized using RPM, which can be visualized through the RPM Volcano plot of the six samples. Using the miRNA expression levels of each sample, the correlation between samples within each group (Figure 4), as well as the distance between the two groups was analyzed and expressed in differential expression cluster heat map (Figure 3) and volcano plot for the up and down regulated known miRNAs (Figure 2i-j). The samples within each group were highly relevant, and that the distance between the two groups was quite noticeable.

3.3.6. Differential Expression Analysis of miRNAs and Target Genes

Through the process of analyzing differential expression, the miRNAs between six samples were examined, revealing unique and distinct patterns of genetic expression. The differential expression of miRNAs among nine samples with CPM >1, p-value < 0.05, |log2FoldChange| >1 was analyzed using the R package edgeR. Subsequently, a heat map was employed to graphically represent the clustering analysis and differentially expressed miRNAs of X vs D. A total of 1114 miRNAs were identified in two age groups (X vs D) with 683 known miRNAs and 431 novel miRNAs. Specifically, 221 known miRNAs were differentially expressed in X vs D, of which 118 were up regulated and 103 were down-regulated. Whereas 237 novel miRNAs were also found differentially expressed in X vs D, of which 77 were up regulated and 160 were down-regulated. The differentially expressed miRNAs were shown in volcano plot (Figure 2i-j). A total of 1641 up and 1506 down regulated target genes with known miRNAs and 1100 up and 2254 down regulated target genes with novel miRNAs were identified based on the complementarity of miRNA sequences and candidate 3’UTR or mRNA sequences as well as the thermal stability of miRNA and mRNA duplexes. From the sequence analysis, 61 known miRNAs were expressed in the young Holstein bulls with 6 known miRNAs significantly different (p<0.05) and 8 known miRNAs expressed in the mature Holstein bull with one known miRNA significantly (p<0.05) expressed. Furthermore, 21 miRNA were expressed in the young Holstein bulls with 6 novel miRNAs significantly different (p<0.05) and 57 novel miRNAs expressed in the mature Holstein bulls, of which 55 novel miRNAs significantly different (p<0.05). Finally, a merged analysis of X vs D was conducted, where in the common up-regulated and down-regulated miRNAs of Holstein bulls were consolidated (Figure 5).

3.3.7. Target Prediction, Functional Enrichment, and Network Analysis

The target genes of differentially expressed miRNAs in X vs D. were predicted with the help of Miranda software. From the 1114 total known and novel miRNAs, 6501 up and down regulated target genes were predicted with miRanda target gene prediction software. The predicted target genes were significantly (p<0.05) participated in in a broad range of biological processes, molecular function and cellular component (Figure 6 and Figure 7). The regulation of mRNA expression by miRNAs was accomplished through patterns of down regulated mRNAs and up regulated miRNAs. The two age groups were examined individually for up-regulated mRNAs and down-regulated miRNAs for further investigation of the phenomena. Gene ontology analysis was conducted on the miRNA targets, and the results were further revealed in GO terms. Additionally, using GO and KEGG pathway analysis, the potential functions of differentially expressed genes in testicular development was investigated (Table 7, Figure 7 and Figure 8). GO enrichment and KEGG functional enrichment analysis was done to provide insight into the DE target gene function in the miRNA-mRNA relationship pairs. From the Known miRNAs, a total of 8131 GO terms were determined from the GO analysis. From the total GO terms analyzed, 6090 terms involved in biological processes, 729 participated in cellular components and 1312 terms belonged to enriched molecular functions. These are significantly (corrected p-value<0.05) enriched processes. While a total of 7951 GO terms were found involved from the novel miRNAs GO term analysis mainly, 5873 GO terms participated in biological processes, 769 in cellular component and 1309 in molecular functions. The significantly differentially expressed miRNAs also enriched in related GO terms for bovine reproduction such as plasma membrane part, plasma membrane region, cell junction, anatomical structure development, cellular developmental process, system development, cell differentiation, and cellular localization (Figure 8 and Figure 9). A total of 458 DE genes, of which 195 up regulated genes and 263 down regulated genes are controlled by miRNAs, were enriched with Go terms associated with the descriptions of the development of the testes and spermatogenic processes.

The miRNA-mRNA correlation network was extensive; therefore hub analysis was performed in order to extract miRNA datasets capable of controlling more than 152 differential genes from the 1295 miRNA-mRNA interaction pairings based on the findings of co-expression research. The differentially expressed genes that are closely associated with the spermatogenic processes are ANKRD49, PAIP2, ALKBH5, STRBP, LGR4, BCL2L11, TYRO3, TBPL1, RGS2, ACVR2A, AGFG1, HMGA2 and DDX25. The STRB with bta-miR-134 and BCL2L11 with bta-let-7e, bta-miR-181b and bta-miR-214 gene were found significantly correlated in the miRNA target gene correlation and network analysis results showing their significant correlation with the spermatogenesis process. The some of the candidate genes which are closely related to reproductive system development includes NOS3, KITLG, MDFI, LRRC6, BAK1, DHCR24, KRT19, BCL2L11, UBE3A, NKX3-1, TYRO3, UBE2Q1, ACVR2A, HMGA2, ADAMTS1, STAT5A, MED1, MAPK14, PRDM1, TTPA, WNT2, FEM1B, MAP2K1, STAT5B, ASCL2, CYP27B1and PDGFRA. The BCL2L11 with bta-let-7e, bta-miR-181b and bta-miR-214, ADAMTSI1 with 17_15431-5p (eca-miR-8912), MED1 with bta-miR-490, bta-miR-433, and 24_29969-5p (has-miR-7155-5p) and PRDM1 with 17_16057-5p were also found strongly linked in the network hub analysis as well as miRNA target gene correlation analysis indicating their involvement in the reproductive system development. More than eight network module results were identified by the hub analysis using Cytoscape, including 17_16057-5p, X_45311-3p (cpo-miR-12072-3p), bta-miR-34c, bta-miR-375, 17_15431-5p (eca-miR-8912), 10_2547-3p (isc-miR-5313), and bta-miR-449b (Figure 6). The DE miRNAs were significantly (p<0.05) enriched in the top-30 KEGG pathways rich-factors including TNF signaling pathways, cell adhesion molecules(CAM), calcium signaling pathway, PI3K-Akt signaling pathway, Ras signaling pathway, ErbB signaling pathway, cGMP-PKG signaling pathway, MAPK signaling pathway, GnRH signaling pathway, VEGF signaling pathway and Rap 1 signaling pathway.

3.3.8. Functional Genes of Differentially Expressed (DE) in the GO Terms

The GO and KEGG-enriched pathways were implemented for the computation of DEGs within the networks. It’s interesting to observe that the DEGs from the young and mature Holstein bull groups shared likely richer GO terms. The GO terms chosen to be most relevant to bovine reproduction were those that dealt with sexual reproduction, sperm part, lateral plasma membrane, androgen secretion, male genitalia development, male gonad development, Sertoli cell development, germ cell development, fusion of sperm and egg, germ cell migration, spermatid development, reproductive process, and the sperm component. It is notable that a considerable number of DE genes significantly enriched in signaling pathways coincide with to the listed GO terms (Table 9).

3.3.9. RT-qPCR Validation

To ensure the consistency of the expression levels of DE miRNAs and their target genes (FOSB, CMC1, SLC2A5, THAP3, and CCNT2), a qPCR validation on 12 miRNAs(bta-miR-11971,bta-miR-2284c, bta-miR-3956, bta-miR-135b, bta-miR-299-2, bta-miR-433, bta-miR-191, bta-miR-20a,4_37054-5p(cli-miR-9632-3p), 7_42582-5p, 6_40791-3p, 21_23590-5p(mmu-miR-665-3p) and its target genes were conducted. Our findings indicate the significance of relative expression of miRNAs compared to the young and adult Holstein bulls. Furthermore, the expression levels of the target genes and miRNAs were found to be consistent with the RNA sequence results (Figure 9 and Figure 10). The Table S2 and Table S4 contain further information on the expression levels of miRNA target genes in the RNA-sequence and RT-qPCR data from this investigation. Supplementary information on the expression levels of miRNAs in the RNA-sequence and RT-qPCR data in this investigation is provided in Table S1 and Table S3.

3.4. Discussion

The development of the reproductive system in male mammals is established primarily during the early stages of testicular development. During this time, the differentiation of somatic cells to Sertoli cells in the early gonad initiates male-specific development and guides germ cells to the spermatogenic lineage (Mäkelä et al., 2019). Thus, testicular development and spermatogenesis are critical processes that impact the reproductive efficiency of male Holstein bulls. Researchers have identified several functional genes that regulate testicular development and spermatogenesis in cattle breeds and yaks ( La et al., 2022; Liu et al., 2022; Xu et al., 2018). Moreover, recent studies have implicated miRNAs in gene regulation by repressing the translation or degradation of target mRNAs (Fabian et al., 2010). Many scholars found that miRNAs are involved in multiple processes of testicular development and spermatogenesis, such as testosterone hormone synthesis, cell proliferation, spermatocyte meiosis, and spermiogenesis of haploid spermatids (Yao et al., 2015). However, no research has been done to identify the miRNAs associated with the early puberty of Holstein bulls’ testicles.

The study selected two age groups with 8 week and 80 week old Holstein bull testes to investigate the expression profiling of miRNAs and mRNA and their contribution to testes development and spermatogenesis. This study is conducted to discover the molecular pathways and mechanisms that control testis development. According to the most recent findings, miRNAs play an important role in early male germ cell development for spermatogenesis, beginning with spermatogonia and progressing to spermatocytes and spermatids, during postnatal male reproduction, mRNA degradation in germ cell proliferation and development via post-transcriptional processes, and expressing meiotic sex chromosome inactivation resistance (Buchold et al., 2010; Gong et al., 2013). A number of studies have found miRNAs to be significant biomarkers in spermatogenesis, testis development, and primordial germ cell proliferation (Gao et al., 2020; Hayashi et al., 2008). Testis development and spermatogenesis, which are governed by numerous functional genes and miRNAs that mainly participated in the spermatogenesis process and testis development, have an important influence on bull reproductive success, sexual maturity and sexual function development (Ding et al., 2016; Liu et al., 2022).

The cDNA libraries of young and mature Holstein bull testes tissues were constructed with clean reads ranging from 98.17% to 98.33%,Q20 exceeding 97.98% and Q30 exceeding 93.77%. The miRNA reads mapping, filtration, differential expression; target prediction and functional enrichment were performed. Furthermore, the expression level of both novel and known miRNAs of the Holstein bull testes tissues was determined through the utilization of the TPM value. Similar methodologies were employed in previous studies for the miRNAs characterization and DE genes identification, such as in the immature and mature testes of indigenous Wandong cattle (Liu et al., 2022), in three ruminant species testes ( Wang et al., 2018), in the cattle-yak male infertility (Zhao, Sun, et al., 2022), in the cattle yak and yak testis (Xu et al., 2018), in the Duroc and Meishan boars (Ding et al., 2016), in the 2-, 6- and 12-month-old small-tail Han-sheep testes (Bai et al., 2019) and in the Saanen Dairy Goat Testis (Wu et al., 2014).

The comparison analysis of miRNA profiles of young (8weeks) and mature (80weeks) of Holstein bull testes were performed to ascertain the potential role of specific miRNAs in the regulation of testicular development and spermatogenesis. A total of 458 miRNAs of which 195 genes significantly up regulated and 263 were observed to be down regulated miRNAs were found differentially expressed in the comparison of X (young) versus D (mature) Holstein bulls. A similar comparison of different pig breed age groups of DE 16 miRNAs with 28 targeted link genes related to spermatogenesis and 14 miRNAs with 18 targeted genes that have a connection to testes development associated with cell adhesion were identified (Li, Yao, et al., 2016). Furthermore, a total of 246 known miRNAs were co-expressed in testicular and ovarian tissues, of which 21 were testis-specific and 9 were ovary-specific, as described by (Huang et al., 2011). Zhao et al. (2022) also investigated the miRNA expression profiles of three ruminant species testis tissues to determine the regulation of spermatogenesis and their regulatory mechanisms involving in male sterility, which suggests our understanding of the physiology of sexual maturity in the testis tissue (Zhao, et al., 2022).

In this study, the Illumina high throughput sequencing technology was used to sequence small RNAs in testicular tissues of 8 weeks and 80 weeks of Holstein bulls. The analyzed DE miRNAs, predicted novel miRNAs, and performed GO enrichment and KEGG pathway analysis of miRNAs and target genes. The DE genes that are enriched in the GO terms in the study for spermatogenesis includes ANKRD49, PAIP2, ALKBH5, STRBP, LGR4, BCL2L11, TYRO3, TBPL1, RGS2, ACVR2A, AGFG1, HMGA2 and DDX25. During the post natal development, the TYRO3 gene is well expressed in Sertoli cells and lack of this gene in male mice results with no production of mature spermatozoa. The mammalian spermatogenesis is also essentially regulated by the Tyro-3 family receptors (Lu et al., 1999). Known and novel miRNAs the Holstein bull’s testes tissues were predicted from cattle (Bos taurus) miRBase data; therefore, more novel miRNAs were determined compared to known DE miRNAs. These newly identified novel miRNAs can serve to enhance our understanding of the crucial miRNAs involved in Holstein bull testis development, and can also provide new information for the cattle (Bos taurus) miRBase.The known DE Holstein bull miRNAs, such as bta-miR-122, bta-miR-146a, bta-miR-184, bta-miR-34c, bta-miR-93, bta-miR-455-3p, bta-miR-369-5p, miR-34a, miR-34b, miR-34c, bta-miR-382 and bta-miR-411c-5p were previously reported in several studies, and have been found to play vital roles in processes such as plasma membrane part, cell cycle, cellular protein modification process, reproductive processes and target protein coding genes regulating associated signaling pathways (Fang et al., 2021; Wang et al., 2018; Zhao, et al., 2022). The expression of many miRNAs such as miR-34 family (miR-34a, miR-34b, and miR-34c) was also reported in other cattle breeds, yak and their interspecies called the dzo (H. Wang et al., 2018). Many authors reported that the miR-34 family involves in an increase of p53 activation in many cells (Cheng et al., 2014; Navarro & Lieberman, 2015) playing vital role in the downstream effector of p53 (Salzman et al., 2016). The miR-34c inhibition also causes apoptosis of male germ cells in murine targeting activation transcription factor 1(atf1) (X. Liang et al., 2012) that involved in cell transformation and survival mediating various extracellular transcriptional signals and responses (Liang et al., 2012). The Ras signaling pathway, HTLV-I infection, cGMP-PKG signaling pathway, EGFR tyrosine kinase inhibitor resistance, apoptosis - multiple species and oxytocin signaling pathways was assessed in relation to DE source genes associated with spermatogenesis such as BCL2L11, TBPL1, RGS2, ACVR2A, HMGA2 and determined to be linked to bovine male reproduction.

The results of this study may shed light on the molecular mechanisms that regulate testicular development and spermatogenesis in male Holstein bulls. Furthermore, this research could help to identify potential therapeutic targets for improving reproductive efficiency in male Holstein bulls. Ultimately, this study could have significant implications for the Holstein bulls breeding industry.

3.5. Conclusions

This study represents the first comprehensive analysis of miRNA and mRNA expression patterns in the developmental process of Holstein bull testes. Through this investigation, key miRNAs and miRNA-targets in three 8 weeks young Holstein bull testes (X) vs three 80 weeks old mature Holstein bull testes (D) were successfully explored. The insights gained from this study will serve as a solid foundation for elucidating the molecular regulatory mechanisms underlying Holstein bull testes development, as well as identifying biomarkers that can significantly affect the reproductive efficiency of Holstein bull testes.