Abstract Larimichthys crocea is an economically important marine fish in China. To date, the molecular mechanisms underlying testicular development and spermatogenesis in L. crocea have not been thoroughly elucidated. In this study, we conducted a comparative transcriptome analysis between testes (TES) and pooled multiple tissues (PMT) (liver, spleen, heart, and kidney) from six male individuals. More than 54 million clean reads were yielded from TES and PMT libraries. After mapping to the draft genome of L. crocea, we acquired 25,787 genes from the transcriptome dataset. Expression analyses identified a total of 3853 differentially expressed genes (DEGs), including 2194 testes-biased genes (highly expressed in the TES) and 1659 somatic-biased genes (highly expressed in the PMT). The dataset was further annotated by blasting with multi-databases. Functional genes and enrichment pathways involved in spermatogenesis and testicular development were analyzed, such as the neuroactive ligand–receptor interaction pathway, gonadotropin-releasing hormone (GnRH) and mitogen-activated protein kinase (MAPK) signaling pathways, cell cycle pathway, and dynein, kinesin, myosin, actin, heat shock protein (hsp), synaptonemal complex protein 2 (sycp2), doublesex- and mab-3-related transcription factor 1 (dmrt1), spermatogenesis-associated genes (spata), DEAD-Box Helicases (ddx), tudor domain-containing protein (tdrd), and piwi genes. The candidate genes identified by this study lay the foundation for further studies into the molecular mechanisms underlying testicular development and spermatogenesis in L. crocea. Keywords: Larimichthys crocea, RNA-Seq, testicular development, spermatogenesis 1. Introduction Large yellow croaker (Larimichthys crocea) is an economically important maricultured species in southeastern China, and it has received considerable attention because of its high commercial and nutritional value. Due to breakthroughs in artificial propagation technology for L. crocea in 1990, its aquaculture production has increased rapidly. In 2017, the annual output of cultured L. crocea in China was approximately 176,600 tons, which ranked first among maricultured fishes [[38]1]. However, with the rapid development of the culture industry, marine cage-cultured L. crocea have suffered from large-scale sexual precocity, which is extremely detrimental to the normal growth of fish [[39]2]. An increasing number of L. crocea mature when small-sized, and once the fish are sexually mature, energy and nutrients are mostly shifted to the processes of gonad maturation or reproduction, and only a small amount remains for somatic growth; these changes eventually reduce the culture benefit of L. crocea [[40]3] and hinder the conservation of its germplasm resources. Therefore, the molecular mechanisms underlying reproduction regulation must be investigated, and the genes involved in gametogenesis and gonadal development must be identified, to provide a theoretical reference for further studies on the above problems in L. crocea. Directly correlated with sexual maturation and reproduction, in various animal species, including L. crocea, the testis is an essential basic component of the animal reproductive system and responsible for the production of male gametes via spermatogenesis [[41]4]. During the process of spermatogenesis, diploid spermatogonia slowly evolve into many highly specialized spermatozoa through mitosis, meiosis, and spermiogenesis. Stringent spatial and temporal expression of genes during both transcriptional and translational processes are fundamental to ensure the precise processes of spermatogenesis [[42]5]. Previous studies have focused on the reproductive and developmental biology of L. crocea based on anatomical/histological aspects [[43]6,[44]7], whereas few studies have focused on genes related to spermatogenesis and testicular development. To the best of our knowledge, only a small number of reproduction-related genes have been identified in L. crocea [[45]8,[46]9,[47]10,[48]11]. To further comprehensively explore the molecular mechanisms underlying spermatogenesis and testicular development in L. crocea, large-scale genetic resources must be identified and studied. Recently, RNA-Seq technology has been widely used to explain biological processes in almost all organisms. Over the past three years, RNA-Seq analyses of reproductive tissues have been performed in several commercial teleosts, such as Thalassoma bifasciatum [[49]12], Takifugu rubripes [[50]13], and Ictalurus punctatus [[51]14]. Numerous reproductive-related genes and signaling pathways have been screened, which has deepened our understanding of gonadal development and gametogenesis. In this study, we identified the genes and pathways potentially participating in testis development and spermatogenesis in L. crocea. This study could offer fundamental information for further research on the mechanisms underlying reproductive development in male L. crocea. 2. Materials and Methods 2.1. Ethics Statement All fish handling and experimental procedures were approved by the Animal Care and Use Committee of Ningbo University. Prior to sampling, all fish were anesthetized with 0.05% MS-222 (3-Aminobenzoic acid ethyl ester methanesulfonate, Sigma, Saint Louis, MO, USA). 2.2. Sample Collection, Histological Identification, RNA Isolation, and cDNA Library Preparation Healthy male fish (17 ± 0.8 cm, 59 ± 4.1 g) were commercially obtained and sampled from a local farm in May 2016. After 2 weeks of acclimatization in seawater (DO 8.6 ± 0.3 mg L^−1, 23 ± 1 °C), fish were dissected on ice and the tissues (testis, liver, spleen, heart, and kidney) from 20 fish were removed, frozen in liquid nitrogen, and stored at −80 °C. Concurrently, a small part of the testis of each fish was used to perform histological analysis so that we could select individuals with testes at developmental stage IV, in which the process of spermatogenesis was active. The testicular developmental stages were identified based on histological features, as mentioned in a previous report [[52]15]. After validating the testes development, 6 male individuals (stage IV) of the 20 fish were chosen for follow-up experiments. Total RNA was isolated with TRIzol reagent (Invitrogen, Shanghai, China) according to a standard protocol. RNA integrity, purity, and concentration were validated as described by Zhan [[53]16] and Chen [[54]17]. Only samples with an RNA Integrity Number (RIN) > 8.0 were used for RNA library construction. The MixS (minimum information about any (x) sequence) descriptors are presented in [55]Table 1. Table 1. MIxS information for the transcriptome of L. crocea. Item Description Investigation type Eukaryote Project name Transcriptome for L. crocea Collection date May 2016 Lat_lon 29°86′ N, 121°56′ E Geo loc name Ningbo, China Environment Marine water Biotic_relationship Free-living Trophic level Heterotroph Temp 21 °C Salinity 24 PSU Estimated size 16.35 Gb Sequencing meth Illumina HiSeq^TM 2500 Mapping software TopHat2 Annotation source Nr/Nt/Swiss-Prot/KEGG/COG/GO BioProject ID of raw reads PRJNA471154 for TES/PRJNA471574 for PMT Accession number of raw reads SRP148410 for TES/SRP148493 for PMT [56]Open in a new tab To normalize individual differences, an equal amount of RNA isolated from the testes of the six male individuals were pooled and constructed a TES cDNA library. Similarly, an equal amount of RNA isolated from other somatic tissues (liver, spleen, heart, and kidney) were also pooled and constructed a PMT cDNA library. Sequencing library construction and Illumina sequencing were conducted at Beijing BioMarker Technologies (Beijing, China) following the manufacturer’s recommendations. The complete RNA-Seq process and bioinformatics data analysis workflow are presented in [57]Figure 1. Figure 1. [58]Figure 1 [59]Open in a new tab RNA-Sequencing data processing and analysis workflow. 2.3. Sequencing and Mapping to the Reference Genome Clean data (clean reads) were obtained by removing reads containing adapters, ploy-N, and inferior reads from the raw data. Then, we used the L. crocea genome (accession number: PRJNA354443) as a reference for sequence alignment with the Tophat2 software. HTSeq v0.6.1 with −m union was used to count the number of reads mapped to each gene. For normalization, the count for each gene was divided by the number of fragments per kilobase of transcript sequence per million base pairs (FPKM) sequenced in each sample. FPKM were used to estimate the gene expression level. 2.4. Differential Expression Analysis The differential expression analysis of two samples was performed using the EBseq R package. p-values were adjusted using the Q value [[60]18], and Q value < 0.01 and |log[2] (fold change)| ≥ 2 were set as the thresholds for significant differential expression. 2.5. Functional Annotation and Pathway Enrichment Analysis The transcripts sequences in our study were aligned by BLASTx to the Nr, GO, Swiss-Prot, KEGG, and COG databases (E-value < 10^−5) and aligned by BLASTn to the Nt database (E-value < 10^−5). Candidate DEGs potentially associated with spermatogenesis and testicular development in L. crocea were identified by the combination of enrichment analyses, annotations, and manual literature searches. 2.6. Quantitative Real-Time PCR (qPCR) Validation A total of 15 DEGs with significantly different expression in the TES and PMT were randomly selected for validation of the Illumina sequencing data via a qPCR (LightCycler 480 Roche) analysis. cDNA was synthesized using the ReverTra Ace qPCR RT Master Mix with gDNA Remover (TOYOBO BIOTECH CO., LTD. Shanghai, China) with the total RNA used in the RNA-Seq analysis. The primers designed for amplification were based on RNA-Seq data in our library and NCBI nucleotide databases. All primer sequences are listed in [61]Table 2. The reaction was performed in a 20 μL reaction volume containing 4 μL cDNA, 10 μL 2×RealStar Green Fast Mixture (Genstar, China), 1 μL each primer (10 mM), and 4 μL nuclease-free water. β-actin was used as an internal control to normalize the gene expression level. The PCR conditions were as follows: 10 min at 95 °C, 40 cycles at 95 °C for 30 s, 60 °C for 30 s, 72 °C for 30 s, and followed by a cooling step at 4 °C. Every sample was amplified in triplicate to normalize the system and reduce pipetting errors. The 2^−ΔΔCT method was used to analyze the relative expression level. Pearson’s correlation coefficient was used to assess the expression data consistency between the RNA-Seq and qPCR. Table 2. Primers used for qPCR. Gene ID Annotation Primer Sequence (5′→3′) gene15399 tex9 F GCTGTAGACGACTCGGCTGACTT R TGAGCATCTGAAACGCCTGATCCT gene14348 tex36 F GCAAGGAGTTGTCACACTGGCATC R TCCATCGTGGCACAGGCAGAAG gene2626 tex11 F TCGGTGAAGTCTCTGCTCTGGAAG R GGACGCCCTCTGTTGGATTCTCA gene21265 tcte1 F TGATCGGAGACAGAGGAGCCAGAG R CGGTTGAGACGCAGGTTGAGTGA gene20395 stpg2 F GCAGTCTCCAGAACCGCTCCAA R CAGTGTGTCCTCCTCGTAGCCAAA gene3104 spata7 F CTGAGGATGAGTCCAACGGCACAT R AGATTTCCCGCCTTCTGGTGAGAC gene26146 spata20 F GTTCCTGGACGACTACGCCTTCAT R TGGACGCCGACACTGAGTTAGC gene26767 spag5 F GGACATCCAGCAAGCCAATGACAG R CCTCGCCAACTCGTTCTCCATCT gene26247 spag17 F CCAGACGAGGAGGAGGACAGAGAA R TTCAGGATGGTGATGCCGAACTCA gene16701 pabp1l F AGTCCGCTAATGGAGGCTCTGTC R AGTGGTGGTCCTTGTGGTTGATGT gene10808 nme5h F TCCACGGCAGCGAGTCATTTCAT R TCAGCCAGTCAGCAAGCCAGATAC gene21216 meig1 F ACAACTCCAAGCCGAAGTCCATGA R TTGACATCACGGTCCAGGCACTC gene13149 lhcgrl F TGGCATCAAGGAGGTGGCAAGT R TGGTAGGCGGACTCTGCGATCA gene1213 hsp40a1 F AGGTCGTGGGAGTCGGAAAGGA R TGTGGACACTTGCTGGACCATACC gene17675 ggn F GCTGAAGTGCCACCTGAGTCACA R GCCGCTGTTGTATTGCTGCTCTG [62]XM_019257255.1 β-actin F CTGTCCCTGTATGCCTCTGGTC R CTTGATGTCACGCACGATTTCC [63]Open in a new tab 3. Results 3.1. Identification of Testes in Developmental Stage IV We successfully screened L. crocea individuals with testes in developmental stage IV. All germ cell populations from spermatogonia to spermatids were observed in these testes. In addition, many spermatozoa differentiated from spermatids were released into the lobular cavity ([64]Figure 2). Figure 2. [65]Figure 2 [66]Open in a new tab Developmental stage identification of testes by histological analysis. The results show that the testes were in developmental stage IV. The following germ cell populations were observed in the testes: SG, spermatogonia; PS, primary spermatocyte; SS, secondary spermatocyte; ST, spermatid; and SP, spermatozoon. 3.2. Overview of the Transcriptome Profiles An overview of the reads and quality filtering of the two libraries is presented in [67]Table 3. The Illumina HiSeq 2500 platform yielded 26.2 million and 28.78 million clean reads from the TES (testes) library and the PMT (pooled multiple tissues) library, respectively; and 33,390,310 and 38,470,859 reads were mapped to the L. crocea genome, representing 63.72% and 66.82% of the clean reads from the two samples, respectively. The mapped reads represented nearly 70% of the L. crocea genome. The Q30 values (percentage of sequences with a sequencing error rate lower than 0.1%) in the two libraries were 85.81% and 86.27%, respectively, and the GC percentages of the libraries were almost 50%, indicating that the sequencing results were reliable ([68]Table 3). Table 3. Summary of L. crocea transcriptome data. Samples TES PMT Clean reads 26,201,689 28,786,700 Clean bases (Gb) 7.8 8.5 GC percentage 50.73% 50.16% % ≥ Q20 93.38% 93.64% % ≥ Q30 85.81% 86.27% Total mapped 33,390,310 (63.72%) 38,470,859 (66.82%) Uniquely mapped 30,802,638 (58.78%) 35,380,000 (61.45%) [69]Open in a new tab For an overview of the function of all genes from our dataset, the 25,787 L. crocea genes were annotated based on multiple databases using the BLASTx or BLASTn algorithm (E-value ≤10^−5). A total of 25,384 (98.44%) genes were annotated, with 25,378 (98.41%) genes annotated in the Nr database, followed by the Swiss-Prot (17,642; 68.41%), KEGG (15,140; 58.71%), and GO (13,813; 53.57%) databases ([70]Figure 3). To assess their evolutionary conservation, the 25,378 genes mapped to the Nr database were searched against the sequences of other species in the Nr database using the BLASTx algorithm. The results show the highest homology with L. crocea (80.0%) ([71]Figure 4). Figure 3. [72]Figure 3 [73]Open in a new tab Venn diagram showing the annotation by four databases. For homology annotation, non-redundant sequences were searched using public databases, including the Nr, Swiss-Prot, KEGG, and GO databases. Figure 4. [74]Figure 4 [75]Open in a new tab Species that match the annotated sequences of L. crocea. The species distribution is based on the results of the Nr annotation. The percentages of contigs that are homologous with other species are presented. According to the GO annotation result, 13,813 of all genes were mapped to the GO database and classified into three major functional categories (biological process, cellular component, and molecular function) and 59 subcategories. Amongst these genes, 10,237 genes (74.11%) were categorized into the “biological process” functional group; 7396 genes (53.54%) were categorized into the “cellular component” functional group; and 11,697 genes (84.68%) were categorized into the “molecular function” group ([76]Figure 5). Figure 5. [77]Figure 5 [78]Open in a new tab All genes GO cluster distribution. Genes were classified into three main categories: Biological process, cellular component, and molecular function. 3.3. Differentially Expressed Genes (DEGs) of the Two Libraries A total of 21,934 genes did not show significant differences in expression in both groups, while 2194 genes (testes-biased genes) were highly expressed in the TES and 1659 genes (somatic-biased genes) were highly expressed in the PMT ([79]Figure 6a,b). Based on the Q value < 0.01 and |log[2] (fold change)| ≥ 2 threshold, 3853 DEGs were identified between the TES and PMT. Among the 3853 DEGs, 332 were novel genes. Detailed DEGs information is shown in [80]Table S1. The expression of DEGs from the two samples clustered into two distinct groups based on hierarchical clustering ([81]Figure 7). Figure 6. [82]Figure 6 [83]Open in a new tab Comparative results of RNA-Seq and differentially expressed gene distributions between the TES and PMT. (a) Venn diagram showing genes expressed highly in the TES (light purple circle), expressed highly in the PMT (light yellow circle), and expressed without significant difference in both samples (intersection). (b) Volcano scatter plot of differentially expressed genes (TES vs. PMT). Red points represent the testes-biased genes with log[2] (fold change) > 2 and Q value < 0.01, i.e., −log[10] (Q value) ≥ 2.0. Green points represent somatic-biased genes with log[2] (fold change) < −2 and Q value < 0.01, i.e., −log[10] (Q value) ≥ 2. Black points represent genes with no significant differences. Fold change equal to normalized gene expression of the TES/normalized gene expression of the PMT. Figure 7. [84]Figure 7 [85]Open in a new tab Hierarchical clustering of differentially expressed (DE) mRNAs among the TES and PMT. Heatmap of the count data for DE mRNA libraries for the DEGs identified between the TES and PMT. 3.4. DEGs Annotation and Pathway Analysis Functional prediction and classification of the testes-biased gene sequences were achieved by searching against the COG database, and the DEGs were classified into 24 categories ([86]Figure 8). The cluster “general function prediction only” was the major group (242, 28.34%), followed by “signal transduction mechanisms” (106, 12.41%); “replication, recombination, and repair” (105, 12.30%); “transcription” (102, 11.94%); “cytoskeleton” (46, 5.39%); “posttranslational modification, protein turnover, chaperones” (38, 4.45%); and “cell cycle control, cell division, chromosome partitioning” (30, 3.51%). No DEGs were sorted into the “extracellular structures” and “nuclear structure” categories ([87]Figure 8). Figure 8. [88]Figure 8 [89]Open in a new tab Clusters of orthologous groups (COG) classification; 541 testes-biased genes with Nr hits were grouped into 23 COG classifications. According to the GO analysis, 776 testes-biased genes were mapped to the GO databases and classified into three major functional categories (biological process, cellular component, and molecular function) and 38 subcategories ([90]Figure 9 and [91]Table S2). Amongst these, 634 genes (81.70%) were categorized into the “biological process” functional group; 438 genes (56.44%) were categorized into the “cellular component” functional group; and 589 genes (75.90%) were categorized into the “molecular function” group ([92]Table S2). Moreover, among the 38 significantly enriched GO subcategories (Q < 0.05), 19 were associated with “biological process”, 11 were associated with “cellular component”, and 8 were associated with “molecular function”. For biological processes, the most frequently occurring GO term was “cellular process”. For cellular components, the terms “cell” and “organelle” showed the highest frequency. In the category of molecular function, the term “binding” accounted for the largest proportion of annotations, followed by “transporter activity” ([93]Table S2). Some of these enriched GO subcategories are potentially involved in spermatogenesis and testicular development, such as “reproduction”, “reproductive process”, “cellular process”, “biological regulation”, and “developmental process”. Figure 9. [94]Figure 9 [95]Open in a new tab Gene ontology (GO) classifications of testes-biased genes. The X-axis presents the number of genes; the Y-axis presents 38 function subcategories. In addition, 32 third-level GO functional categories closely related to spermatogenesis and testis development were significantly enriched from the present dataset with KS < 0.05 ([96]Table 4). These categories covered multifaceted physiological and biological processes of spermatogenesis and testicular development and offered many candidate genes for further investigation. Table 4. Biological processes involved in spermatogenesis and testicular development. GO.ID The Third-Level Functional Categories All Gene Number DEG Number KS ^a GO:0001539 cilium or flagellum-dependent cell motility 18 12 4.60E-06 GO:0007018 microtubule-based movement 61 31 7.90E-05 GO:0007283 spermatogenesis 20 11 0.00017 GO:0003006 developmental process involved in reproduction 70 23 0.0035 GO:0048515 spermatid differentiation 10 5 0.00452 GO:0007126 meiotic nuclear division 22 10 0.00474 GO:0044702 single organism reproductive process 89 30 0.00857 GO:0071695 anatomical structure maturation 6 1 0.00991 GO:0007286 spermatid development 9 4 0.01166 GO:0046903 secretion 64 8 0.01212 GO:0090068 positive regulation of cell cycle process 12 2 0.01577 GO:0033143 regulation of intracellular steroid hormone receptor signaling pathway 6 1 0.01586 GO:0007281 germ cell development 29 11 0.01627 GO:0035265 organ growth 11 4 0.01688 GO:0000910 cytokinesis 25 7 0.018 GO:0051495 positive regulation of cytoskeleton organization 23 4 0.01847 GO:0048729 tissue morphogenesis 289 50 0.01997 GO:0002009 morphogenesis of an epithelium 236 38 0.02012 GO:0051653 spindle localization 12 1 0.02114 GO:0048608 reproductive structure development 38 9 0.02247 GO:0061458 reproductive system development 38 9 0.02247 GO:0007010 cytoskeleton organization 223 50 0.02305 GO:0045010 actin nucleation 14 4 0.02346 GO:0007548 sex differentiation 33 9 0.02645 GO:0060070 canonical Wnt signaling pathway 77 12 0.03034 GO:0008406 gonad development 28 8 0.03043 GO:0060029 convergent extension involved in organogenesis 12 1 0.03087 GO:0051493 regulation of cytoskeleton organization 63 11 0.03124 GO:0007017 microtubule-based process 170 58 0.03224 GO:0030029 actin filament-based process 132 30 0.0363 GO:0046330 positive regulation of JNK cascade 11 1 0.04252 GO:0030036 actin cytoskeleton organization 130 30 0.04805 [97]Open in a new tab ^aKS < 0.05 indicates statistical significance of enrichment. To further identify molecular interaction networks within cells, the DEGs were mapped to the KEGG Pathway Tools and assigned to a total of 152 pathways ([98]Table S3). The top 50 pathways are presented in [99]Figure 10. Among them, several signaling pathways were identified and documented as essential in gonadal development and maturation, including the following pathways: Neuroactive ligand–receptor interaction pathway (65 genes), MAPK signaling pathway (60), focal adhesion (50), purine metabolism (38), regulation of actin cytoskeleton (38), cell cycle (25), adherens junction (24), oocyte meiosis (20), pyrimidine metabolism (18), pyrimidine metabolism (18), ubiquitin-mediated proteolysis (16), GnRH signaling (13), steroid hormone biosynthesis (12), spliceosome (5), mismatch repair (3), DNA replication (2), and RNA polymerase (2) ([100]Figure 10). Figure 10. [101]Figure 10 [102]Open in a new tab KEGG classifications of DEGs. The X-axis presents number and percentage of annotated DEGs; the Y-axis presents top 50 pathways. Based on the literature review results, the role of the cell cycle signaling pathway in spermatogenesis was highlighted. There were 23 significantly TES-biased genes (cdc14a, cyclin b3, e2f5, cyclin d2, espl, fzr1, cdkn2d, bub1, bub1 beta, cyclin-a1, cyclin-b2-like, wee 2, cdk 2, orc5, rad21, cyclin b2, cdc20, cyclin-a2, cdk 4, stag1, cdc25-1-b, cdca20, and smc1b) and 2 significantly PMT-biased genes (cdk6 and myc-2) found in the cell cycle ([103]Figure 11). Figure 11. [104]Figure 11 [105]Open in a new tab Cell cycle pathway. Green, significantly decreased expression; blue, proteins encoded by both up- and downregulated genes; red, significantly increased expression. 3.5. Validation of RNA-Seq Data by qPCR To verify the reliability of the transcriptome data, we measured the mRNA expression levels of 15 DEGs by qPCR. The results show that the data from the qPCR and the transcriptome sequencing tend to be consistent ([106]Figure 12), and the Pearson’s correlation coefficient of 0.814 between the transcriptome and qPCR gene expression results confirmed the reproducibility and reliability of the Illumina sequencing. Figure 12. [107]Figure 12 [108]Open in a new tab Validation of RNA-Seq data using qPCR. (a) Consistency of log[2]FoldChange (log[2]FC) between RNA-Seq (X-axis) and qPCR assay (Y-axis) is high (R^2 = 0.814) based on the 15 selected genes. (b) Log[2]FC of the selected DEGs by qPCR compared with RNA-Seq data. The log[2]FC values of qPCR are shown as the means ± standard error (SE). 4. Discussion Although L. crocea is an economically important marine fish in China, comprehensive studies on its gonadal development and gametogenesis molecular mechanisms are lacking. Descriptive and quantitative RNA-Seq analyses are crucial for elucidating functional elements of a genome and uncovering the molecular constituents of cells and tissues. RNA-Seq technology has emerged as a powerful tool for studying the molecular mechanisms of certain biological processes in all organisms. To identify the genes and pathways related to spermatogenesis and testicular development of L. crocea, we performed, for the first time, a comparative transcriptome analysis between the testes (stage IV) and other somatic tissues. The testes in stage IV were chosen for the testicular transcriptomic analysis because we focused primarily on genes involved in spermatogenesis, and the process of spermatogenesis in this stage was active [[109]19] and presented typical features, such as the emergence of spermatogonia, primary and secondary spermatocytes, spermatids and spermatozoa, and the release into the lobular cavity of many spermatozoa that had developed from spermatids ([110]Figure 2). A total of 25,787 genes were obtained from the two libraries, and 3853 DEGs were screened using a rigorous set of thresholds. Moreover, we identified numerous important functional genes and signaling pathways involved in testicular development and spermatogenesis in L. crocea by combining the results of the enrichment analysis, annotations, and manual literature searches. In this study, we mainly focused on the most important functional genes and pathways, including the neuroactive ligand–receptor interaction pathway, GnRH and MAPK signaling pathways, cell cycle pathway, and dynein, kinesin, myosin, actin, hsp, sycp2, dmrt1, spata, ddx, tdrd, and piwi genes; and we attempted to combine our results with the references and explain the probable role of the