Abstract

Background

Understanding the molecular mechanisms underlying the adaptation of woody plants to local environmental conditions remains a big challenge in biology due to their long and perennial life history. However, woody plants represent a large proportion of biodiversity on the earth and harbor many different phenological traits that herbaceous plants do not feature (https://www.worldwildlife.org). One such trait is the continuous flowering behavior of rose, an important crop of high importance in human society. Guaranteeing a constant supply of raw materials for cut flowers and related products, continuous flowering becomes one of the most important biological and economical traits for roses [1]. Therefore, the genetic control and related gene-regulatory-network for continuous flowering regulation attracts efforts for many years not only from scientists but also from breeders [1, 2]. The number of QTL regulating continuous flowering remains disputable [1–3]. RoKSN, a homolog of Arabidopsis TFL1-like gene in roses, is the only known gene responsible for continuous flowering [4, 5].

Domestication of cultivated roses mainly involves hybridization among more than a dozen species [2, 6–8]. Frequent inter-species crossing/backcrossing and polyploidization of roses has made the classification of roses very difficult [9–12]. A set of high-quality and well-characterized genomic tools/resources are necessary for understanding the biology and domestication of modern roses that encompass more than 30,000 cultivars [13]. Recently, several genetic mapping populations have been developed (see reviews [1, 2]) and determination of rose genomes is carried out [14] with the sequence released very recently for a doubled-haploid of Rosa chinensis ‘Old Blush’ [15]. However, due to high-level heterozygosity caused very likely by inter-species crossing and polyploidization, achieving an accurate and complete rose genome seems not so easy. Alternatively, a comprehensive gene expression atlas can be constructed with multiple tissues in different species.

The first sets of gene expression atlas were constructed using microarrays containing about 350 (tetraploid R. hybrida) [16] and later with about 4800 selected ESTs (R. chinensis, R. wichuriana, and R. hybrida) [17]. A more comprehensive database containing about 80,714 transcript clusters for R. chinensis ‘Old Blush’ was constructed from 13 tissues/organs at different developmental stages or under different abiotic and biotic stresses with Illumina and 454 sequencing platforms [18]. Several recent studies have also been done with various Rosa species for different purposes [19–22] even for a single-nucleotide-polymorphism (SNP) array [23, 24]. Though all these studies have promoted significantly the understanding of the rose biology, the quality of these transcriptomes is normally relatively poor with low N50 value, poor completeness, and short average length. Very recently the publishing of genome sequences of a doubled-haploid for R. chinensis ‘Old Blush’ set a milestone for rose research [15, 25]. However, missing half genome information of roses, most of which often feature high heterozygosity due to frequent intra−/inter-species hybridization and polyploidization, might cause lacking of power in understanding the genetic bases of roses traits associated with recessive markers [3].

To identify the molecular components underpinning rose specialty and to find molecular markers for clarifying phylogenetic relationship of Rosaceae, we generated a set of high quality reference transcriptomes for two rose species, which featured at least six pair of contrast traits [3], by sequencing three tissues at different developmental stages and by integrating published datasets. We identified about 4447 transcripts conserved in Rosaceae plants, among which 405 were under significant selection pressure, and 164 transcripts present only in roses.

Results

Transcriptome sequencing and assembling

RNA samples from shoot materials and young leaves were profiled with high-throughput sequencing (Figs. 1 and ​and2).2). This resulted in 421.8 million and 427.5 million clean reads for Rosa chinensis ‘Old Blush’ (OB) and R. wichuriana ‘Basye’s Thornless’ (BT), respectively (Table 1). The raw sequence files have been uploaded to the National Center for Biotechnology Information Sequence Read Archive (http://www.ncbi.nlm.nih.gov/sra/) under the accession numbers SAMN07808857–07808870.

Open in a separate window

Fig. 1

Leaf (a) and shoot (b) materials used for RNA-seq in this study. For each panel, left for Rosa wichuriana ‘Basyes’ Thornless’ (BT), and right for R. chinensis ‘Old Blush’ (OB). Bars=1cm

Open in a separate window

Fig. 2

Working flow for assembling of reference trancriptomes and identification of Rosacaeae-common and Rosa-specific transcripts. Main steps were shown in boxes with key transcripts numbers given. Major tools used in these analysis were marked in blue. Dashed arrows and boxes indicated the data generated from this study could be explored in these applications

Table 1

Summary of sequencing strategies and sequences obtained

Species	Sample	Repetition	Reads number	Reads bases (nt)	Q20 (%)	Accession code
Rosa wichuriana ‘Basye’s Thornless’ (BT)	SAMa	Total	156,722,692	14,104,983,240	95.8	_
		Rep1	51,669,692	4,650,272,280		SAMN07808870
		Rep2	53,341,618	4,800,745,620		SAMN07808871
		Rep3	51,710,726	4,653,965,340		SAMN07808872
	leaf_novb	Total	131,025,452	18,923,248,232	95.8	_
		Rep1	42,294,066	6,100,424,928		SAMN07808867
		Rep2	42,725,012	6,170,189,794		SAMN07808868
		Rep3	46,006,374	6,652,633,510		SAMN07808869
	leaf_marb	Total	134,110,978	19,389,814,312	96.1	_
		Rep1	44,236,512	6,418,203,960		SAMN07808864
		Rep2	46,651,880	6,762,697,476		SAMN07808865
		Rep3	43,222,586	6,208,912,876		SAMN07808866
Rosa chinensis ‘Old Blush’ (OB)	SAMa	Total	159,834,774	14,385,129,660	95.3	_
		Rep1	52,443,532	4,643,438,236		SAMN07808861
		Rep2	53,284,740	4,800,745,624		SAMN07808862
		Rep3	54,106,502	4,940,945,800		SAMN07808863
	leaf_novb	Total	137,399,212	19,719,505,590	95.8	_
		Rep1	44,950,340	6,501,086,190		SAMN07808858
		Rep2	45,260,556	6,426,857,278		SAMN07808859
		Rep3	47,188,316	6,791,562,122		SAMN07808860
	leaf_marb	Total	130,341,518	18,669,073,728	96.7	_
		Rep1	42,299,982	6,021,996,256		SAMN07808855
		Rep2	43,001,278	6,170,274,204		SAMN07808856
		Rep3	45,040,258	6,476,803,268		SAMN07808857
	all other datac		550,108,308	63,356,156,640	95.6	_

Open in a separate window

Data are sum of three biological replications. ^aand ^b, samples sequenced via Illumina pair-end methods (PE100bp for ^a and PE150bp for ^b);^c, data from references (see Table ​Table22)

Next we generated the assemblies for BT (68,612 transcripts) and OB (74,975 transcripts) with transcripts N50 about 2099/1732bp (BT/OB) and average length about 1338/1170bp (BT/OB). This was much better than previously published data for these two species (Table 2, Additional file 1: Table S1, Additional file 2: Table S2; Fig. 3a) [18, 20, 26, 27] and other species/materials in Rosa [19, 21, 22, 28]. A further assembly with all available data for OB produced even better transcript N50 (2092bp) and average length (1359bp). Comparing the assembled OB transcripts with the genome assembly (v2.0_a1) that containing 45,469 coding and 4918 non-coding genes [15] revealed that almost 100% of these transcripts could be mapped. On the other hand, by looking for the genes included in the genome assembly in the assembled transcripts we identified about 97% of these genes in our OB assemblies. Mean GC content of all assemblies (44.2–46.4%) was comparable to that of published (45.8–46.5%) roses (Table ​(Table2).2). A BUSCO analysis revealed a high proportion of complete (C) and single copy (S) from 54.4 to 68.8%, and complete (C) and duplicated (D) from 24.5 to 28.2%. Fragmented (F) and missing (M) BUSCO items occupied about 5.6–17.4% (Fig. ​(Fig.3b;3b; Additional file 1: Table S1). The prediction of high number of transcripts could be correlated with high frequency of intra−/inter-species hybridization and polyploidization in roses. Taken together, these results suggest the high quality, completeness, and coverage of the assembled BT and OB transcriptomes.

Table 2

Statistics of final assemblies for this study and published data

Assembly components	Contig number	Transcript number	Transcript N50	GC content (%)	Total assembled bases	Average length (bp)	Data sources
BTa	86,642	68,612	2099	46.4	92M	1338	This study
OBa	99,456	81,389	2092	44.2	111M	1359	This study
OB	na.	80,714	na	na.	36M	444	Dubois et al. 2012 [18]
OBb	na.	68,565	na.	46.46	61M	887	Yan et al. 2016 [20]
OB	na	85,663	na	na.	70M	814	Guo et al. 2017 [26]
OBc	208,039	111,954	1997	45.8	231Mc	1111	Han et al. 2017 [27]
Core set 1 BT vs. OB; 90% identity	na.	10,773	2282	na.	20M	1863	This study
R. multiflora	78,676	61,864	1907	46.03	75M	1216	Zhang et al. 2016 [21]
R. jacq cv. gold medal	na.	80,226	na	na.	60M	743	Gao et al. 2016 [22]
R. roxburghii	na.	106,590	na	na.	37M	343	Yan et al. 2015 [28]
R. hybrida	93,947	na.	1589	na.	na.	na.	[22]
R. chinensis ‘pallida’	na.	89,614	Na.	na.	38M	428	Yan et al. 2014 [19]
Rosa transcriptomed	60,944	na.	314	na.	18M	302	Fei lab
Core set 2 All samples; 80% identity	na.	5959	2326	na.	13M	2161	This study

Open in a separate window

^aassembly based on data produced from this study; ^b, assembly based on data from this study and references Yan et al. [20] and Han et al. [27]; ^c, conceptual confusion in original text; ^d, data from Fei lab (http://bioinfo.bti.cornell.edu/cgi-bin/rose_454/index.cgi) with transcript N50 and average length recalculated

Open in a separate window

Fig. 3

The assembly of high quality transcriptomes for roses. a Length distribution in proportion of assembled unigenes for the two species, Rosa chinensis ‘Old Blush’ (OB, bars filled in grey color), and R. wichuriana ‘Basyes’ Thornless’ (BT, open bars). Bars filled in black color mark the length distribution of shared transcripts between the two species (coreset1; see below and main text). b BUSCO analysis shows the completeness of assemblies and coreset1. c Annotation results of the assembled unigenes and core-sets for Rosa. The coreset1 is between the two species, while coreset2 is for the unigenes shared among Rosa (see Fig. ​Fig.2)2) based on published and newly collected data from this study. For each category (Nr_plants, GO, Uniprot, Swissprot and COG databases), total unigene counts annotated in different databases besides the proportion (in brackets) are given. Shared and total unigenes annotated by all databases are also given. d Venn diagram shows the results of coreset1 identification. About 10,773 transcripts were identified at the 95% sequence identity level between the two species

Identification of the Rosaceae-common and Rosa-specific transcripts

One of the purposes of this research was to find transcripts/markers that could be potentially used for phylogenetic reconstruction of Rosaceae plants. We then compared the coreset2 transcripts with genes from M. domestica, R. occidentalis, P. avium, and Fragaria veseca, the closest relatives of Rosa (Fig. 4a). About 4447 transcripts were found to present in all five genus (Rosaceae-common; 74.6% of coreset2 transcripts), while only 164 transcripts were specific for Rosa plants (Rosa-specific; 2.75% of coreset2). Most of the Rosaceae-common transcripts were annotated (4228 or about 96%; Table ​Table22 and Additional file 5: Table S4) and could encode for proteins (4341; 97.6%). About 504 Rosaceae-common transcripts belong to the 1440 BUSCO single-copy genes (35%; Additional file 5: Table S4, Additional file 6: Table S5). GO enrichment assay revealed that most of these transcripts were involved in very important and basic function (Fig. ​(Fig.4b,4b, Additional file 5: Table S4). Phylogenetic analysis showed that most of Rosaceae-common transcripts supported the clustering of P. avium with M. domestica (2812 or 65%; Fig. ​Fig.4c4c upper panel and Fig. ​Fig.4d4d black dots) or R. chinensis with F. veseca (1436 or 33%; Fig. 4c lower panel and Fig. 4d blue dots) (Additional file 7: Figure S2).

Open in a separate window

Fig. 4

Identification and characterization of Rosaceae-common potential coding gene. a Venn diagram shows the Rosaceae-common and Rosa-specific transcripts. Note that, except Rosa, transcripts specific for other genera were not identified (marked with na). For that we are not interested in other share sets. b GO enrichment analysis of the 4447 Rosaceae-common transcripts (http://bioinfo.cau.edu.cn/agriGO). X-axis shows the enrichment fold of specific GO terms in comparison with the background. BP, CC and MF mean biological process, cellular component and molecular function separately. The area indicates gene counts. c Representive phylogentic topologies based on 4447 Rosaceae-common genes. Upper panel indicates about 65% topologies (2812) supporting the clustering of Prunus with Malus, while the topology in lower panel is supported by 33% genes (1436). Numbers on branches indicate distance. d Three-dimensional plots for the genetic distances of the 4447 transcripts between Rosa and Fragaria (X-axis) /Malus (Y-axis) /Prunus (Z-axis). Black and blue dots mark the genes supporting the topologies in C (Black for upper panel and blue for lower panel), while gray dots show genes supporting other topologies. e Distribution and GO enrichment analysis of the 409 selected Rosaceae-common transcripts. Y-axis shows the enrichment fold of specific GO terms in comparison with the background. Only four GO items are significantly enriched (marked in orange color). f Clustered heat map comparing scaled expression values for the 409 selected Rosaceae-common transcripts. Yellow indicates higher while purple marks lower expression. Blue and red bars indicate membership in the identified transcription clusters

Characteristics of the Rosa-specific transcripts

Among the 164 Rosa-specific molecules, 147 had protein-encoding potential with 136 of them having annotation (65 were uncharacterized previously; Additional file 12: Table S9). About 72 of these Rosa-specific genes were expressed at relatively lower level in both leaf and SAM materials (Fig. 5a, green and black bars), while around 92 of them (~56.1%) were expressed higher in SAMs than leaf materials (Fig. ​(Fig.5a,5a, blue and yellow bars). No GO item was significant enriched. However, transcripts encoding for putative F-box family members were strongly enriched (24x; Fig. ​Fig.5b;5b; Additional file 12: Table S9). F-box proteins represent one of the largest super-families in eukaryotic organisms, and contain F-box motif to recognize substrate proteins in ubiquitin-proteasome pathway [32, 33]. F-box proteins play a pivotal role in many physiological activities such as cell-cycle progression, transcriptional regulation, hormone signaling, programmed cell death and cell signaling transduction. F-box proteins can be distinguished to three classifications: F-box/kelch-repeat type, F-box/LRR-repeat proteins, and other F-box proteins. At least seven Rosa-specific transcripts belonged to F-box/LRR-repeat type, which might play essentials roles in adaptation of biotic stresses (Additional file 12: Table S9). Coincidently, four TMV resistant genes, which are essential for defenses against various plant viruses [34, 35], were also significantly enriched in these Rosa-specific genes (Fig. ​(Fig.5b).5b). These indicate that roses might have evolved new defense-related proteins to against its specific biotic stresses. Further experiments will be required to address the biological roles of these Rosa-specific transcripts.

Open in a separate window

Fig. 5

Identification and characterization Rosa-specific transcripts. a Heat map comparing scaled expression values for the 164 Rosa-specific transcripts. Yellow indicates higher while purple marks lower expression. Blue, green, yellow, and black bars indicate membership in the identified transcription clusters. b F-box and TMV resistance protein were significantly enriched in Rosa-specific transcripts. X² tests were performed online (http://www.quantpsy.org/chisq/chisq.htm) by comparing the Rosa-specific transcripts number with those from Rosaceae-common, coreset1 and coreset2 genes. P values were corrected with Bonferroni correction

Discussion

As one of the most important horticultural plants, rose has its special biology. Continuous flowering, fragrance, flower shape, thorn and many traits not presenting in Populus and other woody model plants could be found in roses, hence roses are now becoming a model woody species for understanding the molecular mechanisms regulating these traits. Interestingly, the breeding of modern roses often involves frequent hybridization and polyploidization among species, which often feature stronger diseases resistance and cold resistance, better fragrance and lack of prickles [36–38]. On the other hand, tracing the processes and history of modern roses domestication and breeding remains a challenge [2].

In this report, we produced high quality transcriptome assemblies for R. chinensis ‘Old Blush’ (OB) and R. wichuriana ‘Bayes’ Thornless’ (BT) with transcripts N50 above 2kb and mean length about 1.3kb. By incorporating published data for OB, we generated an improved assembly with mean transcript length longer than 1.3kb. We identified 10,073 highly conserved transcripts (coreset1) between OB and BT with BUSCO assay confirming the high level of completeness with these assemblies. As coreset1 transcripts were based on a relatively stringent level of sequence identity (90%), they could directly be used to evaluate the differential expression of orthologous genes between species.

These assemblies were then explored to identify about 5959 transcripts shared by Rosa sp. Finally we detected 4447 transcripts shared among Rosaceae, and about 164 transcripts only present in the Rosa transcriptomes. Since comparative selection pressure analysis provides an efficient way for understanding the molecular signatures underlying morphological trait variation and domestication [39–41], we next identified the genes under positive selection in roses among these Rosaceae-common transcripts. This analysis detected 409 rose genes significantly selected with around 40% of them being highly expressed in SAM tissues, suggesting that they might play important roles in the development of young organs/tissues of roses. It’s noteworthy that nine of these positively selected genes belong to GO items related to DNA repair and stress adaptation (Fig. ​(Fig.4)4) [29–31, 42]. This strong enrichment seems be associated with a strong requirement in the biological processes upon genome damage/confliction due to switch of environmental conditions and/or frequent intra- and inter-species hybridization. The Rosaceae-common dataset contains about 504 single-copy transcripts overlapped with the BUSCO set of genes. These single-copy transcripts could be directly used to clarify the Rosaceae phylogenetic relationship (Fig. ​(Fig.4),4), a challenge likely caused by frequent hybridization, rapid radiation, polyploidization and domestication [43–47].

In contrast to the Rosaceae-common transcripts, the identification of Rosa-specific transcripts seemed very interesting. Although the domestication processes had been documented [2], the evolutionary history and molecular mechanisms controlling traits special for roses are still not clear [1]. About half of the 164 transcripts are uncharacterized or without known GO annotation. These Rosa-specific transcripts might be related to the phenotypes that have not been characterized in other species. It’s possible that sequences of these transcripts have significantly diverged but remain similar functions like the homologs in Arabidopsis and roses. Interestingly, we observed a strong enrichment of F-box and TMV resistant proteins. These proteins might play a role in processes involved in biotic stress adaptation of roses.

Conclusions

In this study, we provided better quality transcriptome assemblies for roses, and pinned out the genes might make rose special. In addition, we identified transcripts common in Rosaceae, and these should help us to clarify the phylogenetic relationships of Rosaceae plants.

Methods

Additional files

Acknowledgements

Abbreviations

Authors’ contributions

JH, XD, and SL conceived the research together with DL, SD, SW, and KT. SL, MZ, XJ, YS, and FC performed experiments and helped in data analysis. MZ, XD, YX, and JH analyzed data. JH wrote the manuscript with contributions from all authors. All authors have read and approved the final manuscript.

Notes

Contributor Information

Shubin Li, Email: moc.361@asornc.

Micai Zhong, Email: nc.ca.bik.liam@iacimgnohz.

Xue Dong, Email: nc.ca.bik.liam@euxgnod.

Xiaodong Jiang, Email: nc.ca.bik.liam@gnodoaixgnaij.

Yuxing Xu, Email: nc.ca.bik.liam@gnixuyux.

Yibo Sun, Email: nc.ca.bik.liam@obiynus.

Fang Cheng, Email: nc.ca.bik.liam@gnafgnehc.

De-zhu Li, Email: nc.ca.bik.liam@lzd.

Kaixue Tang, Email: moc.liamtoh@gnatxk.

Siqing Wang, Email: moc.anis@gniqisgnaw.

Silan Dai, Email: moc.anis@iadnalis.

Jin-Yong Hu, Email: nc.ca.bik.liam@gnoynijuh.

References