Transcriptome analysis and identification of key genes involved in 1-deoxynojirimycin biosynthesis of mulberry (Morus alba L.)

Plant materials

Mulberry leaf samples were obtained from a mulberry plantation in Jiangsu University (Zhenjiang, China). Mulberry tender leaves from top of branches of the same tree were collected on July 31, 2016 (M7), November 15, 2016 (M11) for the first time, and on July 25, 2017 (Ma7), November 15, 2017 (Ma11) for the second time—DNJ contents were significantly different between the two samples in each collection time, 2017 samples were collected from a different tree than 2016 samples. The collected samples were immediately frozen in liquid nitrogen, and then stored at −80 °C until transcriptome analysis. Fresh mulberry leaves were collected from April 15, to November 25, from the same location on each branch (from the top 6–10 leaves); the samples were collected every 10 days, and each independent sample was made up with ten mulberry leaves, the leaves were dried for DNJ content determination.

Determination of DNJ content

1-Deoxynojirimycin was extracted as previously described (Ouyang, Chen & Li, 2005), with minor modifications. Fresh leaves were dried at 60 °C and ground to powder prior to analysis. A total of 200 micrograms of the powder was added to 10 mL aqueous 0.05M HCl, extracted by vortexing for 20 s, and centrifuged at 9,000×g for 15 min. The supernatant was collected in a 50 mL volumetric flask, the sediment was extracted a second time, and the final supernatant was added to the volumetric flask and diluted to 50 mL with water. A 10 μL aliquot of the extracting solution was mixed with 10 μL of 0.4M potassium borate buffer (pH 8.5), and 20 μL of 5 mM fluorenylmethyloxycarbonyl chloride (FMOC-Cl) in acetonitrile was added to the mixture and allowed to react at 20 °C for 20 min. Then, 10 μL of 0.1M glycine was added to neutralize the excess FMOC-Cl and 950 μL of 0.1% (v/v) aqueous acetic acid was added to stabilize the reaction compounds. The reactant solution was filtered through a 0.22 μm Millipore filter (Millipore Corp., Billerica, MA, USA), and 10 μL of the filter liquor was used for subsequent analysis by high performance liquid chromatography (HPLC) with a fluorescence detector (FLD). An Agilent 1260 Infinity HPLC system (Agilent Technologies Inc., Palo Alto, CA, USA) with a HIQ Sil C18-10 column (size 4.6 × 250 mm, no. 00V00194; Kya Technologies, Tokyo, Japan), and a 1260 FLD (excitation 254 nm, emission 322 nm). The mobile phase consisted of acetonitrile and 0.1% aqueous acetic acid (1:1, v/v) at one mL/min for 20 min, and the column was maintained at 25 °C. The data were analyzed using an Agilent 1260 Open Lab workstation. All the determinations were performed in triplicate.

Total RNA extraction, cDNA library construction, and RNA-Seq

Total RNA was extracted using TRIzol reagent (Sangon Biotech, Shanghai, China). RNA quality and quantity were analyzed using 1% agarose gel electrophoresis and a Qubit 2.0 RNA kit (Life Technologies, Carlsbad, CA, USA). Independent total RNA extracts were prepared from multiple leaves (more than three leaves) and equal amounts of RNA from each sample were used to generate RNA-Seq libraries. The library was constructed using the VAHTSTM mRNA-seq v.2 Library Prep kit for Illumina (Vazyme Biotech Co., Nanjing, China). Briefly, one μg total RNA was purified using poly-T oligo-conjugated magnetic beads, and cleaved RNA was used for cDNA synthesis with random primers. The cDNA fragments were purified, the ends blunted, and a poly-A tail added, followed by adapter ligation. Polymerase chain reaction (PCR) amplification of cDNA fragments was carried out using PCR Primer and Amplification Mix (Vazyme Biotech Co., Ltd., Nanjing, China). The reaction conditions were 15 cycles (98 °C for 30 s, 98 °C for 10 s, 60 °C for 30 s, 72 °C for 30 s, 72 °C for 5 min, and holding at 4 °C). Reaction products were purified using a Qubit 2.0 DNA kit, and double-stranded cDNAs were used for HiSeq 2500 paired-end sequencing.

De novo transcriptome data processing and assembly

