De novo transcriptome analysis of Cnidium monnieri (L.) Cuss and detection of genes related to coumarin biosynthesis

Sample preparation and RNA extraction

Whole C. monnieri plants were collected from Tongcheng city in Anhui Province, China in May 2018, with verbal permission of the Manager Zhongquan Fang (Anhui Bowen agriculture and Forestry Development Co., Ltd), and authenticated by Professor Dequn Wang (Anhui University of Chinese Medicine). Prior to the experiment, the plants were grown at a temperature of 22−26 °C/14−18 °C (day/night) and relative humidity of 65–80%. All samples were rinsed in ultrapure water, and leaves, roots, flowers and stems, which were harvested from three individual C. monnieri plants, were placed in centrifuge tubes, frozen in liquid nitrogen immediately and stored at −80 °C.

Total RNA was extracted from each tissue using RNA plant Plus Reagent (Tiangen, Beijing, China) according to the manufacturer’s instructions. The extracted RNA was checked using a NanoDrop 2000 (Thermo, CA, USA), and the concentration of the isolated RNA, 28S/18S and integrity (RIN) were estimated using an RNA Nano 6000 Assay Kit and the Agilent Bioanalyzer 2100 system (Agilent, CA, USA) (Table S1).

Determination of osthole content

Dried C. monnieri samples from flowers, roots, leaves and stems were used for isolation of osthole as previously reported (Jialei et al., 2014). Dried powder (0.1 g) from each sample was mixed with ethanol (25 mL), incubated for 2 h at room temperature (25 °C) and subjected to ultrasonic extraction for 60 min (100 W, 50 kHz). The supernatant was collected, and detected using UV-spectrophotometry at 322nm (SHIMADZU Coporation, kyoto, Japan). Osthole (purity >99.5%) (Aladdin Industrial Corporation, Shanghai, China) was used to construct a standard curve of the relationship between concentration and absorbance (Fig. S1A). The yield (percentage) of osthole was calculated as follows: ${\displaystyle \text{Yield}\left(\text{\%}\right)=\frac{\text{Osthole content of extraction}\left(\mathrm{g}\right)}{C.monnieri\text{tissue powder weight}\left(\mathrm{g}\right)}\times 100\text{\%}}$

cDNA library construction and sequencing

Total RNA was treated with RNase-free DNase I (TaKaRa, China) to remove DNA residues and then mixed with magnetic beads containing oligo (dT) to purify mRNA. After purification, the mRNA was fragmented into small pieces under elevated temperature. The cleaved RNA fragments were copied into first-strand cDNA using reverse transcriptase and random primers. Second-strand cDNA was synthesized using DNA Polymerase I and RNase H. These cDNA fragments had the addition of a single adenine at the 3′ end for subsequent ligation of adapters. The products were then purified and enriched using PCR amplification. The quality of each sample library was evaluated using an Agilent 2100 Bioanalyzer system (ABI, New York, NY, USA). Each cDNA library was sequenced using the BGISEQ-500 system at Beijing Genomics Institute (BGI) (Shenzhen, Guangdong province, China) with paired-end (PE) sequencing length of 100 bp.

Analysis of RNA-Seq data and unigene functional annotation

To ensure the reliability of transcript assembly, clean reads were obtained by discarding low-quality reads, reads containing poly-N and adapter sequences using SOAPnuke software (version 1.5.2) (parameters: -l 15 -q 0.2 -n 0.05 − i − A 0.25). In the absence of a reference genome, clean reads were subjected to transcriptome assembly using Trinity software (version 2.5.1) applying the following parameters: –min_contig_length 150 –CPU 8 –min_kmer_cov 3 –min_glue 3 –bfly_opts ’-V 5 –edge-thr =0.1 –stderr’ (Grabherr et al., 2011). Assembled transcripts were clustered to discard redundant sequences using the TGI clustering (TGICL) tool (parameters: -l 40 -c 10 -v 25 -O‘-repeat_stringency 0.95 -minmatch 35 -minscore 35′) (Pertea et al., 2003), leading to the identification of non-redundant sequences, called unigenes. The quality of assembled transcripts was assessed using a single copy orthologous database (BUSCO) (version 1.22).

