Whole-genome sequencing and analysis of the Chinese herbal plant Gelsemium elegans.

BACKGROUND: Gelsemium elegans (G. elegans) (2n = 2x = 16) is genus of flowering plants belonging to the Gelsemicaeae family.
METHOD: Here, a high-quality genome assembly using the Oxford Nanopore Technologies (ONT) platform and high-throughput chromosome conformation capture techniques (Hi-C) were used.
RESULTS: A total of 56.11 Gb of raw GridION X5 platform ONT reads (6.23 Gb per cell) were generated. After filtering, 53.45 Gb of clean reads were obtained, giving 160 × coverage depth. The de novo genome assemblies 335.13 Mb, close to the 338 Mb estimated by k-mer analysis, was generated with contig N50 of 10.23 Mb. The vast majority (99.2%) of the G. elegans assembled sequence was anchored onto 8 pseudo-chromosomes. The genome completeness was then evaluated and 1338 of the 1440 conserved genes (92.9%) could be found in the assembly. Genome annotation revealed that 43.16% of the G. elegans genome is composed of repetitive elements and 23.9% is composed of long terminal repeat elements. We predicted 26,768 protein-coding genes, of which 84.56% were functionally annotated.
CONCLUSION: The genomic sequences of G. elegans could be a valuable source for comparative genomic analysis in the Gelsemicaeae family and will be useful for understanding the phylogenetic relationships of the indole alkaloid metabolism.

Introduction

Gelsemium is a genus of flowering plants belonging to the Gelsemicaeae family. The genus comprises 3 species: the Asian Gelsemium elegans (G. elegans) and two North American species, Gelsemium sempervirens and Gelsemium rankinii. G. elegans (National Center for Biotechnology Information Taxonomy ID: 427660) is also known as Gou Wen, Da Cha Yao or Duan Chang Cao in Chinese1, 2. This herb is extremely famous for murdered the Yan Emperor of China who was also called Shen Nong in the Chinese mythology “Shen Nong tastes hundreds of grasses”. In this myth, Shen Nong diligently tasted all manner of flora for people to eat or used for medicine. But one day, he tasted Duan Chang Cao which has the yellow flower, and this poison was so terrible that he died quickly. He sacrificed himself to save humanity, so people call him the “Bodhisattva of medicine”, and people forever commemorate him. Just as the myth described, this species is widely distributed in the Fujian, Guangxi, Hunan and Guizhou provinces of China and in southeastern Asia (Fig. 1). It has been used as an herbal medicine for the treatment of rheumatoid arthritis, neuropathic pain, spasticity, skin ulcers and cancer for many years3, 4 and the whole plant has been widely added to animal feed for livestock. To date, more than 200 compounds, including indole alkaloids, iridoids, and steroids, have been isolated and identified from G. elegans5, 6.

Figure 1

Example of the Gelsemium elegans (Gou Wen or Duan Chang Cao). (A) Natural habitat of Gelsemium elegans (image from Qi Tang). (B) Gelsemium elegans (image from Yisong Liu). (C) The flower of Gelsemium elegans (image from Qi Tang).

