Background & Summary

Davidia involucrata Baill., also called dove tree or handkerchief tree, is the sole species in the genus Davidia (Davidiaceae1 or Nyssaceae2) and is listed as a “first-grade” nationally protected plant in China3,4. It is a Tertiary paleotropical relic plant species that is rare in China and usually regarded as a “botanic living fossil”5. Its distribution is limited to the subtropical mountains of central to southwestern China; natural populations are often found in deciduous or evergreen broad-leaf forests at elevations of 1100–2600 m6. D. involucrata is not only an endangered and rare relic species but also famous as an ornamental plant by virtue of the pair of large white bracts that surround the small flowers and create the appearance of doves perching among its branches, giving the tree its common name7. The most unusual characteristics of D. involucrata are its floral organs, of which the inflorescences contain either a mixture of hermaphrodite and many male flowers or entirely male flowers; the flowers are without petals but have large, unequal, paired paper-like bracts instead. These intriguing bracts originally appear small and green, resembling leaves, but they increase in size and turn white as the flowers mature, and then, finally, turn to brown and yellow before being shed8,9. Bracts are special organs that appear in the reproductive development of plants, and D. involucrata is undoubtedly an ideal subject for the study of bracts and the developmental mechanisms of specific flower organs.

To date, studies of D. involucrata have mainly focused on the macroscopic and phenotypic levels, such as taxonomy, morphology, physiology, ecology, reproduction and so on9,10,11,12,13. However, research at the molecular level continues to progress slowly, which could result from the distribution of this species, which is intermittently scattered throughout southwest China, and its growth characteristics; D. involucrata is slow-growing, occurs in harsh growth conditions and has a low survival rate. Therefore, deeper research at the molecular level is critical for this endangered plant.

Studies of the molecular mechanisms underlying the growth and development of D. involucrata as well as its unique bracts and floral organs could further reveal the evolution of floral development, which remains unclear due to a lack of high-throughput data. To date, many molecular studies have focused on exploring and analysing traditional markers, such as simple sequence repeats (SSRs)14,15,16, random amplified polymorphic DNA (RAPD) markers17, microsatellites18 and chloroplast genes19, or utilized functional genes or factors, such as a MYB transcription factor from D. involucrata (DiMYB1)20, a cold-induced gene (DiRCI)21, a clathrin adaptor complex gene (DiCAC)22 and so on. At the same time, the transcriptome of the seed23 and the chloroplast genome24 have been reported. However, transcriptomic and sRNA data from some crucial tissues, for example, the floral organs and, specifically, the bracts, have not been examined in depth23,24,25.

Here, we established a complete gene set from multiple tissues of D. involucrata by means of next generation sequencing technology. This gene set would be useful for further studies. For example, it could be used as a reference for gene characterization, such as expression analysis, gene cloning and phylogenetic analysis, and it could also be used for gene model annotation once the genome is sequenced. In addition, the transcriptome annotations could be used to explore the crucial genes in flower development and bract development, which are of great significance to reveal the evolution of the floral organs in angiosperms. Furthermore, the small RNAs of each sample were sequenced and annotated in parallel to provide more information about sRNA-related regulatory mechanisms during floral organ development. The whole transcriptomes established herein lay a foundation for plant molecular marker-assisted breeding, evolutionary and developmental analysis, and even plant protection in the future.

Methods

Sample collection

The floral organ samples were collected from 6 blooming individuals at two developmental stages (3 individuals for each stage) in a wild population in the county town of Yingjing, Yaan, Sichuan Province, in April 2014 (Fig. 1a,b). Floral organ growth stages were defined based on the colour changes of the bracts. One stage was young inflorescences with small, green bracts (expanded to approximately <8 cm in length) resembling leaves, and the dark purple anthers were immature (Fig. 1a). The other stage was mature inflorescences with large white bracts (expanded to approximately >15 cm in length), and the anthers were mature (Fig. 1b). Each tissue sample collected from each plant was approximately thumbnail-sized and was immediately stored in RNA Fixer solution (Bioteke, China) for further use.

Fig. 1
figure 1

Schematic pipeline illustrating analysis of the whole transcriptome dataset. (a) The young inflorescences of D. involucrata. (b) The mature inflorescences of D. involucrata. (c) After library construction, mRNAs were sequenced with a PE90 strategy, and sRNAs were sequenced with an SE50 strategy. After low-quality read removal, the transcriptome of each sample was assembled using SOAPdenovo-trans, KGF and GapCloser. After clustering using TGICL, the final assemblies were obtained, and addition evaluation and annotation were performed.

Full size image

