Transcriptomic Analysis of the Developing Testis and Spermatogenesis in Qianbei Ma Goats

Reproductive competence in male mammals depends on testicular function. Testicular development and spermatogenesis in goats involve highly complex physiological processes. In this study, six testes were, respectively, obtained from each age group, immature (1 month), sexually mature (6 months) and physically mature (12 months old) Qianbei Ma goats. RNA-Seq was performed to assess testicular mRNA expression in Qianbei Ma goats at different developmental stages. Totally, 18 libraries were constructed to screen genes and pathways involved in testis development and spermatogenesis. Totally, 9724 upregulated and 4153 downregulated DEGs were found between immature (I) and sexually mature (S) samples; 7 upregulated and 3 downregulated DEGs were found between sexually mature (S) and physically mature (P) samples, and about 4% of the DEGs underwent alternative splicing events between I and S. Select genes were assessed by qRT-PCR, corroborating RNA-Seq findings. The detected genes have key roles in multiple developmental stages of goat testicular development and spermatogenesis. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed to determine differentially expressed genes (DEGs). GO analysis revealed DEGs between S and P contributed to “reproduction process”, “channel activity” and “cell periphery part” between I and S, and in “ion transport process”, “channel activity” and “transporter complex part”. KEGG analysis suggested the involvement of “glycerolipid metabolism”, “steroid hormone biosynthesis” and “MAPK signaling pathway” in testis development and spermatogenesis. Genes including IGF1, TGFB1, TGFBR1 and EGFR may control the development of the testis from immature to sexually mature, which might be important candidate genes for the development of goat testis. The current study provides novel insights into goat testicular development and spermatogenesis.


Introduction
The Qianbei Ma goat is a unique goat breed in the mountainous area of Guizhou Plateau. It is large in size, strong in constitution and delicious in taste. It is a goat breed whose skin and meat are both valuable, and has many advantages such as strong adaptability, high fertility, stable genetic inheritance and strong disease resistance. However, the research on its reproductive performance mainly focuses on lambing and estrus in she-goat and the molecular biology understanding of testicular development or spermatogenesis in Qianbei Ma goat is very limited, deserving further investigation. The breeding performance of male animals is also important to the total fertility of the population. Additionally, the most important guarantee for the reproductive performance of rams is the normal development of testis and spermatogenesis. Testicular development and spermatogenesis are very complex and precise processes which are regulated by multiple genes.
The testis represents a critical organ of the male reproductive system in mammals, producing spermatozoids and androgens. Spermatogenesis constitutes a developmental event, which produces haploid spermatozoa from diploid spermatogonial stem cells during meiosis in the testis [1]. Spermatogenesis comprises three steps, including spermatogonial proliferation and differentiation, meiotic division of spermatocytes and spermatozoid maturation [2]. Aside from spermatogenic cells, spermatogenic mechanisms involve many somatic cells in the testicle, including Sertoli and Leydig cells. In the testis of mammals, spermatogenesis occurs in seminiferous tubules, where germ cells are associated with Sertoli cells. Genes associated with the latter cells also have critical functions in precise steps of spermatogenesis [3]. The proper development of the testicles ensures the spermatogenesis process.
RNA sequencing (RNA-Seq) allows an expression profiling of genes and might help map and quantify the transcriptome [4,5]. This approach offers many advantages compared to other transcriptomic tools, including high resolution/sensitivity, a broad dynamic range of gene expression, and the identification of new transcript sequences and splice isoforms of previously reported genes [6]. Ramsköld [7] evaluated multiple tissues from mammals by RNA-Seq and reported that most genes were specifically expressed in testicular samples. Additionally, considering RNA-Seq-based expression patterns, Djureinovic [8] categorized 20,050 putative human genes, which showed specific expression in the human testicle. Their evaluation revealed the testicular tissue had by far the largest quantity of tissue-specific genes. Using microarray analysis, Anand detected differentially expressed genes (DEGs) in testicular samples compared to other tissues, identifying 2868 upregulated transcripts and 2011 downregulated mRNAs [9]. The testicle appears to have a higher degree of metabolic activity relative to other normal tissues. Most current reports assessing the association of testicular development with spermatogenesis have been performed in the human or mouse species, and the goat is scarcely examined.
The present work aimed to perform transcriptome profiling of immature and mature Qianbei Ma goat testis specimens by RNA-Seq and bioinformatic analysis. The findings provide novel insights into the mechanisms regulating goat testicular development and spermatogenesis.

