Full-length transcriptome analysis of Zanthoxylum nitidum (Roxb.) DC.

Plant materials

All analytical methods were carried out in accordance with relevant guidelines and regulations. Fresh roots, leaves, flowers, fruits, old stems, and young branches were all collected from the same Z. nitidum (3-years old) plant in the Wuxu planting base (Nanning, Guangxi Province, China), immediately frozen in liquid nitrogen, and then stored at −80 °C until RNA extraction.

RNA extraction

The total RNA of each sample was isolated using the Trizol RNA extraction kit (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s instructions, then treated with RNase-free DNase I (TianGen, Beijing, China) to remove DNA contaminants.

The integrity and concentration of the RNA were evaluated using the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA) and a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). High-quality RNA of each sample was equally mixed as one pool for full-length transcriptome sequencing.

Library construction and transcriptome sequencing

Full-length cDNA was synthesized from the purified total RNA using the SMARTer PCR cDNA Synthesis Kit (Takara Clontech Biotech, Dalian, China) following the manufacturer’s protocol, and large-scale PCR was conducted to produce more double-stranded cDNA templates. Size selection was then performed to generate SMRTbell™ libraries using a PacBio Template Prep Kit (PacBio, Menlo Park, CA, USA). Subsequently, full-length transcriptome sequencing of Z. nitidum was performed using the Pacific Sequel platform.

SMRT sequencing data processing