The floral organ samples (Table 1) were named YB (young bracts: a pool of similar samples collected from the growing green bracts of three individuals), LB (mature bracts: a pool of similar samples collected from the mature white bracts of three individuals), LX (mature stamens: a pool of similar samples containing the complete, mature stamens from three individuals), LC (mature pistils: a pool of similar samples containing the complete, mature pistils from three individuals), YR (young mixed samples: a pool containing flowers at the early blooming stage including equal amounts of complete stamens and complete pistils from three individuals), ZR (mixed samples: a pool of all collected floral organ samples with added leaves, shoots, and roots at two periods from all six individuals).

Table 1 Summary of transcriptome reads.
Full size table

Total RNA extraction

Total RNA extraction was performed using the TRIzol Reagent (Invitrogen, Carlsbad, CA) according to the manufacturer’s instructions, followed by quality assessment on an Agilent 2100.

Transcriptome sequencing and filtering

The mRNA was extracted from the total RNA with oligo (dT)-attached magnetic beads, and a cDNA library with an insert size of 250 bp was constructed using the TruSeq RNA Sample Preparation Kit according to the manufacturer’s protocol (Illumina Inc., San Diego, CA). Library generation yielding 2 × 90 bp paired-end reads and sequencing (Illumina HiSeq 2000) were performed at BGI-Shenzhen. Low-quality reads matching one or more of the following criteria were filtered out using SOAPnuke (v1.5.6)26: reads containing adaptor contamination; reads including more than 5% of the unknown base “N”; reads including more than 20% bases with quality values lower than 10.

sRNA sequencing and filtering

RNA segments of different sizes from 18–30 nt were separated from the total RNA by 15% denaturing PAGE and selected for small RNA library construction using the Illumina TruSeq Small RNA Sample Preparation Kit according to the manufacturer’s protocols (Illumina Inc.). Library generation yielding 50 bp single-end reads and sequencing (Illumina HiSeq. 2000) were also performed at BGI-Shenzhen. The contaminant tags and low-quality tags were removed using SOAPnuke (v1.5.6)26: (1) tags with 5′ primer contaminants, oversized insertion, poly-A; (2) tags without 3′ primer or insert tags; (3) tags shorter than 18 nt. After the removal of low-quality reads, clean reads (Table 2) were retained and used in subsequent analyses (Fig. 1). The quality of all clean reads was assessed with FastQC (http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc) and MultiQC (v1.8)27.

Table 2 Summary of sRNA reads.
Full size table

Transcriptome assembly

Transcriptome assembly was completed using a combined assembly strategy. Briefly, the RNA-Seq reads were assembled using SOAPdenovo-Trans (version 1.01)28 with the following settings: “-K 31-i 20 -e 2 -M 3-L 100”. The gaps were filled using KGF (v1.19) and the GapCloser tool (v1.12)29. All assemblies were merged into one large dataset using TGICL (v2.0.6)30 with the parameters “-l 40 -c 10 -v 25 -O ‘-repeat_stringency 0.95 -minmatch 35 -minscore 35’”.

To remove unreliably assembled transcripts, clean reads were aligned to the assembly using Bowtie231 with the parameters “–sensitive–score-min L,0,-0.1 -I 1 -X 1000–mp 1,1–np 1–no-mixed–no-discordant”, and the fragments per kilobase of exon model per million reads mapped (FPKM) values of the transcripts were calculated using the tool rsem-calculate-expression in RSEM (v1.2.21)32. The sequences with FPKM values of zero were removed from the assembly. Then, the lowest common ancestor (LCA) algorithm was applied to filter out contaminants. First, the assembly was used to search against the NCBI non-redundant protein database (Nr) using Diamond with e-value < 1e-5. The tool Blast2lca in MEGAN (v6.15.2)33 with the parameter “–minScore 75” was used to apply the LCA alignment and produce taxonomic classifications. Non-Viridiplantae sequences were removed according to their classifications. All unreliably unigenes of each individual were also filtered. Last, we obtained the final transcriptome assembly (Table 3).

Table 3 Summary of unigenes and the final assembly.
Full size table

Transcriptome Annotation

Functional annotation

Functional annotations were performed using a sequence-based search method (Table 4). First, the final assembly was annotated by searching against the NCBI non-redundant nucleotide (Nt) database34 using BLASTn (v2.4.0)35 with e-value < 1e-5, Viridiplantae-related non-redundant proteins in the Nr database, Swiss-Prot in UniProtKB36, the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database37 and clusters of euKaryotic Orthologous Groups (KOG)38 using BLASTx (v2.4.0) with e-value < 1e-5, and the InterPro database39 using InterProScan (v5)40 with the default parameters36. Gene Ontology (GO)41 annotation was performed using Blast2GO (v2.5.0)42 based on the Nr annotation results.

Table 4 Summary of functional annotation.
Full size table

