Single-molecule long-read sequencing of the full-length transcriptome of Rhododendron lapponicum L.

Rhododendron lapponicum L. is a familiar ornamental plant worldwide with important ornamental and economic value. However, a full-length R. lapponicum transcriptome is still lacking. In the present study, we used the Pacific Biosciences single-molecule real-time sequencing technology to generate the R. lapponicum transcriptome. A total of 346,270 full-length non-chimeric reads were generated, from which we obtained 75,002 high-quality full-length transcripts. We identified 55,255 complete open reading frames, 7,140 alternative splicing events and 2,011 long non-coding RNAs. In gene annotation analyses, 71,155, 33,653, 30,359 and 31,749 transcripts were assigned to the Nr, GO, COG and KEGG databases, respectively. Additionally, 3,150 transcription factors were detected. KEGG pathway analysis showed that 96 transcripts were identified coding for the enzymes associated with anthocyanin synthesis. Furthermore, we identified 64,327 simple sequence repeats from 45,319 sequences, and 150 pairs of primers were randomly selected to develop SSR markers. This study provides a large number of full-length transcripts, which will facilitate the further study of the genetics of R. lapponicum.

RnA extraction. Total RNA was extracted using TRIzol LS reagent (Invitrogen, USA) following the manufacturer's instructions. RNA degradation and contamination were monitored using 1% agarose gels. The purity, concentration and absorption peak of RNA were measured using a NanoDrop 2000 spectrophotometer (Thermo Scientific, USA). RNA quality was determined with the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent, USA). The total RNA samples from four tissues were mixed together for the following experiments.
Library construction and SMRt sequencing. To construct a full-length transcript sequencing library, 5 μg of mixed total RNA was reverse transcribed into cDNA using the Clontech SMARTer cDNA Synthesis Kit (Takara Clontech Biotech, Dalian, China) following the manufacturer's protocols. Size fractionation and selection (1-2 kb, 2-3 kb, and 3-6 kb) were performed using the BluePippin Size Selection System (Sage Science, USA). Three SMRT sequencing libraries containing fragments of 1-2, 2-3, and 3-6 kb in length were constructed using the Pacific Biosciences DNA Template Prep Kit 2.0. Finally, 1 (1-2 kb), 1 (2-3 kb) and 1 (3-6 kb) SMRT cells were sequenced on the Pacific Bioscience RS II platform.
Quality filtering and error correction. Raw reads were processed by removing polymerase reads (<50 bp in length). The obtained clean reads were processed into error-corrected ROIs with the following parameters: full passes ≥0 and predicted consensus accuracy >0.75. By identifying the 5′ and 3′ adapters and the poly (A) tail, full-length and non-full-length reads were determined from the ROIs. A full-length read containing both the 5′ and 3′ primer sequences and a poly (A) tail was considered to be a full-length transcript. ROIs with all three elements that did not contain any additional copies of the adapter sequence within the DNA fragment were referred to as full-length non-chimeric (FLNC) reads. Corrected FLNC reads were clustered into transcripts using the ICE algorithm in the PacBio SMRT Analysis (v2.3.0) software. Full-length transcripts with a post-correction accuracy >99% were used for further analysis. prediction of oRfs, lncRnAs, tfs and AS events. To predict ORFs in transcripts, the TransDecoder v2.0.1 tool was used to find potential coding sequences. Based on the obtained transcripts with redundancy removed, we predicted AS events with the software AStalavista 19 . TFs were predicted from the putative protein sequences using the Plant Transcription Factor Database v4.0 tool 20 . We identified unique transcripts without protein-coding potential as candidate lncRNAs with four analysis tools: the coding-non-coding index (CNCI) 21 , the coding potential assessment tool (CPAT) 22 , the coding potential calculator (CPC) 23 , and Pfam protein structure domain analysis 24 . functional annotation. All transcript sequences were analysed for homology via searches against the non-redundant nucleotide database (Nr) 25 , Swiss-Prot protein 26 , protein family (pfam) 27 , evolutionary genealogy of genes: non-supervised orthologous groups (eggNOG) 28 , clusters of orthologous groups of proteins (COG) 29 , eukaryotic ortholog groups (KOG) 30 , gene ontology (GO) 31 , kyoto encyclopedia of genes and genomes (KEGG) 32 databases with BLAST alignment (E-value ≤ 10 −5 ). qRt-pcR analysis. Samples of flowers at four flower developmental stages were collected, and RNA was isolated from them using TRizol regent according to the manufacturer's instructions. The cDNA was synthesized using AMV reverse transcriptase XL for RT-PCR according to the manufacturer's instructions 33 . The qRT-PCR was performed under the following conditions: 95 °C for 2 min, followed by 40 cycles of 5 s at 95 °C, 30 s at 55-60 °C, and a final melting curve step. Three biological replicates were performed in a Roche 480 LightCycler. Threshold values (CT) were used to quantify relative gene expression using the comparative 2 −ΔΔCt method 34 . The information of primer used for qRT-PCR analyses is shown in Table S7.
Development of SSR markers. For SSRs analysis, transcripts longer than 500 bp were selected and MISA software was used. The parameters were set for identifcation of mono-, di-, tri-, tetra-, penta-, and hexa-nucleotide motifs with a minimum of ten, six, five, five, five, and five, respectively. SSR primers were designed by Batch Primer 3 tool 6 . PCR amplifications were performed using the DNA template extracted from R. lapponicum. The PCR products were separated in 8% polyacrylamide gels.