Raw reads were processed into circular consensus sequencing (CCS) reads by adapting the PacBio SMRT analysis software v2.3.0 (https://www.pacb.com/products-and-services/analytical-software/smrt-analysis/). Full-length nonchimeric (FLNC) transcripts were determined and generated by searching for both the 5′ and 3′ cDNA primers and the poly A tail signal in CCS.

Consensus isoforms and full-length (FL) consensus sequences were determined using an iterative clustering for error correction (ICE) clustering analysis of FLNC. High-quality FL transcripts (identity >0.99) were acquired by removing redundant sequences using CD-HIT (Li & Godzik, 2006).

Structure analysis and lncRNA prediction

Candidate coding regions of non-redundant transcript sequences were identified by TransDecoder (https://github.com/TransDecoder/TransDecoder/releases).

LncRNAs were identified based on the threshold of transcripts with lengths >200 nt using the predictor of long non-coding RNAs and messenger RNAs based on an improved k-mer scheme tool (PLEK), the coding potential calculator (CPC2), and the coding potential assessment tool (CPAT).

Functional annotation

Non-redundant transcripts were functionally annotated using the following databases: nonredundant protein sequence database (Nr), Swiss-Prot database, TrEMBL database, Kyoto Encyclopedia of Genes and Genomes (KEGG), euKaryotic Ortholog Groups (KOG), Protein family (Pfam), and Gene Ontology (GO).

Identification and characterization of SSRs and transcription factors

MIcroSAtellite (MISA) software (http://pgrc.ipk-gatersleben.de/misa/) was used to identify SSRs within the transcripts, and the characteristics of the repeated motif types were analyzed using the methods described by Feng et al. (2021). The TFs of Z. nitidum were identified using hmmsearch based on the Protein Family (Pfam) search results of the TF family.

SSR validation

A set of 100 primer pairs were randomly selected and synthesized (Table S1). Leaf samples of thirty Z. nitidum trees from 10 different regions (n = 3), consisting of all of four types of trees, were collected from the Guangxi autonomous region and Guangdong Province, China. Detailed information about the samples collected are listed in Table S2. Genomic DNA was extracted using the cetyltrimethylammonium bromide (CTAB) extraction method as reported by Cheng et al. (2009). Polymerase chain reaction (PCR) was initiated with 1 μL template DNA (20 ng), 0.5 μL forward primers, 0.5 μL reverse primers, 5 μL 2× TaqPCR Master Mix and 3 μL sterile distilled water. PCR assay was conducted with the following condition: initial denaturation 95 °C for 5 min, denaturation 95 °C for 30 s, annealing 10 cycles 62–53 °C for 30 s, denaturation 95 °C for 30 s, annealing 30 cycles 52 °C for 30 s, and final extension at 72 °C for 20 min. Additionally, 21 tail sequences (5′-GAAGGTGACCAAGTTCATGCT-3′) of the forward primer of these SSR loci were added for detection by the ABI 3730XL DNA Sequencer. The polymorphism estimation, including polymorphism information content (PIC), observed heterozygosity (Ho), and expected heterozygosity (He) was then performed, followed by a cluster analysis based on the validated SSRs.

Validation of transcript assembly

A total of 10 assembled transcripts were randomly selected to conduct validation using reverse transcription PCR (RT-PCR) amplification. The primers were designed based on the sequence of the de novo assembled transcripts (Table S3) by primer3. PCR amplification included initial denaturation at 95 °C for 5 min, 35 cycles of denaturation at 95 °C for 30 s, annealing at 60 °C for 30 s, extension at 72 °C for 30 s, final extension at 72 °C for 5 min, and storage at 4 °C. Finally, the amplified products were analyzed using agarose gel electrophoresis.

PacBio single molecule long-read sequencing data analysis

The full-length transcriptome of Z. nitidum was obtained using PacBio SMRT sequencing technology. A total of 56.27 Gb of sequencing data were generated. After removing the adapter sequences, approximately 31,517,099 subreads remained with an average read length of 1,785 bp. To provide more accurate and reliable sequences, 456,109 CCSs with a mean length of 2,216 bp (Table 1), and 353,932 FLNC were generated.

A total of 16,555 high-quality FLNCs were obtained after clustering and removal of redundant sequences, and a subsequent analysis revealed 16,163 consensus transcripts for further annotation (Table S4). The length ranged from 147 bp to 7,917 bp, and the N50 and N90 were 1,617 bp and 621 bp, respectively. The mean transcript length was 1,171 bp (Table 2) with 51.49% of the transcripts less than 1,000 bp and 44.24% between 1,000–3,000 bp in length (Fig. 1).

Functional annotation of transcripts

Of the 16,163 transcripts, 14,231 (88.05%) were successfully annotated against the Nr, Swiss-Prot, KEGG, KOG, GO, TrEMBL, and Pfam databases. The annotation rates were 14,199 (87.85%) in TrEMBL, 14,169 (87.66%) in Nr, 11,887 in Swiss-Prot (73.54%), 11,590 in Pfam (71.71%), 9,182 in KOG (56.81%), 7,684 in GO (47.54%), and 6,830 in KEGG (42.26%; Table 3).

The GO analysis revealed that 7,684 unigenes were clustered into 33 GO terms, of which 52.22% unigenes were enriched in binding for molecular function, 47.39% unigenes were enriched in metabolic process, 27.51% unigenes were enriched in membrane parts.

Furthermore, the KEGG enrichment analysis classified 13,829 high-quality transcripts into 34 signaling pathways, which were involved in five categories: cellular processes (1,226, 7.58%), environmental information processing (1,525, 9.44%), genetic information processing (1,736, 10.74%), metabolism (5,048, 31.23%), and organismal systems (2,081, 12.88%; Fig. 2).

Transcript validation

Reverse transcription PCR (RT-PCR) amplification was used to validate the transcriptome assembly. The agarose gel analysis of the target products showed that the corresponding cDNA fragments were approximately the same as the expected size based on the transcript assembly, which demonstrated that the transcripts generated using SMRT technology could be reliably used for further gene identification and functional analysis (Fig. 3).

Identification and characteristic analysis of lncRNAs

LncRNAs are RNA molecules those are more than 200 nucleotides in length which are not translated into protein. By filtering and excluding isoforms with lengths 200 nt, 3,631; 3,681; 1,186; and 4,589 lncRNAs were evaluated based on the CPC, CPAT, PLEK and Pfam databases, respectively. Only 687 of 6,255 total lncRNAs were found in all four computational approaches (Fig. 4).

Transcription factor detection

A total of 3,482 TFs were identified. The most abundant TF families were the bHLH, MYB_related, ERF, and NAC families (Fig. 5).

SSR detection

A total of 22,780 SSRs were identified using the MISA tool, including 13,800 mononucleotides (60.58%), 928 dinucleotides (4.07%). A total of 1,144 trinucleotides (5.02%), 53 tetranucleotides (0.23%), 14 pentanucleotides (0.06%), and 39 hexanucleotides (0.17%; Table 4).

The number of repeat SSR motifs ranged from 5 to 496, with a mean of 147.432. We found that SSRs with 10 motif repeats were the most common and accounted for 7.04% (1,603) of all SSRs, followed by SSRs with 11 repeats (1,028, 4.51%), six repeats (764, 3.35%), and 12 repeats (728, 3.20%); 17,060 SSRs had motif repeat numbers ≥12, accounting for 74.89% of all SSR loci identified (Table S5).

The statistical analysis of all SSR loci showed that the repeat motif types with the highest numbers were: ATC/GAT (64), GAA/TTC (62), and AGC/GCT (51). In Z. nitidum, A/T was the most common mononucleotide repeat motif, accounting for 87.19% (19,861) of all mononucleotide repeats, while C/G represented only 0.79% (179; Fig. 6A). Of the dinucleotide repeats, the AG/CT motif was the most frequent (327, 1.44%), followed by GA/TC (299, 1.31%; Fig. 6B).

SSRs validation

A total of 100 randomly-selected SSR primer pairs were chosen for verification. A total of 36 SSR primer pairs were successfully amplified and showed expected product size in all tested samples; 20 loci showed allelic polymorphism. The observed heterozygosity (Ho) ranged from 0.069 to 0.958, with an average of Ho = 0.478 (Table 5). A total of 126 alleles were obtained from 20 SSRs, and the number of alleles ranged from three to 16 per locus. Dominate alleles were found at the LMZ 51 locus, followed by LMZ73, LMZ96 1903, and LMZ 97. The average value of He was 0.674, but ranged between 0.494 and 0.916. The calculated PIC for each locus ranged from 0.382 to 0.916, with an average of 0.627.

A genetic correlation analysis was performed based on the verified SSR. The results showed three main clusters comprised of four total branches (Fig. 7). Cluster I had two branches, with branch one including nine genotypes belonging to type III and branch II including five genotypes collected from Baise, Guangxi province. Cluster II included 12 genotypes belonging to type I. Cluster III was comprised of three individuals belonging to type IV.