A perl script was written to remove vector sequences and PolyA (T) tails from the raw sequence data using cutadapt software (http://pypi.python.org/pypi/cutadapt/1.2.1, version: 1.2.1). Low-quality reads (<20 bases) and those with N-portions >35 bp were removed using Prinseq software (http://prinseq.sourceforge.net/, version: 0.19.5). Reads were merged after eliminating repeats. The high-quality reads were assembled using Trinity software (http://trinityrnaseq.github.io/, version: r20140717) to construct unique transcripts (Grabherr et al., 2011).

Function annotation and classification

Unigenes were compared using the basic local alignment search tool (BLAST) against the National Center for Biotechnology Information non-redundant nucleotide (NT, ftp://ftp.ncbi.nlm.nih.gov/blast/db/nt.*tar.gz, using BlastN) and non-redundant protein (NR, ftp://ftp.ncbi.nlm.nih.gov/blast/db/nr.*tar.gz, using BlastX) databases (January 2013; E-value ≤ 1e⁻⁵). SWISS-PROT (downloaded from the European Bioinformatics Institute on January 2013; E-values ≤ 1e¹⁰, using BlastX) was used for annotation and classification. The Clusters of Orthologous Groups of Proteins database (COG; E-values ≤ 1e¹⁰, using rpsBlast), (KEGG, release 58; E-values ≤ 1e¹⁰), and InterProScan Release 36.0 annotated protein domains and functional assignments were mapped onto Gene Ontology (GO, http://www.geneontology.org/) by using the BlastX algorithm.

Identification of differentially expressed unigenes

Unigenes that were differentially expressed between the two cDNA libraries were identified by DESeq and edgeR software (http://www.bioconductor.org/) according to a previously reported method (Anders & Huber, 2010). DEGs were considered as genes having differential expression with a P-value ≤ 0.01 and an expression ratio of ≥2. DESeq software was used to analyze genes with a q-value < 0.001 and fold change >2.

Quantitative real-time (qRT)-PCR analysis of differentially expressed unigenes

A total of 20 unigenes (Table S1) that were significantly differentially expressed were randomly selected to validate the accuracy of RNA-Seq data using qRT-PCR. Total RNA was reverse transcribed using the RevertAid™ First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, Waltham, MA, USA). The real-time PCR experiments were performed using a Fast Start Essential DNA Green Master kit and a Light Cycle 96 Real-Time PCR system (Roche, Basel, Switzerland). The primers were designed using Primer Premier 5.0 (Premier Biosoft Ltd., Palo Alto, CA, USA) and were obtained from Sangon Biotech (Shanghai, China) (Table S2A; Tables 2B and 2C). The RT-PCR program comprised a reverse transcription step at 48 °C for 30 min and a Taq polymerase activation step at 95 °C for 300 s, followed by PCR: 45 cycles at 95 °C for 15 s, 61 °C for 20 s, and 72 °C for 30 s, followed by a melting cycle. The 2^−ΔΔCt (Livak & Schmittgen, 2001) method was used to calculate the relative changes in target gene expression from the qRT-PCR results, using the mulberry β-actin gene as the reference gene. All experiments were carried out in technical and biological triplicate (independent leaves).

Biological validation of RNA-Seq results

The same strategy was used to perform transcriptome studies on samples (M7, M11) and (Ma7, Ma11). BLASTN was used to compare the coverage of sequences. The number and coverage of DEGs were screened to evaluate the biological repeatability of the experiment.

Statistical analysis

Data are expressed as the mean ± standard deviation (SD). Data were evaluated for statistically significant difference by using Student’s t-test. SPSS 13.0 (IBM SPSS, Armonk, NY, USA) was used for analysis, and P < 0.05 was considered statistically significant. Correlation analysis was conducted using GraphPad Prism 5.0 software (San Diego, CA).

DNJ content in mulberry samples

To investigate the relationship between DNJ content and gene expression, the DNJ contents of different developmental staged mulberry plants (M. alba L.) were determined by HPLC-FLD. The DNJ content in M7 was 3.43 ± 0.15 mg/g dry weight, as compared to 0.19 ± 0.05 mg/g dry weight in M11 (P < 0.001; Fig. 2).

Mulberry transcriptome sequencing and de novo assembly

As shown in Table 1, mulberry leaf samples of two different stages of development were sequenced using the Illumina Hiseq2500 system, which generated 105,189,326 and 65,166,690 raw reads for M7 and M11, respectively. The average length of a single read was 150 bp; after filtering out low-quality reads, 104,988,228 and 64,878,089 reads (99.81% and 99.56% of raw sequences, respectively) were obtained for M7 and M11, respectively, with Q20 values of 97.52% and 94.42%, respectively. High-quality reads were assembled from 112,481 unique transcripts with a mean length of 766 bp and an N50 value of 1,392. Transcript lengths ranged from 201–18,301 bp (Fig. S1).

Functional annotation of the mulberry transcriptome

Given the lack of a reference genome for mulberry, 112,481 transcripts were searched using BLAST against nine public databases (Conserved Domains Database, Eukaryotic Orthologous Groups (KOG), NR, NT, Protein Families database, Swiss-Prot, TrEMBL, GO, and KEGG), and 85,949 transcripts were annotated (Fig. S2). The species distribution map according to annotation by BLAST against the NR database is shown in Fig. 3A, and includes M. notabilis (58.93%), Vitis vinifera (4.06%), Ziziphus jujuba (3.87%), Hordeum vulgare subsp. vulgare (1.45%), Dorcoceras hygrometricum (1.26%), Cajanus cajan (1.14%), Prunus persica (1.03%), Ricinus communis (1.03%), Citrus sinensis (1.02%), and Pyrus bretschneideri (0.95%). A total of 20,725 unigenes were annotated by BLAST against the COG/KOG database (Fig. S3), and 189,838 transcripts were categorized as a biological process, cellular component, or molecular function in the GO classification (Fig. 3B); 8,655 were assigned to pathway branches for cellular processes, environmental information processing, genetic information processing, metabolism, organismal systems, and human diseases, and were further divided into 323 KEGG pathways (Table S3).

Analysis of transcripts differently expressed between samples M7 and M11

The homologies of expressed unigenes from the M7 and M11 datasets are shown in Fig. 4 and Table S4; 68,235 and 66,751 transcripts were expressed in M7 and M11, respectively, with 42,994 homologous unigenes and 25,241 and 23,757 nonhomologous genes in M7 and M11, respectively. Unigenes that were differentially expressed between the two cDNA libraries were compared, with q-value < 0.001 and fold change >2 used as thresholds. A total of 11,318 unigenes were found to be differentially expressed, with 6,606 and 4,712 unigenes being up- and downregulated, respectively. The differentially expressed unigenes were annotated as biological processes, cellular components, and molecular functions, and were among the top 50 terms enriched in the GO classification (Fig. 5). Unigenes showing the greatest difference in expression were associated with biological processes, followed by genes corresponding to cellular component. Genes that fell under the molecular functions category did not exhibit large differences in expression between the two samples, and were mostly associated with sequence-specific DNA binding, nucleic acid binding, and oxidoreductase activity. According to the KEGG annotation, these 11,318 differentially expressed unigenes were assigned to 294 KEGG pathways (Table S5), with the highest distribution category being carbon metabolism and biosynthesis of amino acids, followed by oxidative phosphorylation process.

Identification of candidate genes involved in mulberry DNJ biosynthesis

A DNJ alkaloid biosynthetic pathway biosynthetic pathway was outlined on the basis of differentially expressed transcripts and KEGG pathway assignments. Firstly, lysine was generated using aspartic acid as a substrate, and then the lysine was used as a substrate to generate the DNJ alkaloid; 10 enzymes, encoded by 38 genes, were annotated (Table 2; Table S6). Eight enzymes were annotated in the process of aspartic acid to lysine. The homoserine dehydrogenase (EC: 1.1.1.3) can catalyze L-aspartate 4-semialdehyde to L-homoserine. We did not analyze the enzyme because it does not contribute to the production of lysine. Only two enzymes, LysC (EC: 2.7.2.4) (aspartate kinase) and DapA (EC: 4.3.3.7) (4-hydroxy-tetrahydrodipicolinate synthase), exhibited significantly low expression. LysC (EC: 2.7.2.4), was encoded by 10 transcripts. Two transcripts showed significant downregulation, suggesting that the gene may be related to DNJ alkaloid biosynthesis; four transcripts were annotated as DapA with one showing significant downregulation, suggesting that this transcript may reflect DNJ biosynthesis. Another enzyme, aspartate semialdehyde dehydrogenase (ASD) (EC: 1.2.1.11), was encoded by two transcripts, whereas DapB (EC: 1.17.1.8) was encoded by one transcript. Three transcripts were annotated as ALD/AGD (EC: 2.6.1.83), of which one was downregulated significantly. Three transcripts were annotated as DapF (EC: 5.1.1.7). LysA (EC: 4.1.1.20) was encoded by two transcripts.

In contrast, two downregulated enzymes, LdcC (EC: 4.1.1.18) and AOC2/3 (EC: 1.4.3.21), were assigned to the pyridine ring backbone biosynthetic pathway. LdcC (EC: 4.1.1.20) was encoded by one transcript, showing significant downregulation, suggesting that this transcript may be involved in DNJ biosynthesis. AOC2/3 was represented by 12 transcripts, three of which were significantly downregulated and one of which was significantly upregulated. These four transcripts may thus be involved in the DNJ biosynthesis process. Through the above analysis, nine DEGs were obtained as candidate genes for the DNJ alkaloid biosynthetic pathway. The correlation between the expression level of these genes and the DNJ content was further analyzed to screen the key enzyme genes.

In the conversion of 1-piperidinene to DNJ, CYP450 and methyltransferase genes are speculated to play important roles. In this study, 103 CYP450-related genes in total, including 20 DEGs, were identified (Table S7). A total of 11 DEGs were upregulated and nine were downregulated. In total, 110 methyltransferase genes, including 17 DEGs (Table S8), eight DEGs were upregulated and nine were downregulated. The correlation between DEGs expression and DNJ content was further analyzed.

Validation of the accuracy of RNA-Seq results

A total of 20 transcripts were randomly selected to validate the accuracy of the RNA-Seq data using qRT-PCR. The validation results shown in Table 3 suggest that the RNA-Seq data are reliable, two fold-change values were greater than one or less than one at the same time, indicating that they have the same trend. This method could be used to identify the DEGs involved in specific pathways of metabolism, and thus, that the data could be used for subsequent analyses. The gene IDs and fragment per kilobases per million mapped reads (FPKM) values are shown in Tables S1, S6 and S9.

The biological replication result of mulberry RNAseq was shown in Table 4, nine DEGs were found in M7 and M11, and seven DGEs between Ma7 and Ma11 were identified. The BLAST result showed that seven DEGs matched each other with 99–100% coverage, and the changes of the matching sequence were consistent.

Correlation between candidate gene expression and DNJ content of M7 and M11

According to a previous study by our research group, DNJ content exhibits certain regularity with respect to the different growth stages (Ouyang & Chen, 2004); specifically, from April to November the level of DNJ showed a trend from rise to decline, with the peak level occurring at July and August. In this study, mulberry samples were collected from April 15, 2016, to November 25, 2016 for DNJ content determination. The results showed that the DNJ content was higher in M7 than in M11. DNJ contents, the expression of the nine candidate genes, and the relationship between the two are shown in Fig. 6 and Table 5. The expression of three of the transcripts (c83396_g1, c47185_g1, and c48882_g1) showed a decreasing trend that was consistent with the DNJ content. The correlation results indicated that c83396_g1 expression was significantly positively correlated with DNJ content (P < 0.001, Fig. 6F); c47185_g1 and c48882_g1 expression levels were positively correlated with DNJ content (P < 0.01, Figs. 6G and 6H). The expression of the remaining six transcripts did not correlate with DNJ content. The transcripts c8863_g1 and c44849_g1 showed no significant change in August, but increased in September, and fluctuated during the late stages (Figs. 6B and 6E); c40311_g1 was increased from August to October, and decreased during the late stages (Fig. 6C); similar trends were observed for c43595_g1, c47618_g1, and c32423_g1 (Figs. 6D, 6I and 6J).

The correlation analysis showed that the transcript level of five CYP450 genes was significantly correlated to DNJ content (Fig. 7; Table 6). The level of the transcript c44921_g1 was significantly and negatively correlated with DNJ content (P < 0.001, Fig. 7E). The transcript levels of c39525_g1, c29603_g1, and c33960_g1 showed that the expression was significantly postively correlated with DNJ content (P < 0.001, Figs. 7G, 7H and 7K). c105322_g1 expression levels was positively correlated with DNJ content (P < 0.01, Fig. 7S). The correlation analysis results showed that two methyltransferase trancripts were significantly correlated to DNJ content (Fig. 8; Table 7). The expression of the transcripts c41333_g1 (Fig. 8B) and c43454_g1 (Fig. 8Q) was significantly and postively correlated with DNJ content (P < 0.01, Fig. 8B). The expression of transcript c43454_g1 was significantly postively correlated with DNJ content (P < 0.001, Fig. 8Q). No significant correlations were found between other genes and DNJ content.