Transcriptome Profiling of Lisianthus (Eustoma Grandiflorum) under Drought Stress

An Xia ^1*　Zhu Qiang ¹　Lou Xuping ²　Li Lufeng ²　Chen Jie ³　Liu Tingting ¹　Li Wenlue¹　Luo Xiahong ¹　Zhu Guanlin ¹　Yu Lijun ¹

1 Zhejiang Xiaoshan Institute of Cotton and Bast Fiber Crops Research, Zhejiang Institute of Landscape Plants and Flowers, Hangzhou, 311251; 2 Hangzhou Xiaoshan District Agricultural Science and Technology Research Institute, Hangzhou, 311202; 3 Huazhong Agricultural University, Wuhan, 430070

* Corresponding author, anxia@zaas.ac.cn

Abstract　Lisianthus (Eustoma grandiflorum) is an important ornamental plant native to North America. So far, the number of basic molecular research reports on Lisianthus is relatively small, especially the research related to the molecular mechanism of response to drought stress. In this experiment, transcriptome second-generation sequencing technology was used to study Lisianthus seedlings under drought stress. The results showed that the transcriptome sequencing data had good quality control results. In order to obtain genetic information on this basis, the sequences were spliced using software such as Trinity and Corset, a total of 102 014 non-redundant genes were obtained, including 2 929 coding genes which were found to belong to 79 transcription factor families after analyzing. Annotating all genes, Zhongguo coffee was found to be the species with the highest comparison similarity. Finally, simple sequence repeat (SSR) analysis was performed on the gene sequences, and all SSR information contained in non-redundant genes and transcription factor coding genes had been obtained. The results of this study can provide candidate molecular resources for subsequent studies on Lisianthus 's response to drought stress in the future.

Keywords　Lisianthus; Drought stress, Transcriptome

洋桔梗(Eustoma grandiflorum)别名草原龙胆，原产北美地区，是龙胆科(Gentianaceae)多年生草本植物。洋桔梗是重要的观赏植物，其花器官形态漂亮并且瓶插寿命长，已经成为国际越来越重要的鲜切花之一。洋桔梗对栽培环境十分敏感，国内栽培设施及技术相对不足。对洋桔梗生长发育过程中的相关分子机理进行研究，有助于为洋桔梗栽培新技术开发提供理论支持。

植物在响应干旱胁迫的过程中涉及一些功能基因和调节基因的差异表达，形成了一套复杂的信号调控网络，从而影响植物体内一系列的生理生化反应。干旱能够影响洋桔梗花茎伸长，但是目前为止关于洋桔梗如何响应干旱胁迫的相关研究鲜有报导。本研究通过对干旱胁迫处理下的洋桔梗植株进行转录组测序，共获得6.43 Gb高质量测序数据。对该测序结果进行拼接，获得了102,014条基因。对所获得的基因进行注释，发现和现有其它研究基础较好的物种之间亲缘关系都较远。因此，该项研究不仅能为后续针对洋桔梗的相关研究提供大量分子资源，还能为龙胆科其它物种高通量测序研究提供借鉴。

1结果与分析

1.1转录组数据获得、处理与拼接

首先对本研究中的转录组测序数据进行严格的质量控制，主要去除接头以及低质量序列信息，总共获得高质量测序数据(Clean reads) 6.43 Gb，然后采用Trinity软件(Grabherr et al., 2011)对本研究中获得的这些高质量测序数据进行拼接后，共得到132 929条转录本(Transcripts)。随后，使用软件Corset对获得的高质量测序数据进行分析，将所有reads与转录本进行比对并进行层次聚类，然后得到102 014个非冗余的基因(Unigenes)。这些转录本和基因在序列长度方面具有相似分布规律(图1A; 图1B)，且转录本与非冗余基因之间的数量差异主要体现在较短序列(<500 bp)上(图1C)。该结果表明本次测序结果质量较高，主要表现为在较长序列(>1 000 bp)范围内，转录本与非冗余基因数目没有显著差别。

图 1 转录组序列长度分布 注: A: 转录本序列长度分布; B: 基因长度分布; C: 转录本和基因长度分布统计 Figure 1 Lengths distribution from transcriptomic data Note: A: Lengths distribution of transcripts; B: Lengths distribution of unigenes; C: Statistical results of transcripts and unigenes

