H. manillensis larval (10 days after birth), young (3 months after birth), and adult worms (2 years after birth) were collected from Qinzhou Zhanhong Chinese Medicinal Materials company. The leeches were immediately executed using scissors, and 24 samples were collected from the anterior and posterior segments of larval and young worms, while the adult samples were collected from oral suckers, gonads, caudal suckers, and the remaining worm body. Moreover, each time, five H. manillensis were collected and preserved in water bottles filled with pure water for subsequent histological examination.

The histomorphological characteristics of H. manillensis were observed by conventional paraffin embedding, sectioning, and hematoxylin and eosin staining. A microscope (EVOS FL Auto) was used to examine leech gonads and muscle histomorphology, and images were captured using an EVOS Auto microphotographic system.

Total RNA was extracted from different parts of leeches at different developmental stages using TRIzol reagent (Invitrogen Corp., Carlsbad, CA, United States), and RNA purifications were performed using the RNeasy Mini Kit (Qiagen, Chatsworth, CA, United States). The sequencing libraries were constructed using the NEBNext Ultra RNA Library Prep Kit for Illumina (NEB, United States) following the manufacturer’s recommendations, and the libraries were sequenced on an Illumina HiSeq 4000 platform to generate 150 bp paired-end reads. Each sample quality was controlled with Fastp (v0.20.1) (Chen et al., 2018). Furthermore, we used HISAT2 (v2.2.1) to map the reads from each sample to the reference genome, and Samtools was used to sort and convert the SAM files to BAM; we use the “merge” parameter of the Samtools tool to combine the data from different parts together (Kim et al., 2015). Moreover, StringTie (v2.1.4) was utilized to assemble and merge the transcripts of each sample (Pertea et al., 2015), and those merged transcripts were then subjected to Gffcompare (v0.11.2) to compare with the reference annotation file in GTF, and transcript abundances were estimated using StringTie with options “-e -B -p 20” (Pertea and Pertea, 2020). The results of abundances were stored ended with “balltown,” and prepDE.py was used to compare the folders.

The gene expression was predicted in fragments per kilobase of transcript per million fragments mapped (FPKM), and differentially expressed genes (DEGs) among three stages were detected using the DESeq2 (v1.32.1) package (Love et al., 2014). Genes with FDR (false discovery rate) values < 0.05 and |log²FC| ≥ 1 were considered as DEGs between two stages. Furthermore, the simultaneous mapping of differentially expressed genes was carried out using the ggplot2 (v3.3.5) package (Wickham, 2016).

The time course sequencing data analysis (Tcseq) (v1.18.0) software was used to analyze the DEGs at three time points and clusters and graphically illustrated the different gene expression patterns at different developmental stages (Wu and Gu, 2021). The genes with FDR < 0.05 and |log²FC| ≥ 1 were considered DEGs between the two stages.

The following highly stringent criterion was used to screen potential lncRNAs: 1) transcripts with >200 nt length were selected; 2) the transcripts were annotated as “i”, “u”, and “x” according to the cuffcompare classes; 3) the coding potential calculator (CPC2) score was “noncoding”; 4) transcripts with an open-reading frame (ORF) of more than 50 amino acids were removed; and 5) transcripts with known protein-coding domains were excised by aligning to the Swiss-Prot database. Furthermore, as a cis-regulator, lncRNA has the potential to influence the expression of adjacent genes; thus, a coexpression network of the candidate lncRNAs and their upstream or downstream 10 kb mRNAs was constructed. The connectivity and enrichment analysis was performed with the Position Frequency Matrix (Azlan et al., 2019).