Results
SMRt sequencing data. To obtain a representative full-length transcriptome for R. lapponicum, the total RNA from four tissues (root, stem, leaf and flower) were used for the library construction for SMRT sequencing. In this study, we obtained 957,032 polymerase reads from the sequenced library, and a dataset with 15.37 Gb of clean reads was generated (Supplementary Table S1). A total of 658,338 reads of inserts (ROIs) were generated open reading frame and AS event prediction. Using the software TransDecoder, 72,606 ORFs were predicted. A total of 55,255 complete ORFs were identified, and the length distribution of the complete ORFs was analysed (Fig. 2). Among all transcripts obtained by SMRT sequencing, 7,140 AS events were detected (Supplementary Table S2). Due to the lack of an available R. lapponicum reference genome, further characterization of the types of AS events would be warranted in future studies.

GO classification.
To classify the gene functions of the transcripts, GO analysis was performed. GO analysis showed the enrichment of 33,653 transcripts categorized into 51 functional groups, which could be divided into three major categories: biological process, cellular component and molecular function (Fig. 6). In the biological process group, catalytic activity and binding were the main categories. In the cellular component group, cell part, cell and organelle were the most frequent categories. In the molecular function group, the genes were involved in catalytic activity, binding and other categories.
COG classification. To further study the functional annotation and classification of the R. lapponicum transcripts, all transcripts were subjected to a search against the Clusters of COG database. COG analysis showed that 30,359 transcripts were assigned to 24 categories (Fig. 7). The largest group was general function prediction only (8,967,19.69%), followed by transcription (4,862, 10.68%) and then replication, recombination and repair (4,722, 10.37%). The percentages of six groups were less than 1.00%, including RNA processing and modification, nuclear structure, and cell motility.