1.2基因注释和功能分类

为了更科学系统地分析转录组测序获得的序列所涉及的基因功能信息，本研究对拼接得到的102 014条基因，分别采用不同的公共数据库对获得的基因进行注释。在不同数据库中获得注释基因的统计信息见表1，其中有79.94%的基因在至少一个数据库中有注释结果，而来源于Nr数据库的注释结果最多，占所有基因的76.52%。选取基因在Nr、Nt、pfam、GO和KOG五个数据库的注释结果进行分析表明，特异在Nr数据库中得到注释信息的基因数有13,035条，远远多于Nt数据库中特异注释的741条基因和KOG数据库中特异注释的8条基因，在pfam和GO数据库中没有特异注释的基因(图2A)。因此，转录组数据在Nr数据库中的注释信息更全面。进一步分析Nr数据库注释结果，有接近一半(47.3%)的序列与目标序列具有较高相似度(超过80%, 图2B)，且大量序列(占比60.3%)的比对结果e值小于1e-60 (图2C)。这些序列物种注释结果中，比对到最多的物种为中果咖啡(Coffea canephora, 占比39.7%)，并且有更大比例的序列(44.1%)比对到其它(Other)物种(图2D)。在KOG数据库分类中，有25个KOG类别被不同数量基因所注释，共包含28 696条基因。其中被注释基因最多的两个类别是O: 翻译后修饰：蛋白开关和分子伴侣(3 519条基因)和R: 总体功能预测(3 421条基因) (图3)，该结果与黄麻干旱胁迫转录组结果类似。与此同时，从代谢角度来看，这些基因被更多富集在“遗传信息处理”大类中的“翻译”、“代谢”大类中的“碳水化合物代谢”和“遗传信息处理”中的“折叠、排列和降解”等三个代谢条目中(图4)。最后，由于转录因子往往处于基因表达通路的上游，能够调控下游一系列基因表达从而在更大程度上影响洋桔梗响应干旱胁迫的程度。与之对应，本次转录组测序结果中共有2 929条基因可能编码转录因子，这些转录因子属于79个不同的转录因子家族(图5)。在这些转录因子家族中，包含预测基因数目最多的三个家族分别是bHLH家族(256条基因)MYB_related家族(216条基因)和bZIP家族(191条基因)。

表 1 基因注释成功率统计 Table 1 Statistical numbers on successfully annotated genes against different databases

图 2 转录组注释信息 注: A: 转录组结果比对到不同数据库的基因数目; B: Nr数据库比对结果序列相似度分布; C: Nr数据库比对结果e值分布; D: Nr数据库注释比对到最多的物种 Figure 2 Annotation of transcriptome Note: A: Numbers of unigenes annotated by different databases; B: Distribution of sequence similarities against the Nr database; C: Distribution of e-values against the Nr database; D: The most annotated species from our transcriptome data

图 3 转录组结果的KOG分类 Figure 3 The KOG classification of our transcriptome data

图 4 转录组结果的KEGG分类 注: 共分为五大类; A: 细胞进程; B: 环境信息处理; C: 遗传信息处理; D: 代谢; E: 有机系统 Figure 4 The KEGG classification of our transcriptome data Note: A total of five categories were included: A, Cellular processes; B, Environmental information processing; C, Genetic information processing; D: Metabolism; E: Organismal systems

图 5 转录因子数目统计 Figure 5 Numbers of transcription factors from different families

1.3分子标记开发

通过对转录组测序所获得的序列进行简单重复序列(simple sequence repeat, SSR)分析，共在21 329条基因序列(占所有基因序列的20.91%)中发现了25 468个SSR。这些SSR序列主要包括单碱基至六碱基不同程度的重复，以及复杂重复序列。如图6A所示，除了二碱基重复以外，SSR重复序列平均总长度随着重复单位的复杂性增加而递增，其中复杂重复单元的重复序列长度最长。在编码转录因子的2,929条基因中，有833条序列(占比28.44%)具有不同的SSR，这些SSR序列总长度也具有前述类似规律(图6B)。最后分析SSR可能位于基因的不同位置，表明对于所有基因来说，超过一半(51.46%)的重复序列可能横跨相邻的两个基因结构(5’非翻译区: utr5; 编码区: cds; 3’ 非翻译区: utr3)，而位于编码区的SSR所占比例最少(图6C)。对于转录因子编码基因来说，横跨两个基因结构的SSR大幅度减少(27.25%)，而位于编码区的重复序列所占比例依然最少(图6C)。后续可以通过开发这些SSR的特异引物，对特定基因，或者转录因子进行针对性更强的检测和研究。