Identification of coding gene set

First, the open reading frames (ORFs) of each sequence were predicted using TransDecoder (v3.0.1) as implemented in Trinity43. Then, the ORFs were searched against Swiss-Prot36 using Diamond Blastp (v0.8.31)44, and the output file was searched against Pfam45 using Hmmscan (v3)46. Finally, CDSs were predicted using TransDecoder with the results of the previous step (Table 5).

Table 5 Summary of CDS prediction and coding gene set.
Full size table

Identification of high-quality annotated coding gene set

To identify high-quality coding genes, we used the following pipeline. (1) All transcripts had matches with protein databases, and the best hit for each transcript (in the following order: Nr, Swiss-Prot, KEGG, KOG and InterPro) was extracted. (2) The annotated coding region and functional annotation for each transcript were selected from the best hit. (3) The transcript with the longest CDS was chosen if its annotated coding region could also be identified by TransDecoder. (4) After filtering out sequences with more than 5 Ns in the last 10 bases or more than 15 Ns in the last 50 bases, the eligible CDS of each transcript was extracted as the final coding gene. In total, 51,247 annotated coding genes with a minimal length of 297 nt and a maximal length of 12,852 nt were identified as high-quality coding genes (Table 5).

Transcription factor family classification

The open reading frame of each sequence was classified using getorf (EMBOSS:6.5.7.0)47, and the transcriptome factor (TF) family was identified based on plant TF domains in the plant TF database (PlnTFDB)48 using Hmmsearch46. In total, 3,064 genes were arranged in 58 TF families, and the top three annotated families were the Myb DNA-binding domain families MYB and MYB-related and the basic helix-loop-helix (bHLH) family (Fig. 2a).

Fig. 2
figure 2

Transcription factor family and simple sequence repeat statistics. (a) Transcription factor family classifications of the assembled sequences. (b) SSR statistics.

Full size image

Simple sequence repeat detection

SSRs in the transcriptome sequences were detected using MISA (v1.043)49 with the definition parameters “1-12,2-6,3-5,4-5,5-4,6-4” and an interruption parameter of “100”. The total number of identified SSRs was 44,327, and the number of SSR-containing sequences was 20,739 (Fig. 2b).

sRNA Annotation

Small RNA annotation

Because of the lack of a D. involucrata special sRNA reference dataset, we used the microRNA datasets of Arabidopsis lyrata from miRBase (release 22)50 and RNA datasets from the Rfam database (v12.1)51 as the known small RNA reference databases. Clean reads were mapped to the miRBase database using Bowtie2 (v2.2.5)31 with “–sensitive -L 16” and the Rfam database using cmsearch in INFERNAL (v1.1.2)52 with “–noali” to obtain the annotations. Aligned tags were annotated after filtering reads with more than one mismatch out of each alignment. A total of 96 known miRNAs were annotated in this study.

Novel miRNA prediction

Because no more than 10% of the small RNA tags were annotated, miRA software53 was applied to predict novel miRNAs. First, the clean reads of each sample were mapped to the “genome” reference (the merged transcriptome) using Bowtie2 with the parameters “-L 16–rdg 1,10–rfg 1,10”. After removing the unreliable transcriptome sequences from the alignments, 45.28%, 44.08%, 43.62%, 44.21%, 38.90% and 54.74% tags could be mapped to the transcriptome for ZR, YR, YB, LB, LC and LX, respectively. After filtering out the annotated sRNA reads, the aligned read tags from all samples were combined into one FASTA file. Then, the combined FASTA file was mapped to the genome reference using Bowtie2 with “-f -l 16”, and the alignment file was filtered using SAMtools (v1.3.1)54 with “view -hS -F 4”. Finally, miRA (v1.2.0) with the default parameters was used to predict novel miRNAs. In total, 112 novel miRNAs were identified. The annotation and prediction results are summarized in Table 6 (Table 6, Fig. 3).

Table 6 sRNA annotation and prediction results.
Full size table
Fig. 3
figure 3

Length distributions of sRNA samples and sRNA annotation distribution of sample ZR. (a) Length distributions of sRNA. The peak of each sample was located at 24 nt. (b) sRNA annotation distribution of sample ZR.

Full size image

miRNA target prediction

psRobot (v1.2)55 and TargetFinder (v1.0)56 with the default parameters were used to predict miRNA targets. The intersection of the target genes was extracted as the final prediction target result. The final intersection results included 172 miRNAs and 1,737 target genes.

Data Records