Previous studies on the crude and purified alkaloids of G. elegans have demonstrated that this species possesses anti-inflammatory, immunomodulating, analgesic, anxiolytic, anti-tumor8, 9, and neuropathic pain-relieving properties. Indole alkaloids such as gelsemine, koumine, humantenine, gelsemicine and gelsenicine are the major active components of Gelsemium. Gelsemine and koumine are the principal alkaloids in G. elegans, and their toxicity is relatively weak. Gelsenicine {[LD50 = 0.128 mg/kg, mice (i.p.); 0.26 mg/kg, rat (i.p.); 0.15 mg/kg, rat (i.v.)], which was found in a lesser amount, was the most toxic alkaloid in G. elegans. Gelsenicine was also the most toxic compound in G. sempervirens12, 13. The typical symptoms of gelsenicine intoxication include chest distress, asphyxia, dizziness, tonic convulsions, limb paralysis, and difficulty breathing. Severe gelsenicine poisoning can cause multiple organ failure leading to death. Therefore, the actual bioactive components of G. elegans have attracted attention from chemists, pharmacologists and toxicologists due to their complex structural features and multiple biological effects.

Despite the considerable pharmaceutical importance of G. elegans, the genomic information available for this species is limited, which has hindered its utilization. Former results suggest that the Oxford Nanopore Technologies (ONT) can be used to quickly and cost-effectively generate informative assemblies15, 16, and a combination of sequencing and mapping data often leads to improved assemblies and is potentially more cost effective than sequencing alone. For example, the cottons and human genomes were assembled using a combination of long reads and Hi-C (high-throughput chromosome conformation capture techniques)-based data, have remarkably high quality with long contig (contig N50 of 18.7 and 26.8 Mb, respectively), chromosome length scaffolds (scaffold N50 of 87 and 60.0 Mb) and nearly 100% sequence fidelity. Here we report a high-quality reference genome for G. elegans using ONT technology and Hi-C map to cluster the majority of the assembled contig onto 8 pseudo-molecules, which is expected to facilitate and expand its use.

Materials and methods

Sampling and sequencing

All samples were collected from Liucheng city, Guangxi Province, China (N24°39′15.96″, E109°14′25.37″). Genomic DNA was extracted from leaves of a single plant using the Plant Genomic DNAkit (Qiagen, San Diego, CA, USA). Genomic DNA sample was further purified for ONT sequencing with the Zymo Genomic DNA Clean and Concentrator-10 column (Zymo Research, Irvine, CA, USA). The purified DNA was then prepared for sequencing following the protocol in the genomic sequencing kit SQK-LSK108 (ONT, Oxford, UK). Single-molecule real-time sequencing of long reads was conducted on a GridION X5 platform (Oxford Nanopore Technology, OX4 4DQ, Oxford, UK) with 9 flow cells. A total of 56.11 Gb of genomic data (6.23 Gb per cell) with an average read length of 14.59 kb was generated after quality filtering, from which the longest reading is 153.6 kb (Supporting Information Table S1). Compared with other sequencing platforms, Nanopore platform reading length has more advantages. In addition, a separate paired-end (PE) DNA library with an insert size of 400 bp (amplification by 8 PCR cycles) was constructed and sequenced using the Illumina platform (PE150) to enable a genome survey, and a total of 53.2 Gb of raw data was collected (Supporting Information Table S2).

For RNA-seq, total RNA from 12 samples were extracted from leaf, root, stem and flower of one G. elegans, using the QIAGEN RNeasy Plant Mini Kit (QIAGEN, Hilden, Germany). The cDNA library was prepared using the TruSeq Sample Preparation Kit (Illumina, CA, USA), and paired-end sequencing with 150 bp was conducted on a HiSeq X Ten platform (Illumina, CA, USA). A total of 135.9 Gb clean data were obtained (Supporting Information Table S3).

Genome size and heterozygosity estimation

The genome size of G. elegans was estimated by the k-mer method using sequencing data from the Illumina DNA library. Quality-filtered reads were subjected to 17-mer frequency distribution analysis using the Jellyfish program. The genome size (G) of G. elegans was estimated using the following formula: G = k-mer number/average k-mer depth, where k-mer number = total k-mers-abnormal k-mers. The count distribution of 17-mers followed a Poisson distribution, with the highest peak occurring at a depth of 120 (Supporting Information Table S4 and Fig. S1). The estimated genome size was 338,031,359 bp, and the heterozygosity rate of the G. elegans genome was approximately 0.38%.

Genome assembly

Genome assembly was performed on full ONT long reads using Canu v1.7.1 and WTDBG v1.2.8. Because of a high error rate of Nanopore reads, we first corrected reads by the error correction module of Canu (canu -nanopore-raw -correct -fast genome size = 300 m). Then, the corrected reads independently assembled with WTDBG (wtdbg-1.2.8 --tidy-reads 8000 -k 0 -p 17 -S 2 --rescue-low-cov-edges;wtdbg-cns -k 13 -c 3). Finally, the preliminary genome assembly was approximately 331.8 Mb in size with a contig N50 size of 10.14 Mb (Supporting Information Table S5). Nanopolish calibration uses the Burrow-Wheeler Aligner (BWA, v0.7.12-r1039) default parameter to compare the quality-controlled Nanopore data to the assembled genome. The second-generation data are then compared to the Nanopolish-corrected genome using the BWA default parameter, and the Pilon iteration is used to correct it two times. The ultimate version of genome assembly was approximately 335.13 Mb in size with a contig N50 size of 10.23 Mb (Supporting Information Table S6). A guanine-cytosine (GC) depth analysis was conducted to assess the potential contamination during sequencing and the coverage of the assembly, revealing that the genome had an average GC content of 37% and a unimodal GC content distribution (Supporting Information Fig. S2). The GC depth as well as the sequencing depth of the genome assembly suggested that there was no contamination from other species (Supporting Information Fig. S3). G. elegans genome were performed with mitochondrial database in NCBI, the results showed that the coverage of some sequences was nearly 1, but the identity was low (Supporting Information Table S7). As mitochondrion was cyclic, the sequences might be short after the process of DNA extraction and library construction. When we used the long reads to assemble, the very short sequences were filtered which might include the mitochondrial sequence. So nearly all the assembled genome sequences were nuclear genome sequences.

Chromosome assembly using Hi-C data

Hi-C technology enables the generation of genome-wide 3D proximity maps and is an efficient and low-cost strategy for sequences cluster, ordered, and orientation for pseudomolecule construction. This technology has been successfully applied in recent complex genome projects, including goat, Tartary buckwheat, wild emmer, and barely. To generate a chromosomal-level assembly of the G. elegans genome, Hi-C fragment libraries were constructed. The Hi-C library was prepared followed by a procedure with an improved modification. In brief, freshly harvested leaves were cut into 2 cm pieces and vacuum infiltrated in nuclei isolation buffer supplemented with 2% formaldehyde. Crosslinking was stopped by adding glycine and additional vacuum infiltration. Fixed tissue was frozen in liquid nitrogen and grounded to powder before re-suspending in nuclei isolation buffer to obtain a suspension of nuclei. The purified nuclei were digested with 100 units of HindIII and marked by incubating with biotin-14-dCTP. Biotin-14-dCTP from non-ligated DNA ends was removed owing to the exonuclease activity of T4 DNA polymerase. The ligated DNA was sheared into 300–600 bp fragments, and then was blunt-end repaired and A-tailed, followed by purification through biotin-streptavidin-mediated pull down. Finally, the Hi-C libraries were quantified and sequenced using the Illumina Hiseq platform (Illumina, San Diego, CA, USA). In total, 370 million paired-end reads were generated from the libraries. Then, quality controlling of Hi-C raw data were performed using Hi-C-Pro (v2.8.0) as former research. Firstly, low-quality sequences (quality scores < 20), adaptor sequences and sequences shorter than 30 bp were filtered out using fastp v0.12.6 (fastp, RRID:SCR_016962), and then the clean paired-end reads were mapped to the draft assembled sequence using bowtie2 (v2.3.2) (bowtie2, RRID:SCR_005476) to get the unique mapped paired-end reads. As a result, 107 million uniquely mapped pair-end reads were generated, of which 76.28% were valid interaction pairs. Combined with the valid Hi-C data, we subsequently used the LACHESIS (ligating adjacent chromatin enables scaffolding in situ) de novo assembly pipeline to produce chromosome-level scaffolds. As shown in Fig. 2, the assembled sequence was anchored onto the 8 pseudo-chromosomes with lengths ranging from 36.08 to 52.33 Mb. The total length of pseudo-chromosomes accounted for 99.2% of the genome sequences, with scaffold N50 values of 40.47 Mb (Supporting Information Table S8).

Figure 2

Genome-wide Hi-C map of Gelsemium elegans. Interaction frequency distribution of Hi-C links among chromosomes shows in color key of heatmap ranging from light yellow to dark red indicated the frequency of Hi-C interaction links from low to high (0–10).

Evaluation of the completeness of the genome assembly gene space

To evaluate the coverage of the assembly, we randomly selected the RNAseq reads aligned against the G. elegans genome assembly using HISAT2 (hierarchical indexing for spliced alignment of transcripts2) with default parameters. The percentage of aligned reads ranged from 91.57% to 92.10% (Table S2). We then used Benchmarking Universal Single-Copy Orthologs (BUSCO, RRID:SCR 015008) to search the annotated genes in the assembly for the 1308 single-copy genes conserved among all embryophytes. About 92.9% of the complete BUSCOs were found in the assembly (Supporting Information Table S9). These results suggested that the genome assembly was complete and robust.

Genome annotation

The repeat sequences in the genome consisted of simple sequence repeats (SSRs), moderately repetitive sequences, and highly repetitive sequences. The microsatellite identification tool (MISA) was used to search for SSR motifs in the G. elegans genome, with default parameters. A total of SSRs were identified in this way: 134,047, 29,668, 9336, 1557, 409, and 524 mono-, di-, tri-, tetra-, penta-, and hexa-nucleotide repeats, respectively (Supporting Information Table S10).

To identify known transposable elements (TEs) in the G. elegans genome, RepeatMasker (RRID:SCR 012954) was used to screen the assembled genome against the Repbase (v22.11) and Mips-REdat libraries. In addition, de novo evolved annotation was performed using RepeatModeler v1.0.11 (RRID:SCR 015027). The combined results of the homology-based and de novo predictions indicated that repeated sequences account for 43.16% of the G. elegans genome assembly (Supporting Information Table S11), with long terminal repeats accounting for the greatest proportion of 23.9% (Supporting Information Table S12). The de novo and repbase RepeatMasker analysis of the G. elegans genome assembly are shown in Supporting Information Fig. S4.

Homology-based ncRNA annotation was performed by mapping plant rRNA, miRNA, and snRNA genes from the Rfam database (release 13.0) to the G. elegans genome using BLASTN (E-value ≤ 1 × 10−5). tRNAscan-SE v1.3.1 (tRNAscan-SE, RRID:SCR 010835) was used (with default parameters for eukaryotes) for tRNA annotation. RNAmmer v1.2 was used to predict rRNAs and their subunits. These analyses identified 208 miRNAs, 531tRNAs, 279rRNAs, and 1257 snRNAs (Supporting Information Table S13).

The homology-based, de novo based, and RNA sequences-based gene prediction methods were used to annotate protein coding genes. For homology-based predictions, protein sequences from 6 species (Arabidopsis thaliana, Calotropis gigantea, Camellia sinensis, Nicotiana tabacum, Olea europaea and Oryza sativa) (Supporting Information Table S14) were mapped onto the G. elegans genome; the aligned sequences and the corresponding query proteins were then filtered and passed to GeneWise v2.4.1 (GeneWise, RRID:SCR 015054) to search for accurately spliced alignments. For the de novo predictions, we first randomly selected 1000 full-length genes from the homology-based predictions to train model parameters for Augustus v3.0 (RRID:SCR 008417), GeneID v1.4.4, GlimmerHMM (RRID:SCR 002654), and SNAP. Augustus v3.0, GeneID v1.4.4, GlimmerHMM, and SNAP were then used to predict genes based on the training set. Further, G. elegans RNA-seq data and Iso-seq data were used for gene prediction by PASA (v2.0.2, RRID:SCR 014656). Finally, EVidenceModeler v1.1.1 was used to integrate the predicted genes and generate a consensus gene set (Table S14). Genes with TEs were discarded using the TransposonPSI package. Low quality genes consisting of fewer than 50 amino acids and/or exhibiting premature termination were also removed from the gene set, yielding a final set of 26,768 genes. The final sets average transcript length, average CDS length, average exon number per gene, average exon length and average intron length were 3961.71 bp, 1088.1 bp, 4.98, 218.63 bp and 722.54 bp, respectively (Supporting Information Table S15 and Fig. S5).

The annotations of the predicted genes of G. elegans were screened for homology against the Uniprot (release 2017/10) and KEGG (release 84.0) databases using Blastall and KAAS. Then, the InterProScan (release 5.2–45.0) package was used to annotate the predicted genes using the InterPro (5.21–60.0) database. In total, 22,636 of the total 26,768 genes (84.56%) were annotated with potential functions (Supporting Information Table S16).

After all the above prediction we used Benchmarking Universal Single-Copy Orthologs (BUSCO, RRID:SCR 015008) again to search the predicted genes in the assembly for the 1375 single-copy genes conserved among all embryophytes. About 95.5% of the complete BUSCOs were found in the assembly (Supporting Information Table S17). These results suggested that the genome prediction was complete and robust.

Phylogenetic tree construction and divergence time estimation

To investigate the evolutionary position of G. elegans, we compared its genome to the genome sequences of 8 other plants, which included 3 plants in special order or can produce alkaloids (C. gigantea, C. sinensis, and Macleaya cordata), 3 plants from different orders in the same Eudicots clade (A. thaliana, Brassica rapa and Vitis vinifera), and 2 monocotyledons (O. sativa and Oropetium thomaeum) as an outgroup. We used the OrthoMCL (v2.0.9) pipeline (OrthoMCL DB: Ortholog Groups of ProteinSequences, RRID:SCR 007839) (BLASTP E-value ≤ 1 × 10−5) to identify potentially orthologous gene families within these genomes.

Gene family clustering identified 13,792 gene families containing 20,755 genes in G. elegans (Fig. 3). Of these, 903 gene families were unique to G. elegans (Supporting Information Table S18).

Figure 3

Venn diagram of shared gene families between Gelsemium elegans and 8 other plants. Each number represents a gene family number.

Phylogenetic analysis was performed using 2989 single-copy orthologous genes from common gene families found by OrthoMCL (Supporting Information Fig. S6). We codon-aligned each gene family using MUSCLE (MUSCLE, RRID:SCR 011812) and curated the alignments with Gblocks v0.91b. Phylogeny analysis was performed using RAxML (RAxML, RRID:SCR 006086) v8.2.11 with the GTRGAMMA model and 100 bootstrap replicates.

We then used MCMCTREE as implemented in PAML v4.9e (PAML, RRID:SCR 014932) to estimate the divergence times of G. elegans from the other plants. The parameter settings of MCMCTREE were as follows: clock = 2, RootAge <1.93, model = 7, BDparas = 110, kappa gamma = 62, alpha gamma = 11, rgene gamma = 23.18, and sigma2 gamma = 14.5. In addition, the divergence times of O. sativa (148–173 Mya), V. vinifera (110–124 Mya), and A. thaliana (53–82 Mya) were used for fossil calibration.

The phylogenetic analysis showed that G. elegans is more closely related to C. gigantea than to C. sinensis (Supporting Information Fig. S7), which supports the well-established hypothesis of a close relationship between G. elegans and C. gigantea65, 66. The estimated divergence time of G. elegans and C. sinensis was 97.45 Mya, while that of G. elegans and C. gigantea was about 50.69 Mya (Fig. 4).

Figure 4

Inferred phylogenetic tree across 9 plant species. The estimated divergence time (Mya) is shown at each node.

Genes under positive selection

Studies on the crude and purified alkaloids of G. elegans have demonstrated that this species possesses anti-inflammatory, immunomodulating, analgesic, anxiolytic, anti-tumor8, 9, and neuropathic pain-relieving properties. The ratio of nonsynonymous substitution rate (Ka) and synonymous substitution rate (Ks) of protein coding genes can be used to identify genes that show signatures of natural selection. We calculated average Ka/Ks values and conducted the branch-site likelihood ratio test using Codeml implemented in the PAML package to identify positively selected genes in the G. elegans lineage. The parameter settings of Codeml were as follows: Model A: model = 2, NS sites = 2, fix_omega = 0; Model A1: model = 2, NSsites = 2, fix_omega = 1, omega = 1.

These genes might contribute to the secondary metabolites of adaption to unfavorable environments. 94 Genes with signatures of positive selection were identified (P ≤ 0.05), of which 77 genes could be annotated with potential functions in the Swissprot database (Supporting Information Table S19). One gene is homologous required for transport of secretory proteins from the Golgi complex, which catalyzes the transfer of phosphatidylinositol and phosphatidylcholine between membranes in vitro. This gene could potentially contribute to the adaption of G. elegans to the secondary metabolites of environment. While literature reports are rare, other identified genes might also be associated with the adaption of G. elegans. It should be noted that this is just a preliminary analysis of the functions of these genes, and further studies would be needed to clarify their roles.

Whole-genome duplication and gene family expansion analysis

We used 4-fold synonymous third-codon transversion (4DTv) and Ks estimation to detect whole genome duplication (WGD) events in the G. elegans genome. To this end, paralogous sequences of G. elegans, A. thaliana, Glycine max, O. europaea, and V. vinifera were identified with OrthoMCL. Then, protein sequences for each of these plants were aligned against each other with Blastp (using an E-value threshold of ≤1 × 10−5) to identify conserved paralogs in each species. Finally, potential WGD events in each genome were evaluated based on their 4DTv and Ks distribution. The WGD analysis suggested that G. max and O. europaea may have experienced modern WGD events, and G. max have gone through 2 times of whole genome duplications, while the G. elegans has no modern WGD event, and only experience done ancient whole genome duplications (Fig. 5).

Figure 5

Whole-genome duplication (WGD) events of 5 plants (Gelsemium elegans, Arabidopsis thaliana, Glycine max, Olea europaea and Vitis vinifera) inferred by 4-fold synonymous third-codon transversion (4DTv) estimations. (A) 4DTv; (B) Ks.

The OrthoMCL gene family analysis results were analyzed further by using Computational Analysis of Gene Family Evolution v3.0 to detect expanded gene families. This approach revealed 509 expanded gene families and 1013 contracted gene families in the G. elegans lineage (Supporting Information Fig. S8).

Results and discussion

The current study shows that blueprint of an organism is encoded in its genome and genome mining has become a powerful strategy for botanical studying. To investigate the evolutionary history of the indole alkaloid gene cluster, we performed 2 rounds of synteny analysis with either the “all BLASTp” result as input of blocks with distant homology or the default “top 5 BLASTp” result for blocks with close homology. The top ranked syntenic block for the indole alkaloid (e.g. koumine and gelsemine) pathway genes is found (Supporting Information Table S20). Key steps in monoterpene indole alkaloid (MIA) biosynthesis are shown in Scheme 1, catalyzed by the enzymes tryptophan decarboxylase (TDC), strictosidine synthase (STR), strictosidine glucosidase (SDG), geissoschizine dehydrogenase (GSD) and sarpagan bridge enzyme (SBE). The MIAs comprise approximately 3000 compounds which have different chemical scaffolds. This enormous chemical complexity stand in contrast with only 3 sequenced MIA producers, namely Catharanthus roseus, Rhazya stricta (Apocynaceae, order Gentianales) and Camptotheca acuminate (Nyssaceae, order Cornales). Nevertheless, only a small number of the genes primarily from C. roseus related to the enormous diversity of MIAs are known so far, as a result, the genomic context of MIA biosynthesis is largely unknown and has only been systematically investigated in C. roseus. Given the vast chemical diversity of MIAs, we wondered whether these gene clusters would be conserved in MIA producing plants with different chemical profiles, and whether they might potentially be useful for accelerating biosynthetic gene discovery. As a case study, we selected the MIA producer G. elegans, which produces a wide variety of indole and oxindole alkaloids.

Scheme 1

Key steps in monoterpene indole alkaloid (MIA) biosynthesis, catalyzed by the enzymes tryptophan decarboxylase (TDC), strictosidine synthase (STR), strictosidine glucosidase (SDG), geissoschizine dehydrogenase (GSD) and sarpagan bridge enzyme (SBE). The pathway diverges after strictosidine aglycone and leads to very different alkaloids in the plants Gelsemium elegans. Several representative alkaloids for Gelsemium elegans are shown.

Conclusions

This paper reports the sequencing, assembly, and annotation of the G. elegans genome along with details of its evolutionary history and alkaloids metabolism. Thus, we generate a significantly improved genome sequence than another Gelsemium family plant G. sempervirens that is 244 Mb with an N50 scaffold size of 411,072 bp. The genomic data generated in this work will be a valuable resource for further genetic improvement and effective use of the G. elegans.

Availability of supporting data

The raw data from our genome project was deposited in the SRA (Sequence Read Archive) database of National Center for Biotechnology Information with Bioproject ID PRJNA505365 (Biosample ID from SAMN11089884 to SAMN11089892). We also upload the Hi-C results with the same with Bioproject ID PRJNA505365 (Biosample ID from SAMN12083642 to SAMN12083645). Versions and main parameters of the software used in this study are provided in Supporting Information Table S21.