图 6 简单序列重复(SSR)信息统计 注: A: 转录组中SSR序列长度统计; B: 转录因子unigene中SSR序列长度统计; C: SSR位于不同基因功能区域统计; 重复序列类型包括单碱基重复至六碱基重复(p1-p6)以及复杂重复单元(c); 这些SSR可能位于基因的5’非翻译区(utr5), 编码区(cds), 3’非翻译区(utr3)或者未知位置(undetermined)。 Figure 6 Statistical on simple sequence repeats (SSR) Note: A: Sequence lengths distribution of SSRs amongst the whole transcriptome data; B: Sequence lengths distribution of SSRs from transcription factor coding genes; C: Location of SSRs on varied districts of unigenes; The SSR units included mononucleotides to hexanucleotides (p1 to p6), and complex units (c); These SSRs may locate on the 5’ untranslated regions (utr5), 3’ untranslated regions (utr3), coding sequences (cds), or currently unknown positions (undetermined)

2讨论

洋桔梗是重要的观赏植物。然而，针对该物种的分子生物学研究基础较欠缺。目前为止，仅在早期构建过一个针对盐胁迫的差减文库(王继刚等, 2008)，鉴定了可能的差异性表达基因。然而，差减文库的通量一般较低，和现行高通量测序相比，远远不能满足研究的需求。在分子资源开发方面，也仅有早期针对花期进行的转录组测序(Kawabata et al., 2012)。在该项转录组测序中，所得到的63 401条contig仅有65%在NCBI数据库中得到了比对结果，小于本研究中的76.52% (表1)，说明随着测序技术的发展和拼接方法的成熟，本次测序结果的序列注释情况得到了提高。然而，可能是由于与洋桔梗的近源物种分子研究基础均较薄弱，本次转录组注释结果，比对到最高比例的物种(图2D)是茜草科(Rubiaceae)的中果咖啡(Coffea canephora)。茜草科和龙胆科同属龙胆目，因此该中果咖啡可能是和洋桔梗亲缘关系最近的有一定分子研究资源的物种。除此之外，洋桔梗转录组序列的比对结果零散地分布在其他物种上(图2D)。

曾经有研究对洋桔梗MADS家族基因进行鉴定(Ishimori and Kawabata, 2014)和功能研究(Li et al., 2015)。然而，在没有转录组或者基因组等高通量测序结果的支持下，进行基因功能研究，或者基因家族鉴定往往显得更加困难(Nakano et al., 2011)。本研究对洋桔梗在干旱胁迫下进行了转录组测序和并对测序数据进行初步分析，相关研究结果为后期研究洋桔梗响应干旱胁迫处理的分子机理提供相应数据支撑。与此同时，也有研究者完成了洋桔梗质体组测序(Yan et al., 2019)，相关分子资源为通过质体遗传信息调节花形态等农艺性状(Jin and Daniell, 2015)提供了相应基础。该质体组测序结果(Yan et al., 2019) 与龙胆科其它物种对应测序结果序列比对显示，洋桔梗与这些物种亲缘关系都更远。该结果从侧面印证了本次转录组测序注释结果中，比对到最多的物种是龙胆目茜草科下的中果咖啡(占39.7%)，而其余大部分序列信息均零散地比对到其它物种(图2D)。最后，本次转录组测序结果都将为后续研究，如鉴定挖掘干旱响应相关基因或者干旱胁迫相关的转录因子提供了分子数据。

3材料与方法

3.1植物材料

洋桔梗(Eustoma grandiflorum)品种“雪莱”在市场上购买。植株长至8 cm左右，对其进行干旱胁迫处理。处理36 h后，使用液氮取全株植物样品，用于总RNA提取。

3.2总RNA提取及文库构建

样品在用于RNA提取之前一直保存于-80℃超低温冰箱内。将样品取出并在液氮环境下充分研磨成粉末，使用天根公司的RNA提取试剂盒完成总RNA 提取并用于构建转录组文库。

3.3测序数据处理及转录本拼接

