Single molecule, full-length transcript sequencing provides insight into the TPS gene family in Paeonia ostii

Plant materials and sample preparation

Young leaves, stems, roots, flowers, seeds and fruit pods of three different 3-year-old P. ostii from the germplasm repository of the Horticulture and Plant Protection College, Yangzhou University, Jiangsu Province, P.R. China (32°23′31′N,119°24′50′E) were collected and stored in liquid nitrogen for RNA extraction. They had been grown under the same conditions for three years. Seeds, leaves and stems used for qRT-PCR were collected from the same strain of three different P. ostii at the developmental stages of 30, 50 and 70 days after flowering (DAF) in 2019, because the endosperm of seeds begins to develop after 30 DAF and becomes cellular after 70 DAF (Yao, 2020). Leaves, stems, roots, and flowers in the bloom period of P. ostii were collected in April 2019. Fruit pods and seeds were collected at 70 DAF in June 2019. RNA was extracted separately from leaf, stem, root, flower, seed and pod tissues (Figs. 1A–1E) using the CTAB method. To ensure the accuracy of the sequencing data, all RNA samples were checked using a Nanodrop 2000C (Thermo Scientific). The RNA integrity was checked using an Agilent 2100 bioanalyzer (Agilent Technologies), which included RIN, 28S/18S and 5S peaks. Electrophoresis was used to detect whether the RNA samples contained gDNA contamination and assess RNA quality by identifying the ribosomal bands. gDNA were digested by using the Exonuclease I (NEB, M0293S).

PacBio library construction and sequencing

After RNA quality inspection, we first mixed equal quantities of high-quality RNA from different organs of one P. ostii. Then, the mixed samples from three different P. ostii were mixed as one pool (Fig. S1). Ten micrograms of RNA from the sample pool was manipulated as follows: (1) reverse transcribed into cDNA using the SMARTer PCR cDNA Synthesis Kit to prepare full-length, high-quality cDNAs with VN Primer (2 µM Nanopore) and Strand-Switching Primer (10 µM Nanopore) (Table S2); (2) a BluePippin™ Size Selection System (Sage Science, Beverly, MA) was used to screen full-length cDNA fragments and construct cDNA libraries of different sizes (1–2 kb, 2–3 kb, and 3–6 kb); (3) the amplified full-length cDNA was reamplified by PCR and end-repair of the full-length cDNA; (4) then, each SMRT bell library (500 ng size-selected cDNA) was constructed using the Pacific Biosciences DNA Template Prep Kit 2.0; (5) after digestion by Exonuclease I (NEB, M0293), secondary screening was carried out with BluePippin to obtain a sequencing library. After the completion of the library construction, the quality of the library was tested: Qubit 2.0 was use for accurate quantification. The library size was tested using the Agilent 2100, and the library size was as expected before on-board sequencing (the concentration of Qubit >2 ng/µL and the total amount >50 ng/µL). After the library was purified by AMpure PB Beads, full-length transcriptome sequencing was performed using PacBio RSII according to the target data volume.

Full-length corrected transcript collection in SMART sequencing

Full-length transcriptome sequencing based on PacBio SMRT (Sequencing By Synthesis, SMRT) single-molecule real-time sequencing technology does not require the disruption of RNA fragments, and the full-length cDNA had been obtained by reverse transcription using RACE technology. The long read (median 10 kb) of the platform contains a single complete transcript sequence information without post assembly analysis (Gordon et al., 2015; Thomas et al., 2014). The analysis process for acquiring the full-length transcriptome mainly includes three stages (Sharon et al., 2014): (1) Full-length sequence identification. Convert all original sequences to ROI (Reads Of Insert) sequences based on the adaptor. According to whether there are 3′ primers, 5′ primers, or Poly A (optional), they are divided into two categories: full-length and non-full-length sequences. (2) Isoform level clustering to obtain consensus sequences. ROI sequences from the same transcript were clustered using the ICE (Iteratively Clustering and Error Correction) algorithm, and each cluster obtained a consensus sequence. (3) Consistent sequence polishing. The non-full-length sequence was used to perform the polishing on the obtained consensus sequence to obtain a high-quality sequence for subsequent analysis. The resulting transcript sequences can be directly used for subsequent isoform, homologous gene, gene family, SSR, alternative splicing, lncRNA and other analyses.