Animal Handling and Sample Collection
Permission was granted to castrate eighteen healthy Qianbei Ma goats in Fuxing Husbandry Co., Ltd., Zunyi, Guizhou, China. Goat ages were obtained from goat farming records. There were six immature goats (1 month old, before sexual maturation, i.e., samples I1, I2, I3, I4, I5 and I6), six sexually mature goats (6 months old, after sexual maturation but before physical maturation, i.e., samples S1, S2, S3, S4, S5 and S6) and six physically mature goats (12 months old, after physical maturation, i.e., samples P1, P2, P3, P4, P5 and P6). We surgically collected the right testes from the eighteen goats by castration after anesthesia, followed by storage in RNA/DNA sample protector (Servicebio, Wuhan, China). The testis from each goat was cut longitudinally, and a small amount (3-5 g) of the parenchyma, including seminiferous tubules and Leydig cells, underwent snap freezing in liquid nitrogen and was transported to the lab for further studies. All castrated goats remained in Fuxing Husbandry Co., Ltd. (Guizhou, China) after our study, for fatten feeding.

RNA Quantitation and Quality
Total RNA extraction was carried out from the testicular tissue in groups I, S and P using TRIzol reagent (Servicebio, Wuhan, China) and RNeasy RNA purification kit (Servicebio, Wuhan, China) with DNase as directed by the manufacturer. A NanoDrop™ One spectrophotometer (Thermo Fisher Scientific, Waltham, Ma, USA) was utilized to assess RNA purity and amounts. RNA quality assessment utilized 1% agarose gel electrophoresis. High-quality RNA samples (OD 260/280 of 1.8-2.0, integrity > 7.0 and 28S:18S above 1.0) were sequenced on an Illumina NovaSeq 6000 system, generating 150-bp paired end reads.

Transcriptome Sequencing
Approximately 5 µg RNA/sample constituted the input material for RNA sample preparation. Index-coded samples were clustered with NEBNext ® Ultra™ Directional RNA Library Prep Kit for Illumina ® according to a protocol provided by the manufacturer. Upon clustering, the prepared libraries underwent sequencing on an Illumina NovaSeq 6000 (Illumina, San Francisco, CA, USA). The image data of the sequences yielded by the high-throughput sequencer were converted into sequence data (reads) by CASA V A base recognition to obtain FASTQ files. Raw RNA-Seq FASTQ data next underwent filtration with Fastp v [10], a quality control software that can quickly filter and correct FASTQ data to exclude adapter-containing, N-containing and low-quality (quality score below 20) reads, resulting in clean reads, and were mapped to the goat (C. hircus) (ARS1.2) reference genome [4] using HISAT2, a software uses a graph Ferragina Manzini index to align DNA and RNA sequences [11].

Quantification of Gene Expression
Reads mapped to a given gene were counted with featurerts for estimating the expression of various gene transcripts. Gene expression was determined from million mapped reads per kilobase (FPKM) values [12], the current commonest approach to estimate gene expression [13].

Differential Expression Analysis
The DESeq2 software (1.20.0) was used to analyze differential expression between the treatment and control groups. The Benjamini-Hochberg algorithm was utilized to adjust p values (p-adj) to control for false discovery rate. |log2 (FoldChange)| ≥ 1 and padj < 0.05 was set as the significance threshold for differential expression [14].