Gene Ontology (GO, http://www.geneontology.org/) and KEGG (Kyoto Encyclopedia of Genes and Genomes) analyses of the differentially expressed were performed using clusterProfiler (statistical analysis and visualization of functional profiles for genes and gene clusters). The value of FDR < 0.05 was considered significantly enriched.

The redundant samples from transcriptome sequencing were utilized for qRT-PCR analysis. The first-strand cDNA was synthesized with 1 μg total RNA using a HiScript^® III-RT SuperMix for qPCR (+gDNA wiper) kit (Vazyme, R323-01). Primers of the chosen genes were designed by using oligo 7 software (Supplementary Table S3) (Rychlik, 2007). The expression levels of selected genes were quantified using Universal SYBR qPCR Master Mix (Vazyme, Q711-03), and data were acquired using an RT-PCR instrument (ABI 7500), according to the manufacturer’s instructions. For each detection cycle, three biological and three technical replicates were used. The relative gene expression was calculated using the 2^−△△Ct method using RO21 as the reference gene to normalize the relative mRNA expression levels (Schmittgen and Livak, 2008).

All the obtained data were statistically analyzed by using Student’s t-test and analysis of variance (ANOVA) with DUNCAN’s Multiple Range Test (DMR) in SPSS 17.0 software (IBM Corp., Armonk, NY, United States) (SPSSAU, 2017). The p values < 0.05 were considered statistically significant.

First, the collected H. manillensis worms at larval (L), young (Y), and adult (A) stages were fixed with formaldehyde and then projected as shown in Figure 1A. As presented in Figure 1B, a highly significant difference (p < 0.0001) was observed in the bodyweight of H. manillensis at different growth stages, with the larval stage weighing 0.081 ± 0.031 g, young 1.922 ± 0.563 g, and adults 8.329 ± 2.100 g. Furthermore, we also compared the transverse paraffin sections (Figure 1C) of H. manillensis (larval, young, and adult), and only the adult H. manillensis had a significant muscle bundle tissue, while it was not spotted in the larval and young stages.

We have constructed a total of 24 cDNA libraries for larva, young, and adult H. manillensis worms (HMLB1, HMLB2, HMLB3, HMLF1, HMLF2, HMLF3, HMYB1, HMYB2, HMYB3, HMYF1, HMYF2, HMYF3, HMO1, HMO2, HMO3, HMG1, HMG2, HMG3, HMP1, HMP2, HMP3, HMB1, HMB2, and HMB3) (Supplementary Table S1). A total of 54,639,118 sequences were obtained, with a GC content ranging between 43% and 45% (Supplementary Table S2). Furthermore, the sequenced raw data from different growth periods were combined to obtain three samples for each stage (L1, L2, and L3; Y1, Y2, and Y3; and A1, A2, and A3).

Using RNA-seq, about 18,106 mRNAs transcripts were detected, and their expression is presented in FPKM values and the FPKM distribution of mRNAs is shown in Figure 2A, whereas the expression of different samples is presented in the Box chart (Figure 2B). All the nine samples were divided into three parts by PCA analysis, which demonstrated satisfactory repeatability (Figure 2C) and then, a relationship cluster diagram was constructed to represent the relationship between samples (Figure 2D) intuitively. Sequences exhibited a reliable clustering effect, ensuring the veracity of the subsequent analysis.

We used an FDR value < 0.05 and |log²FC| ≥ 1 as the criteria to screen the DEGs by comparing pairwise differences between larva, young, and adult H. manillensis. We obtained 1,700 DEGs after comparing the L vs. Y, 320 from Y vs. A, and 2,799 in L vs. A (Figures 3A,B). We presented VENN plots of these DEGs from L, Y, and A growth stages (Figure 3C), and a total of 47 intersections of DEGs were assessed with 3,356 DEGs, and their expression patterns were illustrated in a heat map using clustering analysis (Figure 3D). We found that both L vs. Y and L vs. A had high differentially expressed genes during the growth of H. manillensis, while some of the differentially expressed genes were significantly reduced during Y vs. A. Therefore, the growth and development of H. manillensis gradually stabilized when they reached the young stage.

Furthermore, we examined the DEG expression trends and screened the DEGs (gene expression profile indicated by FPKM) using FDR < 0.05, and |log²FC| ≥ 1 to show the development of H. manillensis at various growth phases. We restricted our study to nine clusters of which five had overall upregulated trends and four downregulated trends, and the majority of the genes were enriched in the overall upregulated trend, especially the Cluster 3 which enriched a total of 786 DEGs (Figure 4). These findings can effectively reveal the gene expression patterns of H. manillensis in different growth phases.

In addition, we performed two-by-two comparisons between larval, young, and adult stages of H. manillensis at FDR < 0.05 and |log²FC| ≥ 1, and the genes having the value of log²FC ≥ 1 were considered upregulated genes and downregulated with log²FC ≤ −1 and then performed the GO and KEGG enrichment analyses. Using the GO enrichment analysis, we obtained significant enrichment of differentially upregulated and downregulated genes in three components: biological processes, cellular components, and molecular functions of the larval and adult stages (Figures 5A,B). These upregulated and downregulated genes were found to be enriched in biological processes such as cilium organization, cilium movement, cell cycle G2/M phase transition; in cellular components such as cilium, myofilament, microtubule-associated complex; and molecular functions such as motor activity, oxidoreductase activity, and disulfide oxidoreductase activity. Moreover, the KEGG enrichment analysis of the differentially upregulated and downregulated genes for larval and adult stages is shown in Figures 5C,D, where upregulated genes were mainly enriched in cytoskeleton proteins, dilated cardiomyopathy, adrenergic signaling in cardiomyocytes, and glycolysis/gluconeogenesis, while downregulated genes were mainly enriched in the renin–angiotensin system, maturity onset diabetes of the young, mTOR signaling pathway, fatty acid degradation, cytochrome P450, and other pathways. Also, the GO and KEGG profiles of H. manillensis L vs. Y and Y vs. A differential expressed genes are shown in Supplementary Figures S1, S2.

Likewise, among the biological components examined by GO enrichment, the expression of DEGs in the myofibrillar and contractile fiber pathways exhibited a progressive upregulation during the L-Y-A period (Figure 6A). The gene family expression related to the disulfide oxidoreductase activity pathway in molecular functions tends to be gradually upregulated (Figure 6B), and the genes of the PAI gene family were the majorly detected genes. There is also a progressive downregulation of genes that are enriched in response to nutrient levels, extracellular stimuli, and in cellular lipid metabolism (Figures 6C, D). The downregulated genes were also linked with development-related pathways such as gland development, respiratory system development, and compound eye development (Figure 6E). The KEGG enrichment analysis findings at the L-A stage revealed that the significantly downregulated genes were associated with the cytochrome P450 family pathway, which tended to be progressively lower in the H. manillensis’s three growth stages (Figure 6F), notably for the CP2U1 gene. Furthermore, we observed a significantly higher expression of the HIRM2 and KREM1 genes throughout the adult stage and an overall lower expression of the HIRM1 gene family across the three developmental periods (Figures 6G,H). The GPRL1 gene family, which is associated with spermatogenesis, was also found to be highly expressed during the L-Y period, while decreased expression was seen during the Y-A period (Figure 6I).

The DElncRNAs were screened by comparing pairwise differences from three different growth periods of H. manillensis using FDR < 0.05 and |log²FC| ≥ 1 as a criterion, and a total of 107 DElncRNAs were detected of 19 DElncRNAs from L vs. Y, 21 from Y vs. A and 67 from L vs. A (Figure 7A). VENN plot presentation of DElncRNAs for these different growth periods is shown in Figure 7B. We obtained a total of 78 DElncRNAs and were shown in a heat map by clustering analysis (Figure 7C). Moreover, GO enrichment analysis showed that these DElncRNA target genes were mainly enriched in genitalia development, actin cytoskeleton, and actomyosin (Figure 7D) and might target mRNA to affect mitochondrial biogenesis, cyanoamino acid metabolism, and Parkinson’s disease (Figure 7E). The 491 lncRNAs were predicted by cis prediction, and there were 14 lncRAs that could target the eight mRNAs (Figure 7F).

The differentially expressed genes (ERO1A, HAND2, HS12A, and GPRL1) were selected for Q-PCR quantification in order to confirm their differential expression in different growth periods of H. manillensis obtained by RNA-seq. The expression of these genes was confirmed by Q-PCR (Figures 8A,B).