The PacBio RS II used for SMRT sequencing contained 150,000 Zero-Mode Waveguides (ZMWs) per cell. The reads were sequenced into ZMW wells, and one ZMW contained reads as valid data. The ROI sequences satisfying the conditions of full passes ≥0 and accuracy ≥0.75 were extracted. To obtain the iteratively clustered sequences, the iterative clustering for error correction algorithm was performed by the SMRT Analysis (v2.3.0). Combined with a non-full-length sequence, high-quality transcripts (HQs) with accuracies greater than 99% were obtained. In the process of full-length transcript sequencing, the 3′ end has a poly-A structure, and it can be determined that the 3′ end is relatively intact, and the 5′ end sequence may be degraded, resulting in different copies of the same transcript being assigned to different clusters. CD-HIT (Li & Godzik, 2006) software was used to combine sequences with high similarity and to remove redundant sequences from among the high-quality transcripts.

Functional annotation, ORFs, SSR analysis and lncRNA identification

We used BLAST to combine the obtained nonredundant transcript sequences with the NR (http://www.ncbi.nlm.nih.gov) (April, 2021), SwissProt (http://www.ebi.ac.uk/swissprot) (April, 2021), GO (http://www.geneontology.org) (April, 2021), COG (http://www.ncbi.nlm.nih.gov/COG) (March, 2021), KOG (https://mycocosm.jgi.doe.gov/help/kogbrowser.jsf) (April, 2021), Pfam (http://pfam.xfam.org/) (March, 2021) and KEGG (http://www.genome.jp/kegg) (April, 2021) databases to obtain annotation information of the transcripts (Altschul, 1999).

TransDecoder software 2.0.1 (https://transdecoder.github.io/) is based on the length of the open reading frame (ORF), log-likelihood score and the comparison of amino acid sequences and protein domain sequences in the Pfam database, which can identify reliable coding sequences (CDSs) from the transcript sequence.

Simple repetitive sequences were identified through MISA software. Seven types of SSRs were identified by analyzing the sequences of the transcripts: mononucleotide, dinucleotide, trinucleotide, tetranucleotide, pentanucleotide, hexanucleotide and compound SSRs.

We filtered out the transcripts with coding potential to obtain the predicted lncRNAs. In this study, the most widely used coding potential analysis method was used to predict the lncRNAs among the transcripts, including four methods: CPC analysis (Kong et al., 2007), CNCI analysis, Pfam protein domain analysis, and CPAT analysis.

Quantitative Real-Time PCR anlaysis

Each organ was represented by three biological replicates with three technical replicates. Total RNA from the seeds, leaves and stems from three P. ostii was extracted by a MiniBEST Plant RNA Extraction Kit (TaKaRa, Japan). Then, the cDNA was synthesized by the PrimeScript RT reagent Kit with gDNA Eraser (TaKaRa, Japan) (Zhao et al., 2020). Before qRT-PCR, we used water instead of cDNA and ran the PCR. The agarose gel electrophoresis showed clear objective band, so there was no pollution in the PCR system. A BIO-RAD CFX Connect Optics Module (Bio-Rad, Des Plaines, IL, U.S.A.) was used to analyze the PoTPS expression levels. Expression values were calculated by the 2^−ΔΔCt comparative threshold cycle (C_t) method (Livak & Schmittgen, 2001). SYBR Premix Ex Taq (Perfect Real Time) (TaKaRa, Japan) with 12.5 µL 2 × SYBR Premix Ex Taq, 2 µL cDNA solution, 8.5 µL ddH₂O, and 2 µL solution of mix primers were the system to perform qRT-PCR. Amplification conditions: 95 °C for 30 s, 40 cycles at 95 °C for 5 s, 52 °C for 30 s, and 72 °C for 30 s. The ubiquitin gene (JN699053) was used as an internal reference for this experiment and the expression level of this reference gene was stable in all organs of Paeonia suffruticosa (Wang et al., 2012). All primers were listed in Table S2.

Analysis of the TPS genes family in Paeonia ostii

To classify the PoTPS genes in P. ostii, Cluster X 2.0.12 software (http://www.cluster-x.org/) was applied for multiple sequence alignment by using protein sequences of Arabidopsis. To construct the phylogenetic tree, neighbor-joining (NJ) method was used by MEGA 7.0 software, and bootstrap values were set as 1000 bootstrap replicates (Saitou & Nei, 1987). The conserved motifs of the PoTPS sequences were identified by the MEME program (https://meme-suite.org/meme/), and the parameters were set as a maximum of 10 motifs and an optimum motif width of 6-200 amino acid residues (Zhao et al., 2018). The conserved domains were visualized using the TBtools software.

Paeonia ostii transcriptome sequencing with SMRT

The transcriptome of the pooled samples was sequenced and analyzed with the PacBio Sequel platform to accurately capture the full-length sequences. A total of three libraries divided by cDNA lengths of 1–2 kb, 2–3 kb and 3–6 kb were constructed. Then, 3,568,378 subreads (8.49 Gb) and 343,264 polymerase reads were obtained after filtering. These data have been deposited in the National Center for Biotechnology Information (NCBI) (BioProject accession: PRJNA688625). The ROI sequence was extracted from the original sequence according to the condition that full passes ≥0 and the sequence accuracy ≥0.75. To evaluate the offline date, the number of ROIs, the number of bases of the ROIs and the mean length of the inserted sequence in each library were counted. In total, we obtained 230,736 ROI sequences and the mean read quality of the insert of each size-selected library was 0.92 (Table S3). The ROI read length distribution of each size-selected library had three consistent peaks (Figs. 2A–2C). By filtering out short fragments <300 bp, the sequences containing both the 3′primers and 5′ primers at the same time and having a poly A tail before the 3′ primers were defined as full-length sequences. After further filtering and analysis, 114,215 full-length nonchimeric (FLNC) reads were obtained. The FLNC read length distribution of each size-selected library is shown in Figs. 2D–2F. The trends of the three peaks were consistent, which was in line with the expectations. The percentage of FLNC reads in all ROI sequences was 49.5% (Fig. 3A), accounting for half of the number of all ROI sequences, combined with the ROI and FLNC sequences, which showed that the quality of the transcriptome was better.

By using IsoSeq of SMRT Analysis software to perform cluster analysis on the full-length sequence and combining the non-full-length sequences, 45,006 identical transcripts were obtained (Fig. 4). The consensus isoform read length distribution of each size bin (<1, 1–2, 2–3, 3–6 and >6 kb) is shown in Figs. 3B–3F, and the average consensus isoform read lengths were 912, 1399, 2346, 3665 and 9297 (Table S4). After correcting the consistent sequences in each cluster by quiver program, only 8,239 low-quality isoforms (LQs) and 36,767 high-quality isoforms (HQs) were obtained. Among them, 0–1 kb high-quality transcripts accounted for the largest proportion (91.58%) of consensus isoform reads (Table S4), indicating that the shorter the cDNA length, the more high-quality reads were obtained. Finally, by using CD-HIT to remove the redundant sequences of the high-quality sequences, 30,100 transcript sequences were obtained for subsequent analysis (Fig. 4).

ORF and transcription factor prediction, SSR analysis and lncRNA identification

TransDecoder software was used to predict the coding protein sequences and to identify full-length ORFs. In our study, a total of 28,415 ORFs were obtained, of which 17,904 were complete ORFs. Sequences with start and stop codons were defined as complete coding sequences. The predicted length distribution of the complete ORFs is as follows (Fig. 5A). Transcripts >500 bp in length were selected. By using MISA software to analyze 30,086 transcripts, 9,789 SSR- containing sequences from among the 13,632 evaluated SSRs were identified. Large numbers of SSR containing sequences were mono-, di-, or trinucleotide repeats (Table S5). Currently, long noncoding RNAs (lncRNAs) play an important role in many life activities, such as dose compensation effects, epigenetic regulation, cell cycle regulation and cell differentiation regulation, and they have become a hotspot of genetic research. Therefore, it is very meaningful to analyze lncRNAs through transcriptome analysis. Four methods were used to make lncRNA predictions for the transcripts CPC, CNCI, Pfam and CPAT. To visually display the analysis results, the noncoding transcripts identified by the above four analysis methods were taken into intersection for subsequent lncRNA analysis (Fig. 5B). A total of 827 lncRNAs were predicted by filtering out the transcripts <300 bp in length.

Functional annotation of the transcripts

We have submitted the final polished consensus mRNA sequences to NCBI (accession number GJFW00000000). Blast software was used to compare the nonredundant transcripts with the NR, SwissProt, GO, COG, KOG, Pfam, and KEGG databases. A total of 28,850 transcript annotation information points were obtained (Table S6). We searched for homologous species through sequence alignment. The 28,790 transcripts in Nr were aligned as shown below, showing that the largest number of transcripts were distributed in Vitis vinifera (37.64%) (Fig. 6A). The GO annotation system includes three main branches, namely biological process, molecular function and cellular component. Among them, biological processes involved metabolism, cellular and single-organism processes, responses to stimulus and biological regulation. Cell components mainly involved cells, cell parts, organelles, membranes and other parts. Molecular functions involve catalysis, binding, transport, transcription activities and other activities. The statistical results of 20,052 transcripts of GO are shown in Fig. 6B. The COG (Cluster of Orthologous Groups of proteins) database was constructed based on the phylogenetic relationships of the bacteria, algae and eukaryotes. The COG database can be used to classify gene products as orthologues. The COG classification statistics of 12,892 transcripts are shown in Fig. 6C.

Identification of the TPS gene family

Ten TPS family members were selected and identified from the transcriptome sequencing with the SMRT database for P. ostii, and the accession numbers for the 10 TPS genes are MW700299–MW700308. According to the percentage of similarity between these TPS family members and AtTPS1-11 in Arabidopsis, we named these genes PoTPS1, PoTPS3-PoTPS11 (Table S6). To clarify the related protein information of the TPS family in P. ostii, ProtParam was used to analyze the physical and chemical properties and secondary structure elements of the TPS family members in this study. The number of amino acids was approximately 429-959, the isoelectric point ranged between 5.61−6.37 and the protein molecular weight was between 48.0–108.5 kD. The secondary structure is mainly composed of a helix, extended strand and random coil. Specifically, the random coil structure of 10 TPS members was above 29%, which is beneficial to maintaining the stability of TPS proteins (Table 1). Using Pfam and SMART to analyze protein domains, it was found that these TPS proteins all have conserved domains, which contain two typical characteristics of TPS (Pfam: Glyco-transf-20) and TPP (Pfam: Trehalose-PPase) structural domains. In particular, PoTPS1 contains only one TPS domain, and the remaining nine TPS family members all contain one TPS domain and one TPP domain. In addition, it was also found that the distribution of the domains of the TPS family members is roughly the same as that of Arabidopsis, rice and soybean, showing that the TPS and TPP domains of this family are positionally conserved among different species. However, whether these family members are functionally similar is not yet known. To analyze the genetic relationship and evolution characteristics of the TPS members in P. ostii, a phylogenetic tree was constructed by using MEGA 7.0 software, finding that the PoTPS1 protein has a high similarity to AtTPS1 and that the PoTPS7 protein has a high similarity to AtTPS7 (Fig. 7). In addition, the motif is a small conserved sequence fragment. In biology, a mathematical statistical model based on data can help predict the reliability of phylogenetic analysis, so MEME was used to analyze the motifs of all TPS proteins and it was found that only PoTPS5, PoTPS7 and PoTPS8 had ten conserved motifs, while the other proteins had 5-9 motifs each (Fig. 7).

Gene expression analysis

We analyzed the expression levels of 10 TPS family members on the stems, leaves and seeds at 30, 50, and 70 DAF to explore whether the TPS gene expression in different tissues and different periods follows certain expression patterns and whether TPS genes were specifically expressed in different tissues. TBtools software was used to analyze the obtained PCR results and the red color from light to dark represents the expression level from low to high. The results showed that in seeds the relative expression levels of PoTPS5 and PoTPS6 were relatively high at 30 DAF, and the relative expression levels of PoTPS5 and PoTPS11 were high at 50 DAF. The relative expression levels of PoTPS1 and PoTPS5 were relatively high at 70 DAF. Obviously, the relative expression level of PoTPS1 in seeds tended to increase with the development of the seeds. In leaves, the expression of PoTPS11 was higher at 30 DAF, the expression of PoTPS5 and PoTPS8 was higher at 50 DAF, and the expression of PoTPS5 was highest at 70 DAF. In stems, the relative expression level of PoTPS5 was high regardless of the period (Fig. 8). In summary, the relative expression of PoTPS5 was always been at a high level in different tissues at different periods.