Assembled unigenes were annotated against the NCBI nucleotide sequences (NT) (https://www.ncbi.nlm.nih.gov/nucleotide), NCBI non-redundant protein sequences (NR) (update date: 23 June 2019) (NR at https://www.ncbi.nlm.nih.gov/, update date: 23 June 2019), clusters of euKaryotic Orthologous Groups (KOG) (http://www.ncbi.nlm.nih.gov/KOG), Kyoto Encyclopedia of Genes and Genomes (KEGG) (update date: 1 July 2018) (http://www.genome.jp/kegg/) (Kanehisa et al., 2019), and a manually annotated and reviewed protein sequence database (SwissProt) (update date: 25 May 2018) (http://www.uniprot.org/). Assembled unigenes were annotated against NR database using the software BLAST (version 2.2.23, E-value ≤ 1e−5). And the results of NR annotation were used for Gene Ontology (GO) (update date: 1 January 2019) annotation using Blast2GO (version v2.5.0) program (Conesa et al., 2005). Pfam (update date: September 2018) annotation was carried out using HMMER (version 3.0).

Quantitative real-time PCR (qRT-PCR) analysis of key genes in coumarin biosynthesis

Total RNA from each sample was subjected to reverse transcription using SuperScript™ III First-Strand Synthesis SuperMix (Invitrogen, US). qRT-PCR analysis of each gene was performed on a QuantStudio^® Real-time PCR system (Life Technologies, US) using Power SYBR^® Green PCR Master Mix (Roche, China). The actin gene (CL3748.Contig7) of C. monnieri was used as a reference. Primers for all selected genes were designed using Premier 6.0 and are listed in Table S2. Each qRT-PCR contained 1µl of diluted cDNA, 0.5 µl of each primer (10 µM) and 10.0 µl of Power SYBR^® Green Master Mix. Reactions were conducted under the following conditions: denaturation at 95 °C for 1 min, followed by 40 cycles of 95 °C for 15 s and 63 °C for 25 s. Each qRT-PCR was performed using three biological replicates. The 2^−ΔΔCt method was used to calculate the relative expression level of genes (Livak & Schmittgen, 2001).

Analysis of differentially expressed genes (DEGs)

After assembly, clean reads were mapped to unigenes using Bowtie2 (version 2.2.5) (Langmead & Salzberg, 2012), and then gene expression levels were calculated with RSEM (version 1.1.12) (Dewey & Li, 2011). Since we had no RNAseq bioreplicates, we used a Poisson distribution method to identify candidate genes that might be differentially expressed. Both fold change (≥2.0) and false discovery rate (≤0.001) are used as criteria to narrow down the list of differentially expressed genes between flowers and other tissues (Bullard et al., 2010; Lee et al., 2008; Li et al., 2013). qRT-PCR with bioreplicates was used to confirm differential expression of a subset of these. GO and KEGG analyses were again performed on the DEGs following the method described by Audic (Audic & Claverie, 1997).

Identification of transcription factors

The open reading frames (ORFs) of unigenes were initially detected using getorf (version EMBOSS:6.5.7.0) (Rice, Longden & Bleasby, 2000) and further aligned to the plant transcription factor database (PlntfDB) through comparison with Pfam 23.0 using the hmmsearch method (Mistry et al., 2013). The function of unigenes was identified according to the characteristics of each transcription factor family described by PlantfDB.

Statistical analysis

The contents of osthole and all data from qRT-PCR were presented as mean ± standard deviation. Statistics were done by GraphPad Prism 6.0. Heatmap and boxplot were performed using R software (version 3.4.4). The one-way analysis of variance (ANOVA) was used to compare the data between groups with Duncan t test. P value <0.05 was considered to be statistically significant.

Osthole content determination in different tissues of C. monnieri

We extracted osthole from dried leaves, stems, flowers and roots of C. monnieri. Osthole content was highest in flowers (1.372%) and lowest in stems (0.283%) (Fig. S1B). Osthole quantification revealed significant differences among different tissues, with this active ingredient mainly distributed in flowers and secondly in leaves.

RNA sequencing and de novo assembly

Four mRNA libraries were extracted from leaf, stem, flower and root tissue of C. monnieri using BGISEQ-500 sequencing technology. After eliminating low-quality reads, we obtained a total of 10.6 Gb, 10.47 Gb, 10.81 Gb and 10.38 Gb of clean reads from 83.42 Mb, 81.65 Mb, 84.27 Mb and 82.14 Mb sets of raw data from leaf, stem, flower and root tissue, respectively, with all Q30 values greater than 85% (Table S3). The clean reads were assembled using Trinity and the TGI clustering tool (TGICL). After assembly, a total of 173,938 unigenes were obtained from the four libraries with a mean length of 1,272 bp, GC content of 38.79% and N50 value of 2,121 bp (Table S4). Of these unigenes, 45.80% (79,658) and 64.13% (111,562) exceeded 1,000 bp and 500 bp in length, respectively (Fig. S2). The quality of the assembled transcripts was assessed using a single copy orthologous database (BUSCO), and 98% of the unigenes showed a perfect match (Fig. S3).

Unigene functional annotation and expression overview

Identification of assembled unigenes was carried out using the Basic Local Alignment Search Tool (BLAST) (e <10⁻⁵) against seven public databases (NR, NT, Swissprot, KEGG, KOG, Pfam and GO), and a total of 110,992 (63.81%), 95,603 (54.96%), 80,216 (46.12%), 86,568 (49.77%), 87,523 (50.32%), 79,335 (45.61%) and 56,429 (32.44%) unigenes were aligned, respectively (Table 1). We annotated 56,429 (32.44%) unigenes using GO terms classified into 58 subcategories within three standard categories (molecular functions, biological processes and cellular components) (Fig. S4). We focused on the categories of molecular functions and biological processes. “Catalytic activity” and “binding” were the most enriched terms among the molecular functions, while “cell process” and “metabolic process” were the top subcategories of biological processes. According to Venn diagram analysis, 54,846 (31.53%) unigenes were co-annotated in five databases (Fig. 1). Additionally, 110,992 unigenes were annotated in the NR database and 55.48% of these annotated unigenes were mapped to Daucus carota subsp. sativus (82.2%) and others (17.8%) (Fig. S5).

We determined expression values of transcripts with fragments per kilobase of transcript per million mapped reads (FPKM) >1 for each tissue and found 55,421, 52,787, 48,287 and 44,717 expressed unigenes in flower, leaf, stem and root tissues, respectively (Fig. 2A). The overall expression level of unigenes in flowers was higher than that in leaves, stems or roots (Fig. 2B).

Identification of genes involved in coumarin biosynthesis

To identify the most significant biological processes in the leaf, stem, flower and root tissues of C. monnieri, 86,568 (49.77%) unigenes were annotated and assigned to 135 KEGG pathways (19 subcategories) (Fig. S6, Table S5). Fourteen pathways were associated with the biosynthesis of other secondary metabolites, of which phenylpropanoid biosynthesis pathway was most enriched with annotated genes (Fig. 3). Phenylpropanoids comprise a large group of natural products including coumarins, flavonoids and lignins. In our transcriptome data, 1,575 unigenes were assigned to the phenylpropanoid biosynthesis (ko00940) pathway based on KEGG pathway classification. Of these, 853 unigenes encode enzymes participating in the coumarin biosynthesis pathway: phenylalanine ammonia-lyase (PAL, 14 unigenes), β-glucosidase (BGA, 507 unigenes), cinnamate 4-hydroxylase (C4H, 2 unigenes), 4-coumarate–CoA ligase (4CL, 49 unigenes), shikimate hydroxycinnamoyl transferase (HCT, 70 unigenes), 5-O-(4-coumaroyl)-D-quinate 3′-monooxygenase (C3′H, 95 unigenes), caffeoyl-CoA O-methyltransferase (COMT, 33 unigenes) and feruloyl-CoA ortho-hydroxylase (F6′H1, 83 unigenes) (Table 2). These data enabled us to outline coumarin biosynthesis based on genes encoding enzymes with FPKM >1 (Fig. 4).

Expression analysis of key enzyme genes

We selected six unigenes encoding five key enzymes and used quantitative real-time PCR (qRT-PCR) to determine their relative expression levels in the leaf, stem, flower and root tissues of C. monnieri (Fig. S7). Relative expression of the CL3766.Contig1 (C4H), CL16061.Contig1 (4CL) and CL8379.Contig6 (F6H1) genes was highest in root tissue, whereas relative expression of the CL11351.Contig3 (PAL) and CL2358.Contig9 (BGA) genes was highest in flower tissue. Unigene 20867 (BGA) had a high expression level in leaf tissue. The expression levels determined by qRT-PCR were consistent with the FPKM values.

Identification of DEGs

We identified 12,635 unigenes uniquely expressed in flowers and 77,762 shared unigenes exhibiting expression in all four tissues (Fig. 5A). Total DEGs were identified among the four tissues using gene FPKM values. A comparison between flowers and leaves revealed a total of 46,460 DEGs, of which 26,944 were up-regulated and 19,516 were down-regulated in flowers (Fig. 5B). A further comparison of flowers with roots revealed 47,901 DEGs, of which 29,754 were up-regulated and 18,147 were down-regulated in flowers (Fig. 5B). Comparison among flowers and stems revealed 36,666 DEGs, of which 21,683 were up-regulated and 14,983 were down-regulated in flowers (Fig. 5B).

All DEGs were annotated in the KEGG database to further describe and evaluate their biological functions. This analysis revealed 32,060 DEGs in flower versus leaf associated with 137 pathways; these DEGs were mainly enriched in the “plant hormone signal transduction,” “plant-pathogen interaction,” “circadian rhythm - plant” and “photosynthesis - antenna proteins” pathways (Table S6). In flower versus root, 34,786 DEGs were annotated to 137 pathways and primarily enriched in “biosynthesis of secondary metabolites,” “photosynthesis,” “photosynthesis - antenna proteins” and “phenylpropanoid biosynthesis” pathways (Table S7). In flower versus stem, 27,513 DEGs were annotated to 137 pathways and primarily enriched in “plant hormone signal transduction,” “indole alkaloid biosynthesis,” “glucosinolate biosynthesis” and “phenylpropanoid biosynthesis” pathways (Table S8).

A total of 18,881 up-regulated DEGs showed flower-specific expression with |log2 (fold changes) | >2. We inferred the nature of these DEGs via GOSlim functional analysis. Sequence homology revealed that each of the 18,881 flower-specific DEGs mapped to at least one ontology term, including 2,814 for biological processes, 3,232 for cellular components and 4,508 for molecular function; many genes were enriched in “cellular process” and “catalytic activity” subcategories within the biological processes and molecular function categories, respectively (Fig. S8).

To evaluate the biological functions of these DEGs, the 18,881 flower-specific DEGs were annotated in the KEGG database. KEGG enrichment analysis showed that these DEGs were significantly enriched in the “ribosome,” “starch and sucrose metabolism,” “pentose and glucuronate interconversions” and “phenylpropanoid biosynthesis” categories (Fig. 6). Additionally, we identified 129 up-regulated flower-specific DEGs with —log2 (fold change)—>2 related to coumarin biosynthesis, including BGA (99 DEGs), 4CL (2 DEGs), HCT (13 DEGs), COMT (7 DEGs) and F6′H1 (8 DEGs) (Table 3, Table S9).

Identification of transcription factors (TFs)

Transcription factors (TFs) play an important role in regulating secondary metabolites by modulating the expression of genes related to biosynthetic pathways. We identified 3,799 unigenes encoding putative TFs in the C. monnieri transcriptome database, including 1,169, 1,032 and 920 unigenes up-regulated in flowers compared with leaves, stems and roots, respectively (Table 4). The major types of TF identified belonged to the MYB (v-myb avian myeloblastosis viral oncogene homolog) (564 unigenes), AP2-EREBP (APETALA2 and ethylene-responsive element binding proteins) (285 unigenes), bHLH (Basic helix-loop-helix) (274 unigenes), NAC (222 unigenes), C3H (Cys3His zinc finger) (193 unigenes), WRKY (191 unigenes), and C2H2 (C2H2-type transcription factor) (173 unigenes) TF families. Moreover, a total of 15 unigenes encoding 3 TFs (MYB TF, GRAS TF and HAT TF) participated in the biosynthesis of phenylpropanoid (Table S10). Of these, 7 unigenes encoding MYB TF were up-regulated with flower-specific expression; these unigenes are shown in red in Table S10.