The final transcriptome and related data are published under the International Nucleotide Sequence Database Collaboration BioProject PRJNA513477 (https://identifiers.org/ncbi/bioproject:PRJNA513477) and CNGB Nucleotide Sequence Archive (CNSA) project CNP0000260 (https://db.cngb.org/search/project/CNP0000260). The read files were deposited in the NCBI Sequence Read Archive57. The final high-quality coding gene set and the transcriptome assemblies were deposited in NCBI Transcriptome Shotgun Assembly58 and CNSA. Annotation data set was uploaded to figshare59.

Technical Validation

Quality control and data statistics

To control the sequencing quality, the counts of clean reads, total bases, quality scores and GC content were calculated for each sample using FastQC and MultiQC (Tables 1, 2 and Fig. 4). The mean read counts per quality scores and the mean quality scores in each base position were higher than 30. The length distribution of clean small RNA tags showed that the peak of each sample was located at 24 nt (Fig. 3).

Fig. 4
figure 4

Quality assessment. Read count distributions by mean sequence quality of transcriptome reads (a) and sRNA reads (c). Mean quality score distributions of transcriptome reads (b) and sRNA reads (d).

Full size image

Assembly quality control

Different tissues from different trees were collected in this study to comprehensively cover the D. involucrata transcriptome. Because reads were aligned to the merged assembly, we summarized the mapping percentages in Table 1 and then calculated the FPKM value for each sequence. The read mapping results showed that more than 90% of the reads were mapped to the transcriptome, and the sequences with FPKM values of 0 were removed. At the same time, the contaminant sequences were removed according to the LCA-based taxonomy classification method. The N50 and GC content for each assembly were also calculated using NGS QC Toolkit (v2.3.3)60 (Table 3).

We then employed BUSCO (v3)61 to evaluate the completeness of the final assembly using the 1,440 Embryophyta expected genes database (version 2). This analysis showed (Table 7) that 1,285 (89.24%) and 50 (3.47%) of the 1,440 expected Vertebrata genes were identified as complete and fragmented, respectively, while 105 (7.29%) genes were considered missing in the final assembly. BUSCO was also used to evaluate the completeness of the final gene set, and 87.99% of the 1,440 expected genes were identified (Table 7).

Table 7 Summary of BUSCO and PLAZA results.
Full size table

To assess the completeness of the core gene families (CoreGFs) within the green plant lineage, the CoreGF score was calculated using PLAZA (v2.5)62. First, the assembled transcriptome sequences were searched against the CoreGF gene set using BLASTx, and the protein sequences of high-quality CDSs were searched against the CoreGFs using BLASTp (V2.4.0) with e-value < 1e-5. Then, scores were calculated using the script coreGF_plaza2.5_geneset.py in PLAZA based on the BLAST hits. The CoreGF weighted score was 98.85%, and only 33 out of 2,928 CoreGFs were missing from the final assembly, while the score was 97.13%, and only 85 CoreGFs were missing in the high-quality coding gene set.

Contamination screening

The transcriptome sequences were subjected to contamination screening as mentioned in the methods, and the contaminant sequences were removed. The results of the contamination analysis also showed that (1) the main contaminants came from fungi or Arthropoda; (2) the most abundant fungi were Paraphaeosphaeria sporulosa in the Ascomycota; and (3) the transcriptome showed the highest similarity with Prunus avium, followed by Theobroma cacao and Nicotiana attenuata, in Viridiplantae (Fig. 5a).

Fig. 5
figure 5

Functional annotation. (a) Venn diagram of annotations based on the databases NR, KOG, KEGG, Swiss-Prot and InterPro. (b) Species distribution of annotated database NR sequences. (c) KEGG pathway annotations. Bars represent the numbers of unigenes clustered into KEGG Orthology (KO) hierarchies.

Full size image

sRNAs were assembled using velvet (v1.2.10)63 with the parameters “velveth Assemout 15-short -fastq.gz” and “velvetg Assemout”. Assembled contigs with lengths greater than 100 bp were extracted to detect potential symbionts using the virus detection pipeline published previously64. In brief, a total of 115 queries with lengths greater than 100 bp were searched against the NCBI virus database using BLASTn with e-value < 1e-5. All matched sequences were searched against the Nt database using BLASTn with e-value < 1e-5 to identify false positives. After detection, none of the known viruses were detected in the sRNA data. In addition, sRNA reads mapped to the contaminant transcriptome sequences were removed from further analysis.

Annotation quality control

Statistical results of the functional annotation are summarized in Table 4 (Fig. 5b). A total of 65,534 (62.75%) unigenes were annotated. The distribution of KEGG pathway annotations is shown in Fig. 5c.

Statistical results for the predicted CDSs are summarized in Table 5. A total of 58,561 coding regions were detected by TransDecoder. The total, maximum and minimum lengths of the predicted CDSs were 60,633,393, 12,852 and 297. N50 and GC content were also calculated. Statistical results for the final high-quality coding genes are also summarized in Table 5.