Discussion
SMRT sequencing technology is an efficient and reliable approach for obtaining the full-length transcripts of certain species 35 . Recently, long-read SMRT sequencing has been the most reliable and efficient strategy for whole-transcriptome profiling studies, especially for nonmodel plant species without reference genome sequences. In this study, SMRT sequencing technology was applied to investigate the R. lapponicum transcriptome using the PacBio RS II platform. A total of 15.37 Gb of sequencing data were generated, including 658,338 ROIs and 346,270 FLNC reads. The percentage of FLNC reads in all ROIs was 52.59%, and this result was similar findings obtained in alfalfa 36 Table 5. The distribution of SSRs based on the number of repeat units. www.nature.com/scientificreports www.nature.com/scientificreports/ full-length transcripts were obtained. SMRT sequencing can capture the very long nucleotide sequences, where one read usually represents a full-length transcript 13 . The length of the transcripts obtained by SMRT sequencing technology is longer than that of transcripts obtained by next-generation high-throughput sequencing technology. In this study, the average length of the R. lapponicum transcripts was 2,509 bp, which was longer than those obtained in seashore paspalum (970 bp) 38 , sweet potato (581 bp) 39 , and sesame (629 bp) 40 by Illumina sequencing technology. Furthermore, we found that 58.66% of all transcripts were longer than 2,000 bp in this study, and much higher than that in Rhododendron molle (7.23%) 41 and Neottopteris nidus (13.63%) 6 using Illumina sequencing technology. These results indicated that PacBio SMRT sequencing technology is an efficient approach to capture the transcript sequences, especially for long transcript sequences.
Alternative splicing is a major cellular mechanism generating transcriptome diversity and proteome complexity in plants 42 . In this study, 7,140 AS events were detected from the R. lapponicum transcripts. In addition, 3,150 TFs that are key components involved in the transcriptional regulatory system were identified. LncRNAs are a novel class of non-coding transcripts with lengths greater than 200 nucleotides that play important roles in many biological processes 43 . LncRNAs are largely involved in regulating plant development and growth, secondary metabolism, and the plant stress response 44 . Recently, an increasing number of studies have focused on the functions of lncRNAs in plants such as in red pineapple 45 and hot pepper 46 . However, no lncRNAs from Rhododendron have been reported. In this study, we identified 2,011 lncRNAs using four methods, and these lncRNAs will be useful for further research in R. lapponicum.
A total of 71,386 transcripts were annotated by sequence alignment in eight databases, suggesting that this study generated a very large number of R. lapponicum genes. The percentage of annotated transcripts was 95.18%, which was consistent with that in alfalfa 36 and shrimp 47 . The remaining 3,616 transcripts presented no BLAST matches and might represent R. lapponicum-specific genes or unknown genes in R. lapponicum. The systematic classification of proteins in the transcriptome is crucial for maximizing the utilization of transcripts for functional and evolutionary studies. The results of GO and COG classification suggested that a large number of transcripts were involved in transcription, replication, recombination and repair, and catalytic activity. There were 31,749 transcripts assigned to specific pathways, such as metabolism, genetic information processing, cellular processes, environmental information processing, and organismal systems pathways. The results of GO, COG and KEGG classification showed that a large number of transcripts had diverse molecular functions and were involved in many biological pathways. Therefore, our data provided abundant genetic information on future molecular survey on the growth and development of R. lapponicum.
Flower colour is one of the most important ornamental characteristics of rhododendrons. The biosynthesis of anthocyanin is critical for a wide range of flower colours. Previous studies have shown that C4H, CHS, F3H, F3'H, F3ʹ5ʹH, DFR and ANS are the key enzymes involved in the biosynthesis of anthocyanin for the determination of different flower colours in plants 48 . In the present study, a total of 96 transcripts were identified coding for the enzymes associated with anthocyanin synthesis. Gene expression analysis by qRT-PCR showed that the expression levels of C4H, F3H, F3 5 H, DFR and ANS genes were low at the early flowering developmental stage and increased as the flowers developed. The increases in the expression of these genes were consistent with the changes in anthocyanin content in the flower petals of R. lapponicum during flower development. In addition, transcription factors such as those of the MYB, bHLH and WD40 families play a key role by regulating the expression of genes in anthocyanin biosynthesis 49,50 . According to the functional annotation results, 3,150 putative transcription factor genes belonging to 64 TF families were identified. Among these genes, 83 and 182 transcripts belonged to the MYB and bHLH families, respectively (Fig. 4). In conclusion, the identification of key enzymes and related regulatory TF genes involved in anthocyanin biosynthesis and metabolic pathways may contribute to the understanding of colour-regulating mechanisms in rhododendrons.
The rapid development of transcriptome sequencing technology has enabled the massive development of SSR markers 51,52 . In total, 64,327 SSRs were identified from 45,319 SSR-containing sequences, and the average frequency of SSRs was one in 2.67 kb. Among the six types of repeat motifs, dinucleotide repeats were the most abundant. In the present study, the most frequent mono-, di-, and tri-nucleotide motifs were A/T, AG/CT and GAA/TTC, respectively, which was consistent with the results of studies in non-heading Chinese cabbage 53 , rubber tree 54 and radish 55 . CT/AG/GA/TC were the most abundant motifs, accounting for 92.11% of the total dinucleotide repeats. CT repeats are typically found in transcribed regions that may be involved in antisense transcription and play a role in gene regulation 56,57 . Furthermore, the most abundant mononucleotide motif was A/T, which is thought to be frequent in the genomic sequences of plants 58 . SSR abundance varies among different plant species in different studies. Repeat units of 10, 6, 11, 7, and 5 in SSR sequences accounted for 51.98% of the total SSRs. A total of 150 pairs of PCR primers were designed, and 127 primer pairs successfully amplified PCR products. The failure of 23 primer pairs to achieve amplification might have resulted from the targeting of amplicons with large introns, primers positioned across splice sites or chimeric primers. These results suggested that the development of SSR markers based on R. lapponicum transcripts obtained from PacBio SMRT sequencing is an effective and feasible approach. The newly developed SSR markers from our study will provide a valuable genetic tool that be used in studies on genetic diversity, comparative genomics, gene mapping, and population genetics and other types of genetic studies in rhododendron.
In conclusion, we analysed the full-length transcriptome of R. lapponicum by using the PacBio SMRT sequencing technology. This study represents the first the third-generation long-read transcriptome sequencing of R. lapponicum. Based on the obtained transcriptome data, 7,140 AS events, 2,011 lncRNAs, 55,255 complete ORFs and 3,150 TF members were identified. A total of 96 transcripts were identified coding for the enzymes associated with anthocyanin synthesis. In addition, 64,327 SSRs were detected, and 150 primer pairs were randomly selected to develop SSR markers. The obtained transcriptome data may facilitate further genetic studies on R. lapponicum.