样品上机测序得到的直接数据为原始读数(Raw reads)，需要进行质量控制，去除不确定碱基测序结果和测序信息由于包含大于10%的接头从而导致质量不佳等，余下读数即为高质量测序信息。对高质量测序信息的处理方式主要为序列拼接(Grabherr et al., 2011)和层次聚类，所得到的基因信息即为非冗余基因。

3.4基因功能注释和分类

对经由上述步骤得到的非冗余基因，与四个序列数据库进行比对以进行注释。

3.5转录因子预测

在对非冗余基因序列进行注释的同时，使用在线工具(http://planttfdb.gao-lab.org/prediction.php)可以对这些基因信息可能编码产物进行预测，若编码转录因子则对其进行分类。

3.6 SSR分析

在基因的碱基序列层面，存在一些规律明确的序列特征。这些序列往往以简单的序列单元(单个到多个碱基, 甚至复杂碱基单元)为基础，重复出现多次，成为简单序列重复(SSR)。对于这些SSR一般使用在线工具(http://pgrc.ipk-gatersleben.de/misa/misa.html)进行预测。对于单碱基单元重复十次及以上，以二碱基为单元重复六次及以上，以及三碱基至六碱基为重复单位重复五次及以上的SSR均包含在统计范围内，而复杂重复序列中如果包含以上所列不同的重复单元，则每个重复单元分别满足上述要求。

作者贡献

安霞是本研究的实验设计和实验研究的执行人，完成数据分析，论文初稿的写作；朱强、楼旭平、李鲁峰、陈杰、柳婷婷、李文略、骆霞虹、朱关林、余利隽是实验设计参与者；安霞是项目的构思者及负责人，指导实验设计，数据分析，论文写作与修改。全体作者都阅读并同意最终的文本。

致谢

本研究由省科技特派员项目(梯田景区农家乐景观提升示范与休闲产品创意)资助。

参考文献

Grabherr M.G., Haas B.J., Yassour M., Levin J.Z., Thompson D.A., Amit I., Adiconis X., Fan L.,  Raychowdhury R., Zeng Q., Chen Z., Mauceli E., Hacohen Nir., Gnirke A., Rhind N., Palma F.D., Birren B.W., Chad Nusbaum., Lindblad-Toh K., Friedman N., and Regev A., 2013, Trinity: reconstructing a full-length transcriptome without a genome from rna-seq data, Nature Biotechnology, 29(7): 644-652

Ishimori M., and Kawabata S., 2014, Conservation and Diversification of Floral Homeotic MADS-box Genes in Eustoma grandiflorum, J. Japan, Soc. Hort. Sci., 83(2): 172-180

Jin S., and Daniell H., 2015, The engineered chloroplast genome just got smarter, Trends in Plant Science, 20(10): 622-640

Kawabata S., Li Y., and Miyamoto K., 2012, EST sequencing and microarray analysis of the floral transcriptome of Eustoma grandiflorum, Scientia Horticulturae, 144: 230-235

Li K.H., Chuang T.H., Hou C.J., and Yang C.H., 2015, Functional analysis of the FT Homolog from Eustoma grandiflorum reveals its role in Regulating A and C Functional MADS Box genes to control floral Transition and flower formation, Plant Mol. Biol. Rep., 33(4): 770-782

Nakano Y., Kawashima H., Kinoshit T., Yoshikawa H., and Hisamatsu T., 2011, Characterization of FLC, SOC1 and FT homologs in Eustoma grandiflorum: effects of vernalization and post-vernalization conditions on flowering and gene expression, Physiologia Plantarum, 141(4): 383-393

Wang J.G., Zhang K., Xu Q.J., and Li Y.H., 2008, Construction and Analysis of Eustoma grandiflorum Subtracted cDNA Library, Yuanyi Xuebao (Acta Horticulturae Sinica), 35(7): 1075 -1080. (王继刚, 张坤, 徐启江, 李玉花, 2008, 草原龙胆盐胁迫差减文库的构建及分析, 园艺学报, 35(7): 1075-1080)

Yan J., Cao Q., Wu Z., Chen S., Wang J., Zhou D., and Xie J., 2019, Complete plastome sequence of Eustoma grandiflorum (Gentianaceae), a popular cut flower, Mitochondrial DNA Part B, 4(2): 3163-3164