Transcriptomic profiles of 33 opium poppy samples in different tissues, growth phases, and cultivars

Opium poppy is one of the most important medicinal plants and remains the only commercial resource of morphinan-based painkillers. However, little is known about the regulatory mechanisms involved in benzylisoquinoline alkaloids (BIAs) biosynthesis in opium poppy. Herein, the full-length transcriptome dataset of opium poppy was constructed for the first time in accompanied with the 33 samples of Illumina transcriptome data from different tissues, growth phases and cultivars. The long-read sequencing produced 902,140 raw reads with 55,114 high-quality transcripts, and short-read sequencing produced 1,923,679,864 clean reads with an average Q30 rate of 93%. The high-quality transcripts were subsequently quantified using the short reads, and the expression of each unigene among different samples was calculated as reads per kilobase per million mapped reads (RPKM). These data provide a foundation for opium poppy transcriptomic analysis, which may aid in capturing splice variants and some non-coding RNAs involved in the regulation of BIAs biosynthesis. It can also be used for genome assembly and annotation which will favor in new transcript identification.


Background & Summary
Opium poppy (Papaver somniferum) is one of the most important medicinal plants in the world and remains the only commercial resource of morphinan-based painkillers 1 . Its main active ingredient, BIAs, also displays potential pharmacological activity in relieving cough, muscle relaxation, anticancer, and so on 2,3 . Although approximately 100,000 hectares (ha) of opium poppy are cultivated annually worldwide, this insufficient to meet the demand for managing moderate or severe pain 1,4 . Engineered microbes could be used to produce BIAs such as opiates and noscapine [5][6][7] , however, major hurdles, such as low yield and unclear biosynthetic pathways, make it different to scale-up this method of production for most BIAs. Hence, how to guarantee the source of BIAs to meet the medical applications, such as in pain relief or palliative care, has been recognized as a major issue that needs to be resolved.
Investigation into the BIA biosynthesis mechanism began in the 1960s with radiotracer technology 3 . The emergence of transcriptomics, proteomics, and metabolomics, coupled with recent genome analysis tools, accelerated the discovery of new BIA biosynthetic genes that could facilitate metabolic engineering reconstitution of commercial source of valuable BIAs in microbes 1,5,[7][8][9][10][11] . However, some key steps in the BIA biosynthesis pathway are yet to be identified, and there has been a lack of investigation into the molecular mechanisms of gene regulation in the pathway 3 . In addition, there are few reports major in the processes of BIAs biosynthesis, such as compound dynamic accumulation, tissue-specific distribution, enzyme interactions and metabolism, compartmentalization and transport 3 . Therefore, investigation into the regulation mechanism and understanding the compound accumulation process is a key way to improve the yield of BIAs.
Of all the strategies in metabolic regulation, over-expression and silencing of genes involved in the metabolism of target compounds are the most widely used method. However, non-coding RNAs, alternative splicing/ translation/polyadenylation (AS/AT/APA), formation of heterodimers, and gene fusion have also been shown to increase the flexibility of the transcriptome, functional complexity of plants, and the trend of metabolic flow, tissue-specific accumulation or product yield 8,9,12,13 . However, little research has been conducted on the regulatory mechanisms in opium poppy despite the recently release of its genome 1 . The reason for this oversight may be the lack of the transcriptional information on opium poppy in different growth periods/status or tissues. In addition, the short sequences created by third-generation sequencing and the gene information in DNA standard could not capture the AS/AT/APA of transcripts 1,3 .
Third-generation single-molecule real-time (SMRT, Pacific Biosciences) sequencing has increasingly been used to detect AS/AT/APA, to identify novel isoforms, to predict non-coding RNA, and for gene fusion studies due to its long reads 14,15 . However, to date, there has been no report on a full-length transcriptome dataset of opium poppy, despite the recent release of its genome 1 . Herein, we construct the full-length transcriptome dataset of opium poppy using twenty-one pooled RNA samples from three different tissues and five different growth phases (Fig. 1). This produced 902,140 post-filter polymerase reads and 660,418 circular consensus sequence (CCS) reads (Table 1 and Fig. 2). After data processing, 566,746 full-length (FL) reads, 180,511 non-redundant isoforms and 61,856 unigenes (59,144 protein-coding unigenes and 2,712 non-coding unigenes) were obtained for functional annotation (Fig. 2). In addition, a total of 1,923,679,864 clean, paired-end short reads were produced ( Fig. 3 and Table 2). Gene expression levels were then determined using RSEM and converted into fragments per kb per million fragments (FPKM) value ( Fig. 4) 16,17 . The dataset reported here, provides an overview of the gene expression Fig. 1 Overview of the experimental design and the data processing pipeline. This study consisted of experimental design, RNA isolation, sequencing and data processing, annotation, and analysis, all of which are marked with different colors. Samples were divided into different growth phases, tissues, and cultivars. The numbers 1-5 indicate the five growth phases of opium poppy. Samples in the filling stage were used for either growth phases or tissues and cultivars. All 33 samples were subject to Illumina sequencing, and only samples from B1 were used for PacBio Sequel sequencing. www.nature.com/scientificdata www.nature.com/scientificdata/ levels, AS/AT/APA, and full-length transcript/unigene/mRNA of the opium poppy during key statuses. It can also be used to analyze the regulation mechanisms of BIAs biosynthesis according to its tissue-specific distribution, dynamic accumulation in different growth phases as well as BIAs diversity in different germplasm resources.