GO and KEGG Enrichment Analyses of DEGs
GO and KEGG analyses of DEGs were implemented with ClusterProfiler, correcting for gene length bias. KEGG is an information database based on molecular findings, particularly via genome sequencing and additional high-throughput techniques to produce large-scale molecular data sets, allowing a deep understanding of biological systems (http://www.genome.jp/kegg/, accessed on 21 July 2022) [15]. GO and KEGG terms with |log2 (FoldChange)| ≥1 and padj < 0.05 were deemed to be DEGs with significant enrichment [10].

Prediction of New Transcripts and Alternative Splicing Analysis
StringTie was utilized to build and identify previously reported and new transcripts from HISAT2 alignment data. StringTie utilizes a network-flow algorithm with optional de novo assembly to splice transcripts. Compared with cufflinks, StringTie has the following classifying AS events, which were assessed in various samples separately.

Gene Expression Profiling during Testicular Development in Qianbei Ma Goats
To determine the genes associated with testicular development and sperm formation, 18 libraries, including 6 each from immature, sexually mature and physically mature testes were sequenced. Table 2 shows an overview of the information pertaining to raw and clean reads for all libraries. Error rates and GC contents for various libraries were determined for their quality control. All 18 libraries were of high quality.
Paired-end clean reads underwent alignment to the goat reference genome ARS1.2 that is commonly utilized for C. hircus (goats) [4] with HISAT2, and the expected number of fragments per kilobase of transcript sequence per million base pairs sequenced (FPKM) of each gene was calculated. Table 3 shows an overview of the information about uniquely mapped clean reads.
Additionally, an overview of the percentages of clean reads mapped to exon borders (junction reads) is shown in Table 4.
All RNA-Seq data in our study had been submitted to NCBI (accession number BioProject: PRJNA879963).

Alternative Splicing Data
In this work, AS events were categorized into five types with rMATS (http://rnaseqmats.sourceforge.net/index.html, accessed on 21 July 2022). By determining the types and amounts of AS events, and analyzing each AS type, a large number of AS events were found in testicular development and sperm formation. Of all the DEGs detected between immature (I) and sexually mature (S) testes in our present work, 544 genes underwent AS events between sexually mature (S) and physically mature (P) testes, with no detected AS events (padj ≤ 0.05). Therefore, approximately 4% of genes showed AS between I and S, and no gene had AS between S and P in this study. The above findings indicated AS was very important in the complexity of gene expression during testis development, especially in the period from immaturity to sexual maturity. Error rate: overall sequencing error rate for the data; Q20 and Q30: percentages of total bases with Phred values above 20 and 30, respectively; GC pct: percentage of C and G among the 4 bases in clean reads.  To determine AS types associated with testicular genes, AS events were compared between the S group and the I group. The five known types of AS events include retained intron (RI), mutually exclusive exon (MXE), alternative 3 splice site (A3SS), alternative 5 splice site (A5SS) and skipped exon (SE). All five types of AS were found in the S vs. I group comparison. RI, MXE, A3SS and SE were found in the P vs. S group comparison. These findings showed SE as the most common AS event amounting to 368 in S vs. I. Other identified AS events were RI (48)

Analysis of DEGs
In order to clarify the difference in gene expression in testis of Qianbei Ma goat at different developmental stages, we compared and analyzed the differentially expressed genes in three different developmental stages. DEGs were determined with |log2 (fold−change)| ≥ 1 and padj < 0.05. As a result, a total of 9,724 upregulated and 4,153 downregulated DEGs were detected between the I and S groups; 7 upregulated and 3 downregulated genes were detected between the S and P groups ( Figure 2A1, Figure 2A2). Among all DEGs detected, 1,112 genes were only expressed in the immature group, 379 genes were only expressed in the sexual mature group, 336 genes were only expressed in the physical mature group, and 13,299 genes were expressed in all three groups ( Figure  2B). Figure 2C depicts hierarchical clustering, with DEGs for the eighteen libraries grouped into three clusters, the result suggested the I, S and P groups had differential expression patterns overall, but S and P had similar repetitive expression commonalities.

Analysis of DEGs
In order to clarify the difference in gene expression in testis of Qianbei Ma goat at different developmental stages, we compared and analyzed the differentially expressed genes in three different developmental stages. DEGs were determined with |log2 (fold−change)| ≥ 1 and padj < 0.05. As a result, a total of 9,724 upregulated and 4,153 downregulated DEGs were detected between the I and S groups; 7 upregulated and 3 downregulated genes were detected between the S and P groups ( Figure 2A1, Figure 2A2). Among all DEGs detected, 1112 genes were only expressed in the immature group, 379 genes were only expressed in the sexual mature group, 336 genes were only expressed in the physical mature group, and 13,299 genes were expressed in all three groups ( Figure 2B). Figure 2C depicts hierarchical clustering, with DEGs for the eighteen libraries grouped into three clusters, the result suggested the I, S and P groups had differential expression patterns overall, but S and P had similar repetitive expression commonalities.

GO Analysis of DEGs
GO analysis was carried out to examine the functions of DEGs in testicular development. Totally, 106 GO terms associated with "biological processes", "molecular functions" and "cellular components" were markedly enriched between immature (I) and sexually mature (S) testes. We found 43 GO terms belonged to "biological process", 56 GO terms belonged to "molecular functions", and 7 GO terms belonged to "cellular components". Totally, 36 GO terms associated with "biological processes", "molecular functions" and

GO Analysis of DEGs
GO analysis was carried out to examine the functions of DEGs in testicular development. Totally, 106 GO terms associated with "biological processes", "molecular functions" and "cellular components" were markedly enriched between immature (I) and sexually mature (S) testes. We found 43 GO terms belonged to "biological process", 56 GO terms belonged to "molecular functions", and 7 GO terms belonged to "cellular components". Totally, 36 GO terms associated with "biological processes", "molecular functions" and "cellular components" had significant enrichment in the sexually mature (S) vs. physically mature (P) testes. We found 5 GO terms belonged to "biological processes", 21 GO terms belonged to "molecular functions", and 10 GO terms belonged to "cellular components" (Table 5, Figure 3A,B). Genes 2023, 14, x FOR PEER REVIEW 9 of 16 "cellular components" had significant enrichment in the sexually mature (S) vs. physically mature (P) testes. We found 5 GO terms belonged to "biological processes", 21 GO terms belonged to "molecular functions", and 10 GO terms belonged to "cellular components" (Table 5, Figure 3A, Figure 3B).  Most enriched GO terms between sexually mature (S) and physically mature (P) testes. Abscissas and ordinates represent enriched GO terms and significance levels of GO enrichment, respectively. Purple, "biological processes"; red, "cellular components"; blue, "molecular functions".

qRT-PCR Validation of DEGs
To verify the DEGs in immature (I), sexually mature (S) and physically mature testes, we selected 6 DEGs, including TGFBR1, TGFB1, EGFR, IGF1, MAPK3 and SMA to validate RNA-Seq data by qRT-PCR. qRT-PCR data corroborated RNA-Seq findi suggesting the reliability of RNA-Seq data ( Figure 5).

Discussion
With the current development of detection technology, more and more mRNAs on testis related to testis development and spermatogenesis have been reported. RNA-Seq has emerged as a tool for efficiently and inexpensively detecting new transcripts and genes. RNA-Seq methods have been broadly utilized to determine DEGs or gene expression patterns, new transcripts, AS events and SNPS, and have empowered studies examining porcine [24,25], cattle [26,27] and mouse [28,29] testicular development. In goats, the profiles of ovarian [30,31], uterine [32,33] and testicular [34,35] tissues under different conditions were recently compared by RNA-Seq. However, limited data on testicular development in goats were available. Breed and age represent major factors affecting testicular development. Here, RNA-Seq was performed to build a complete dataset that explains the spatiotemporal transcriptome of the testicular tissue in Qianbei Ma goats. Testicular growth and development constitute the key factors affecting goat reproduction. Therefore, identifying genes regulating testicular growth and development is critical. In this study, 13,887 genes were assessed by RNA-Seq in six immature, six sexually mature and six physically mature testes. Totally, 9724 genes were upregulated and 4153 were downregulated between immature and sexually mature testes; seven genes were upregulated and three were downregulated between sexually mature and physically mature testes which was similar to the hierarchical clustering analysis result in our study. The reason for this result is related to the development process of the testis of Qianbei Ma goat. We speculate that the testis develops from immature to sexual mature with various reproductive functions and also starts to produce sperm, which is a very complex process that requires the co-regulation of a large number of different genes. When the goats continue to develop from sexual maturity to physical maturity, the testes already had all the functions related to reproduction, so the differences in development were relatively small, the number of related DEGs was significantly reduced, which was similar to the results of Gong

Discussion
With the current development of detection technology, more and more mRNAs on testis related to testis development and spermatogenesis have been reported. RNA-Seq has emerged as a tool for efficiently and inexpensively detecting new transcripts and genes. RNA-Seq methods have been broadly utilized to determine DEGs or gene expression patterns, new transcripts, AS events and SNPS, and have empowered studies examining porcine [24,25], cattle [26,27] and mouse [28,29] testicular development. In goats, the profiles of ovarian [30,31], uterine [32,33] and testicular [34,35] tissues under different conditions were recently compared by RNA-Seq. However, limited data on testicular development in goats were available. Breed and age represent major factors affecting testicular development. Here, RNA-Seq was performed to build a complete dataset that explains the spatiotemporal transcriptome of the testicular tissue in Qianbei Ma goats. Testicular growth and development constitute the key factors affecting goat reproduction. Therefore, identifying genes regulating testicular growth and development is critical. In this study, 13,887 genes were assessed by RNA-Seq in six immature, six sexually mature and six physically mature testes. Totally, 9724 genes were upregulated and 4153 were downregulated between immature and sexually mature testes; seven genes were upregulated and three were downregulated between sexually mature and physically mature testes which was similar to the hierarchical clustering analysis result in our study. The reason for this result is related to the development process of the testis of Qianbei Ma goat. We speculate that the testis develops from immature to sexual mature with various reproductive functions and also starts to produce sperm, which is a very complex process that requires the co-regulation of a large number of different genes. When the goats continue to develop from sexual maturity to physical maturity, the testes already had all the functions related to reproduction, so the differences in development were relatively small, the number of related DEGs was significantly reduced, which was similar to the results of Gong studies on testis of mice at three different developmental stages [36]. Using next-generation platforms, we determined most upregulated genes were associated with protein coding and might have functions in testicular development and sperm formation.
Alternative splicing (AS) represents a commonly encountered phenomenon in eukaryotes, which could lead to the production of various protein forms at different times under different circumstances, increasing species/body fitness. Although AS research is a known subfield of molecular biology, only in recent years has this subfield attracted sufficient attention. AS is critical for the complex proteomes and functions found in higher organisms. AS is an important mechanism in the regulation of gene expression and promotes proteome diversity [37]. It is estimated that about 95% of human multiple-exon gene expression is associated with AS events [38]. In metazoans, AS is critical for the production of various protein forms with functions in different cell events such as cell growth, differentiation and death [39]. Here, five AS events were observed, mostly involving ES. The effects of AS events on the functions of related genes can be predicted by a comprehensive analysis of AS events, and GO and KEGG analysis data [40,41]. In the current study, Sec insertion sequence binding protein 2 (SECISBP2) was the gene with the highest number of SE events, i.e., a total of ten SE events. SECISBP2 underwent the most AS events during the development, indicating that SECISBP2 protein synthesis is complex and important to Qianbei Ma goats' testis development. Mutation of the gene altered thyroid hormone metabolism [42,43]. Thyroid hormones can modulate semen quality under physiological conditions by regulating testosterone and changing some semen indexes [44][45][46]. Thyroid hormone levels in the early stages of testicular development can also influence testicular development and spermatogenesis by regulating the length of time during which supporting cells proliferate [47]. Based on previous studies, we hypothesized that there is active synthesis of SECISBP2 during testicular development in Qianbei Ma goat, which is essential for testicular size and spermatogenesis.
Combining previous relevant reports, and KEGG and GO data in the current study, the genes involved in the regulation of testis development and sperm formation through protein phosphorylation were mainly TGFB1, EGFR and IGF1, which have critical functions in testis growth, hormone secretion, spermatogenesis and Leydig cell differentiation.
Transforming growth factor β-1 (TGFB1) plays multiple biological roles, including the control of proliferative and differentiation potentials of cells in rams [48]. In male lambs, TGFB1 regulates tight junctions in Sertoli cells and controls spermatogenesis. It modifies the blood-testicular barrier (BTB) by downregulating tight junction proteins [49]. TGFB1 might play an important role in testicular development because of its high expression in the immature testis and markedly reduced expression in sexual maturity, as spermatogenesis begins. A comparable expression pattern was found for TGFB receptor type 1 (TGFBR1) [50]. In addition, loss of TGFB1 resulted in lower testosterone levels in the testis and serum, and decreased the ability to mate with females [51]. Based on previous studies, we speculated that TGFB1 may not only directly regulate goat testicular development and sperm formation, but also ensure the normal development of male external genitalia and affect fertility.
Epidermal growth factor receptor (EGFR) represents a receptor gene for EGF and controls testicular function in the human, mouse, rat and livestock species as well as in alpacas [20]. In yak and bovine, EGF and EGFR were critical paracrine and/or autocrine modulators of testicular development and sperm formation and regulate testosterone production by testicular interstitial cells [52,53]. We speculate that EGF and EGFR may also be expressed in various goat testicular cells, stimulate testosterone secretion, and regulate testis development and spermatogenesis.
Insulin-like growth factor I (IGF1) contributes to the regulation of testicular function [54]. Pitetti et al. indicated growth factors of the insulin family play essential roles by controlling SC number, testis size and daily sperm production [55,56]. Both IGF1 and its receptor IGF1R were expressed in testicles, and their hormones act directly on male gonads [57,58]. In immature testes, IGF1 promotes the development of sustentacular cells, Leydig cells and gonocytes. In mature testis, the IGF1 gene induces spermatogenesis and regulates Leydig cell function [59,60]. IGF1 may act as an autocrine/paracrine or endocrine signal to regulate testicular steroid production as well as germ cell and Sertoli cell functions [61]. IGF1 plays different roles in testicular function at different stages of testicular development [54,62]. We speculate that high IGF1 and IGF1R protein amounts in the immature testis may suggest they highly promote the development and differentiation of sustentacular cells, Leydig cells and gonocytes in goat testis during sexual maturity.
Transcriptome data revealed the MAPK pathway is implicated in goat testicular development, while TGFBR1, TGFB1, EGFR and IGF1 were enriched in this pathway and downregulated during sexual maturation, as key genes that regulate testis development and spermatogenesis [63,64]. Multiple reports suggest MAPK signaling is a critical regulator of testis growth and development, testis cell proliferation, differentiation and apoptosis and testosterone secretion, thus affecting male fertility [65][66][67]. One of the key downstream target genes of MAPK signaling is IGF1, which together with other genes in the pathway controls testis cell proliferation, testis volume development, hormone secretion and spermatogenesis, and is often reported to be associated with male fertility [55,68,69]. We speculate that MAPK signaling is a critical regulatory pathway in goat testis development and spermatogenesis. As an essential male fecundity-related gene in the MAPK signaling pathway, IGF1 modulates goat testis growth and development and can affect the functions of various cells of goat testis by regulating other downstream genes in the signaling pathway. During male goat sexual maturation, IGF1 regulates the development of the testis, spermatogenesis and hormone synthesis through the MAPK signaling pathway and other cooperative genes.

Conclusions
This study first used RNA-Seq to profile the expression of genes during testicular development in Qianbei Ma goats. We identified 544 genes with AS events between I and S, which suggests that AS of differential genes might be critical in the regulation of testicular development in goats. Totally, eight KEGG pathways were upregulated and 90 were downregulated between I and S. Totally, seven KEGG pathways were upregulated between S and P. Among the six screened DEGs (TGFBR1, TGFB1, EGFR, IGF1, MAPK3 and SMAD4) TGFBR1, TGFB1, EGFR, IGF1 and MAPK3 belonged to the "MAPK signaling pathway", corresponding to the GO term "protein phosphorylation". However, the specific function and pathway mechanism of the selected genes still need to be further studied at the cellular level. The findings should further our understanding of gene regulation during testicular development and sperm formation. Our study also suggests that RNA-seq has great potential in the study of the transcriptome during the development of goat testis.

Institutional Review Board Statement:
The animals used in our study were consented by their owners. All procedures involving animals were approved and authorized by Guizhou University. The laboratory animal castration protocol for this study was approved by the Laboratory Animal Ethics of Guizhou University (No. EAE-GZU-2021-P024, Guiyang, China; 30 March 2021).

Informed Consent Statement: Not applicable.
Data Availability Statement: The raw data from sequencing are available at NCBI under BioProject ID: PRJNA879963.