Methods
Plant material and experimental design. All the opium poppy for this study was cultured in our experimental plot. The original cultivar was named B1 and represents the original state of opium poppy, known as wild opium poppy. Other germplasm resources resulting from directed breeding in our institute were named B2, T, G and Z. Their characteristics and major BIAs content are summarized in Tables 3 and 4, respectively. Samples www.nature.com/scientificdata www.nature.com/scientificdata/ were divided into different tissues, cultivars, and growth phases ( Fig. 1). Samples from different tissues in B1 (root, stem, and leaves) and germplasm resources (B1, B2, T, G, ZB) were collected during filling stage on Jul 4, 2017. Samples in different growth phases (seedling, jointing, bolting, florescence, and filling) were collected from March to July of 2017. For each sampling point (eleven in all), three independent samples were collected. Unless otherwise mentioned, samples collected from opium poppy are all leaves. This produced 33 samples in all and their definitions of sampling are also listed in Table 3. The leaves from B1 at the filling stage were used as a shared sample for subsequent analysis: they were used as the filling stage sample when studying the growth rhythm of opium poppy, and also as a leaf sample in different tissues or a B1 sample in different cultivars. After washing and cleaning, the samples were immediately inundated with liquid nitrogen for 10 minutes and then stored in −80 °C freezers until use.
RNA extraction, library preparation, and sequencing. For each sample, TRIzol Reagent (Tiangen Biotech Co., Ltd., Beijing, China) was used to extract RNA following the protocol provided by the manufacturer. After cDNA synthesis, samples were subjected to phosphorylation, "A" base addition, and end-repair according to library construction protocol. Sequencing adapters were then added to both sizes of the cDNA fragments. After PCR amplification of cDNA fragments, the 150-250 bp targets were cleaned up. We then performed paired-end sequencing on an Illumina HiSeq X Ten platform (Illumina Inc, CA, USA) (Fig. 1). For PacBio Sequel sequencing, RNA from 21 samples of B1 was mixed in equal amounts for reverse transcription using the Clontech SMARTer PCR cDNA Synthesis Kit (TaKaRa, Dalian, China). In order to determine the optimal amplification cycle number for the downstream large-scale PCR reactions, PCR cycle optimization was employed (PrimeSTAR ® GXL DNA polymerase). Then the optimized cycle number was used to generate double-stranded cDNA. Large-scale PCR was performed for SMRTbell library construction (Pacific Biosciences). This include DNA damage repair, end repair, ligating sequencing adapters and removing fragments that failed to connect. Finally, the SMRTbell template was annealed to the sequencing primer, bound to polymerase, and sequenced on the PacBio Sequel platform using V2.1 chemistry (Pacific Biosciences) with 10-hour movies (Fig. 1). www.nature.com/scientificdata www.nature.com/scientificdata/ Data filtering, processing and yield. After sequencing, the raw reads were classified and clustered into a transcript consensus using the SMRT Link 5.1 pipeline (http://www.pacb.com/products-and-services/ analytical-software/smrt-analysis/). After adaptor removal and elimination of low quality regions, we obtain 902,140 post-filter polymerase reads (19.77 GB) with an average length of 22 kb ( Fig. 2 and Table 1). In order to improve the accuracy of sequencing, CCS reads were extracted from the subreads BAM file, which produced a total of 660,418 CCS reads with an average insert length of 2.61 kb (Fig. 2a,b). Briefly, CCS reads were extracted out of subreads.bam file with minimum full pass of 1 and a minimum read score of 0.8. CCS reads were then classified into full-length (FL) non-chimeric (NC), non-full-length (NFL), chimeras (C), and short reads based on cDNA primers and poly-A tail signal. Reads shorter than 50 bp were discarded. CCS reads with 5′ primer, 3′ primer and polyA tails were identified as FL reads, and 566,746 FL sequences ranging from 300 bp to 25,247 bp were obtained (Fig. 2c,d). Subsequently, the full-length non-chimeric (FLNC) reads were clustered by Iterative Clustering for Error Correction (ICE) software to generate the cluster consensus isoforms 18 . NFL reads were used by Arrow software to polish the obtained cluster consensus isoforms to obtain the final 55,114 FL polished high quality consensus sequences (accuracy ≥ 99%) 18 . Lordec was used to correct FL transcripts and CD-HIT was used to remove redundant sequences according to sequence similarity of high-quality transcripts 19,20 . Finally, 180,511 non-redundant isoforms and 61,856 unigenes were obtained for functional annotation (Fig. 2e,f). For Illumina paired-end RNA-seq, the low-quality reads (reads containing sequencing adaptors, reads containing sequencing primers, nucleotide with q quality score lower than 20) were removed. After that, a total of 1,923,679,864 clean, paired-end reads were produced (Table 2 and Fig. 3a). Illumina clean data was mapped onto our SMRT sequencing data using hisat2 v2.05 21 . RSEM was used to identify gene expression levels, which were then converted into an FPKM value 16,17 . DESeq R package was used to analyze differential expression 22 . Fold change ≥2 and adjusted P-value < 0.05 were set as threshold for significance of gene expression differences between the two samples (Fig. 4).  www.nature.com/scientificdata www.nature.com/scientificdata/ Major BIAs content measurement. To analysis the major BIAs content, some 0.5 g dry leaves in each sample of different germplasm lines at filling stag were employed according to our previous publicized method 23 . For details, the sample was sequentially extracted three times in methanol, with the help of ultrasonication, for 30 min at room temperature. Then, methanol extracts were combined and concentrated under reduced pressure conditions to a volume of 2 mL. At last, 10 μL concentrated extracts was subject to analysis using HPLC HPLC equipped with a reversed phase C18 column (XDB-C18, 5 mm; Agilent, USA) according to our publicized  Table 3. Description and definition of the main characteristics of opium poppy in each germplasm lines, growth phase, and tissue.
www.nature.com/scientificdata www.nature.com/scientificdata/ method. The BIAs content was show as weight (mg) in one Kilogram dry leaves (mg/kg) and listed in Table 4, and its raw data of individual measurements are available at Figshare 24 .
Functional annotation. Functional annotations of the novel genes were performed using BLAST searching against public databases such as Swiss-Prot, GO (Gene Ontology), and KEGG (Kyoto Encyclopedia of Genes and Genomics) 25   www.nature.com/scientificdata www.nature.com/scientificdata/

Data Records
The Illumina HiSeq X Ten data (different growth phases, cultivars and tissues) and PacBio Sequel sequencing data have been submitted to the Sequence Read Archive (SRA) of NCBI under accession numbers SRP173551 27 , SRP173565 28 , SRP173546 29 , and SRP173728 30 , respectively ( Table 5). The functional annotation and gene expression (RPKM) information of high-quality transcripts and unigenes are deposited in Figshare and Gene Expression Omnibus (GEO) in NCBI 24,31 . The differential gene expression data among different samples was also deposited in Figshare 24 .

Technical Validation
RNA quality control. The integrity of RNA sample was determined using the Agilent 2100 Bioanalyzer (Agilent Technologies, USA) and agarose gel electrophoresis. The purity and concentration of RNA samples were determined with the Nanodrop microspectrophotometer (Thermo Fisher Scientific, USA). For Illumina sequencing, 33 high-quality RNA samples (OD260/280 = 1.7~2.3, OD260/230 ≥ 2.0, RIN ≥ 7) were used to construct the sequencing library. The OD260/280, OD260/230, and RIN values for all RNA samples are listed in Table 6.
Quality evaluation of raw data. A high quality region finder was used to identify the longest region of a singly-loaded enzyme using a signal-to-noise ratio of 0.8 to filter out low-quality areas. In order to improve the accuracy of sequencing, the same polymerase reads were read multiple times in a closed loop, and the random error correction was then performed on the sequence read multiple times by the same insert fragment. This produced 660,418 CCS reads with a mean of 12.62 passes per read (Fig. 2a,b). FL reads were classified based on the location of and relationships between the 5′ primer, 3′ primer and polyA tail (Fig. 2c,d).
Assessment of sample composition. For Illumina paired-end RNA-seq data, we measured the correlation coefficient and quantitative saturation of gene expression among 33 samples (biological repetition and biological variation, Fig. 3). Correlation of expression levels correlation among samples is an important index to  Table 6. RNA sample quality used in this study.
www.nature.com/scientificdata www.nature.com/scientificdata/ test the reliability of experiments and the rationality of sample selection. Correlation coefficients close to 1 (red) indicate that the samples have a high similarity of expression patterns. If there is biological repetition in the sample, the correlation coefficient of biological repetition is usually higher. By calculating the relationship between the number of different full-length transcripts and the number of genes, we can observe whether the total number of genes has been measured up to the saturation level (Fig. 3b).