基于高通量测序的极小种群野生植物长梗杜鹃转录组分析

李太强 刘雄芳 万友名 李正红 起国海 李钰莹 刘秀贤 和 锐 马 艳 马 宏

(中国林业科学研究院资源昆虫研究所，昆明 650233)

基于高通量测序的极小种群野生植物长梗杜鹃转录组分析

李太强 刘雄芳 万友名 李正红 起国海 李钰莹 刘秀贤 和 锐 马 艳 马 宏*

(中国林业科学研究院资源昆虫研究所，昆明 650233)

为加强特有濒危植物长梗杜鹃资源的评价、保护和鉴定工作，及为今后遗传育种和改善其农艺性状提供有益参考。本研究采用新一代高通量测序平台Illumina HiSeq 4000对其转录组测序，得到的数据过滤后进行de novo组装并聚类去冗余，获得74 092个Unigenes，平均长度、N50、Q20、Q30以及GC含量分别为938 nt、1 616 nt、98.22%、95.20%和43.24%，其中1 Kb以上的Unigenes有23 879条。通过与七大功能数据库比对，分别有39 876(NR:53.82%)、38 065(NT:51.38%)、27 384(Swissprot:36.96%)、16 099(COG:21.73%)、30 401(KEGG:41.03%)、17 518(GO:23.64%)以及29 676(Interpro:40.05%)条Unigenes获得功能注释。长梗杜鹃转录组中的Unigenes根据GO功能大致可分为生物学过程、细胞组分和分子功能3大类56亚类，其中执行生物学过程的基因最多，占41.53%；与COG数据库比对，根据其功能大致可分为25类；以KEGG数据库为参考，可归为6个代谢通路大类、32类代谢途径，并发掘出176条与人类疾病相关的Unigenes，包括内分泌及代谢疾病(167条)和抗药性(9条)。根据注释结果共检测出39 418个CDS，未注释上的Unigenes使用ESTScan预测后获得3 194个CDS。同时,预出到1 488个编码TF的Unigenes，以及检测到57 927个SNP多态位点。该转录组分析为今后长梗杜鹃乃至杜鹃属植物功能基因挖掘与利用、基因克隆、抗性机理分析、遗传资源分类和进化、分子标记开发以及分子辅助育种等研究提供了基础数据和重要参考。

长梗杜鹃；转录组；Unigene；功能注释；编码序列；转录因子；单核苷酸多态性

我国具有世界上最丰富的杜鹃花野生资源，目前记录有580余种(包括Flora of China出版后被描述的一些新种)。因此，欧美至今流传着一句名言“没有中国的杜鹃花,就没有西方的园林”。以长江为界，长江以南种类最多，长江以北种类很少；云南最多，四川次之，西藏第三，离分布中心愈远，种类愈少[1～3]。

杜鹃花是世界著名观赏花卉，年产量和产值巨大。然而，相较于种类众多的野生资源，我国在杜鹃花品种选育方面明显落后，市场上常见的品种均为比利时、美国、英国、日本、德国等国家选育。近年来，我国的杜鹃花育种工作取得一定进展，培育了一些具有自主知识产权的新品种[4～5]。另一方面，由于发展旅游、商业开发等人为破坏的影响，不少野生杜鹃花种类的生存受到不同程度的威胁，其中不乏极度濒危的种类。因此，对杜鹃花野生资源中濒危且具有较高观赏价值的种类优先开展相关研究，已迫在眉睫。长梗杜鹃(Rhododendronlongipedicellatum)就是其中的典型代表。

长梗杜鹃系杜鹃属、杜鹃亚属(subg.Rhododendron)、越桔杜鹃组(sect.Vireya)、类越桔杜鹃亚组(subsect.Pseudovireya)常绿植物。其花朵质地厚实、花色鲜黄，比同属的滇越杜鹃(R.rushforthii)、羊踯躅(R.molle)、鲜黄杜鹃(R.xanthostephanum)、纯黄杜鹃(R.chrysodoron)、硫磺杜鹃(R.sulfureum)等花色更为纯正。国际上，杜鹃花的花色育种趋势为纯色花，特别是黄色和蓝色更为珍贵[6]。更令人称奇的是，与野生杜鹃花期多集中于3～6月不同，长梗杜鹃的自然花期为11月下旬至翌年的2月上旬。此外，本研究组近几年通过多次专门的野外调查，仅发现位于云南省麻栗坡县的4个野生居群，各居群植株数为80～350株，海拔介于1 183～1 316 m。因其生境狭窄、分布区单一、居群数量少，受人为干扰严重，加之种群内实生苗较少，自然更新较弱，为典型的极小种群野生植物(Plant Species with Extremely Small Populations,简称PSESP)，野生资源状况令人担忧，亟待保护。

鉴于此，本文以长梗杜鹃为研究材料，采用Illumina HiSeq 4000高通量测序技术先对其进行转录组测序、数据过滤与组装，再用生物信息学的方法对得到的大量Unigenes进行功能注释、代谢途径、CDS预测及SNP检测等分析，以期了解长梗杜鹃在一定时期的基因表达情况和功能总体特征，为转录组水平上的研究以及分子层面揭示其濒危机制奠定基础，亦为今后合理开发利用这一珍稀杜鹃花种类提供参考。

1 材料与方法

1.1 试验材料

试验材料取自云南省麻栗坡县(104°93′E，23°16′N)，选择5株距离15 m以上的长梗杜鹃样株，各株采集质量大体相同(5～6片)的当年生幼嫩叶片，迅速放入液氮中，带回实验室后于-80℃冰箱中保存备用。

1.2 长梗杜鹃叶片总RNA提取

各单株分别称取0.1 g叶片等量混匀后，在液氮中迅速研磨至粉末状，用RNeasy Plant Mini Kit(Qiagen,Vallencia,CA)试剂盒按照说明书的操作指南提取，提取出的总RNA使用DNaseⅠ将其DNA消化，样品浓度、RNA完整性、28S/18S用Agilent 2100检测，OD260/280用NanoDrop检测，其RIN=8.3>8，28S/18S=1.84>1.8，OD260/280=1.97>1.8，所提取RNA质量较高，完整性较好。

1.3 转录组测序

进行mRNA的富集、cDNA合成、加碱基“A”、连接接头和PCR扩增等步骤，构建cDNA文库；构建好的文库用Agilent 2100 Bioanalyzer和ABI StepOnePlus Real-Time PCR System质检，合格后使用Illumina HiSeq 4000平台进行测序。

1.4 数据过滤和de novo组装

为保证结果的可靠性，先对测序的原始数据进行过滤，包括去除包含接头的reads、未知碱基N含量大于5%的reads以及质量值低于15的碱基占该reads总碱基数的比例大于20%的reads，过滤后的reads称为“Clean Reads”。然后使用Trinity对clean reads进行de novo组装，包含3个独立模块：Inchworm(构建k-mer库，K=25。以k-mer间overlap长度等于k-1对种子进行延伸，直到不能再延伸，形成线性Contig)；Chrysalis(把可能存在可变剪切及其它平行基因的Contigs聚类，用reads验证，看每个Component的reads支持情况)以及Butterfly(合并在de Bruijn图中有连续节点的线性路径，以形成更长的序列)，Trinity的组装结果称为转录本。最后使用Tgicl进行聚类去冗余得到Unigenes。

1.5 Unigene注释、CDS、TF编码能力预测以及SNP检测

使用Blast[7]对Unigenes进行NT、NR、COG、KEGG以及SwissProt注释，使用Blast2GO[8]以及NR注释结果进行GO注释，使用InterProScan5[9]进行InterPro注释。根据功能注释结果，按照NR，SwissProt，KEGG，COG的数据库优先顺序，挑选Unigenes的最佳比对片段作为该Unigenes的CDS。未能注释上的Unigenes使用上一步预测的CDS作为模型进行建模，然后通过ESTScan[10]进行预测。使用getorf[11]检测Unigenes的开放阅读框(ORF)，在通过hmmsearch[12]将ORF比对到转录因子蛋白结构域(数据来源于PlantfDB)，然后根据PlantfDB描述的转录因子家族特征对Unigenes进行能力鉴定。最后使用HISAT[13]把clean reads比对到Unigenes，借助GATK[14]检测SNP。

2 结果与分析

2.1 长梗杜鹃转录组测序数据过滤

测序共获得58.30 Mb的Raw Reads，过滤后得到44.85 Mb的Clean Reads，占原始读长的76.93%，包含6.73 Gb Clean Bases。其中，Q20为98.22%，即质量值大于20的碱基数目有6.61 Gb；Q30为95.20%。从Clean reads的碱基含量分布图可以看出(图1:A)，除RNA-Seq测序导致前6 bp出现正常波动外，其余reads每个位置的碱基含量分布稳定，无AT或GC分离现象。其次，在Clean reads的碱基质量分布图中(图1:B)，颜色深的点质量值集中于30～40，说明达到该质量值范围的碱基数目较多，而低质量(质量值小于20)的碱基比例较低，表明测序质量较好，测序所得reads准确性较高。

图1 Clean reads的碱基含量和质量分布图Fig.1 Distribution of base composition and quality on clean reads

2.2 长梗杜鹃转录组de novo组装

通过组装得到94 906个转录本，组装的连续性N50较大，为1 470 nt，说明组装效果较好；平均长度、GC含量分别为864 nt和43.31%。进一步聚类去冗余得到N50、平均长度以及GC含量分别为1 616 nt、938 nt和43.24%的Unigenes 74 092条，序列信息达69 505 225 nt，其中clusters(同一个clusters里面有若干条相似度大于70%的Unigenes)为51 505条，singletons(单独的Unigenes)为22 587条。从组装的转录本长度分布来看(图2:A)，随着转录本长度的增加，对应的数目呈递减的趋势；且大多数集中于200～2 000 nt，占转录本总数的89.85%，其中≥1 000 nt的占28.51%，仅3.36%的转录本≥3 000 nt。对聚类后Unigenes的长度分布进行统计(图2:B)，有32.23%的序列长度大于1 000 nt，3 000 nt及以上的序列仅占4.12%。表明聚类去冗余的效果较好，所得Unigenes平均长度、N50以及1 000 nt以上的序列数都有所提高，对Trinity组装结果起到了很好的优化作用。

图2 长梗杜鹃转录本和Unigene的长度分布图Fig.2 Transcripts and Unigenes length distribution for R.longipedicellatum

2.3 长梗杜鹃转录组Unigene的NR功能注释

将去冗余后的74 092条Unigenes比对到NR等七大功能数据库，最终被任一数据库注释上的序列共有45 538条，注释率为61.46%。其中，在NR库中(E≤1e-5)获得注释的Unigenes数量最多，有39 876条(占Unigenes总数的53.82%)，且与其它物种已知基因序列具有不同程度的匹配性。如注释物种分布图(图3)所示：注释为葡萄(Vitisvinifera)相关基因的序列最多，达10 261条(占NR库中被注释Unigenes的25.73%)，其次，序列同源性大于4%的近缘物种还有芝麻(Sesamumindicum)、罗布斯塔咖啡树(Coffeacanephora)、可可(Theobromacacao)以及莲(Nelumbonucifera)，其余近3/5的被注释Unigenes分布于其它511个物种中。由于缺乏长梗杜鹃基因组和转录组信息，尚存一定数量的Unigenes在NR库中仍未能获得注释。

图3 长梗杜鹃转录组Unigene的NR注释物种分布图Fig.3 NR annotated species distribution of Unigenes of transcriptome for R.longipedicellatum

2.4 长梗杜鹃转录组Unigene的COG注释及其分类

将获得的长梗杜鹃Unigenes比对到COG蛋白质直系同源数据库(E≤1e-5)，根据注释结果预测Unigenes可能的功能，并对其进行分类统计(图4)。

结果表明，有16 099条Unigenes(占总数的21.73%)被注释到25个COG的功能分类中，其各个分类的基因表达丰度各不相同；共包含26 605个功能注释信息，基本涵盖了长梗杜鹃大部分的生命活动。其中，一般功能基因(General function prediction only)为最大类别，占COG库总功能注释信息的17.08%；其次是转录(Transcription)功能，复制、重组和修复(Replication,recombination and repair)，翻译后修饰、蛋白质折叠和分子伴侣(Posttranslational modification,protein turnover and chaperones)以及信号传递机制(Signal transduction mechanisms)，分别占总功能注释信息的10.22%、8.75%、7.39%和6.76%。而核结构类基因(Nuclear structure)和胞外结构类基因(Extracellular structures)较少，仅有4(0.02%)和9(0.03%)个。另外有1 251个功能注释信息未明确其生物学功能，占4.7%。

图4 长梗杜鹃转录组Unigene的COG功能注释分布统计图Fig.4 COG functional annotation distribution of Unigenes of transcriptome for R.longipedicellatum

图5 长梗杜鹃转录组Unigene的GO功能分类统计图Fig.5 GO functional classification of Unigenes of transcriptome for R.longipedicellatum

2.5 长梗杜鹃转录组Unigene的GO注释及其分类

根据NR注释结果进行GO基因功能分类(参数为Blast2GO软件默认值)，综合描述长梗杜鹃叶片中相关基因的生物学特征，从宏观上了解其表达基因的功能分布情况，对获得相应GO条目的Unigenes进行统计分析(图5)。可知，有17 518(23.64%)条Unigenes得到92 861个GO注释，平均每条5.3个；近3/4 Unigenes在GO库中无同源基因与之匹配，进一步说明当前长梗杜鹃含有的基因信息较为缺乏。

获得的GO注释可划分为生物过程、细胞组分和分子功能三大类，其中执行生物学过程功能的注释最多(占41.53%)，涉及分子功能的最少(占20.73%)。3个大的类别又被划分为更加详细的56个小类别。依据其参与的生物学过程，38 567个GO条目被划分为24个亚类，其中，以细胞过程(cellular process)最多，占该亚类的24.22%；其次是代谢过程(metabolic process，23.90%)和单一有机体过程(single-organism process，16.65%)；而以细胞杀伤(cell killing)最低，仅占生物过程全部的0.002 6%。根据其在细胞中所处的位置，35 048个GO条目被划分为18个亚类，其中，以细胞(cell)、细胞组分(cell part)以及膜(membrane)较多，Unigenes分别占该类型的20.06%、19.90%和17.38%，而以细胞外基质组分(extracellular matrix component)和细胞外基质(extracellular matrix)较少，仅占0.003%和0.014%。分子功能包含的14个亚类中，以催化活性(catalytic activity)和结合(binding)居多，分别占该类别的45.71%和40.86%，而以蛋白标签(protein tag)和金属伴侣蛋白活性(metallochaperone activity)较少，分别占0.016%和0.021%。从这些结果可以看出，长梗杜鹃在细胞活动、代谢活动的基因表达丰度较高，表明长梗杜鹃自身具有较强的代谢能力。

2.6 长梗杜鹃转录组Unigene的KEGG代谢通路分析

将长梗杜鹃的Unigenes序列映射到KEGG代谢库(E≤1e-5)，根据注释信息对其可能参与或涉及的代谢途径进行统计分析。结果表明，30401条Unigenes获得注释(41.03%)；其参与的代谢通路可归为6大类别、21个亚类(图6)。

图6 长梗杜鹃转录组Unigene的KEGG功能注释分布统计图Fig.6 KEGG functional annotation distribution of Unigenes of transcriptome for R.longipedicellatum

由图6可知，6个代谢通路大类中，代谢(Metabolism)和遗传信息处理(Genetic Information Processing)相关的通路所占比例最高，分别为60.79%和23.25%，而以人类疾病(Human Diseases)相关的通路最少，仅占0.57%，环境信息处理(Environmental Information Processing)、生物系统(Organismal Systems)以及细胞过程(Cellular Processes)相关的通路分别占4.84%、5.04%和5.52%。进一步细分为21类代谢途径，其代谢包含的11个亚类中，全局整体映射(Global and overview maps)所占比例最高，为38.36%，其次是碳水化合物代谢(Carbohydrate metabolism，14.65%)和氨基酸代谢(Amino acid metabolism，8.12%)，而萜类化合物和聚酮化合物代谢(Metabolism of terpenoids and polyketides)所占比例最低，为4.10%。细胞过程和生物系统仅分别包括运输和代谢(Transport and catabolism)以及环境适应(Environmental adaptation)1个亚类，其余遗传信息处理、人类疾病和环境信息处理相关的通路分别包括4、2、2个亚类，这些亚类中，翻译，折叠、分类和降解(Folding,sorting and degradation)，所占比例较高，分别占总注释量的9.38%和6.78%；而内分泌及代谢疾病(Endocrine and metabolic diseases，0.537%)以及抗药性(Antimicrobial resistance，0.029%)较少。总的KEGG注释结果表明，该时期长梗杜鹃代谢活动和遗传信息处理能力较强。

2.7 长梗杜鹃转录组Unigene的CDS预测

对编码序列(Coding Sequence,CDS)进行预测分析能为长梗杜鹃的基因功能发掘利用和遗传连锁图谱构建提供重要参考。本研究按NR、Swissprot、KEGG和COG注释库的顺序共检测出39 418个Unigenes的最佳比对片段作为CDS，其总长度为35 487 000 nt，平均长度、N50以及GC含量分别为900 nt、1 350 nt和46.10%；未注释上的Unigenes使用ESTScan预测后获得3 194个CDS，其总长度为1 046 334 nt，平均长度、N50以及GC含量分别为327 nt、321 nt和49.31%，其长度分布如图7所示。从图7A中可看出，根据注释结果检测出的CDS序列较长，1 000 nt以上的占33.34%，其中1 000～2 000 nt的占总数24.61%，2 000～3 000 nt的占6.13%，3 000 nt以上的仅占2.6%。而通过ESTScan预测后获得的CDS序列较短，大部分集中在200～800 nt(占96.87%)，且大于1 600 nt的序列仅6条(图7：B)。

图7 长梗杜鹃转录组Unigene的CDS长度分布图Fig.7 CDS length distribution of Unigenes of transcriptome for R.longipedicellatum

2.8 长梗杜鹃转录组Unigene的TF编码能力及SNP检测

本研究共预测出60个基因家族共1 488个编码转录因子(Transcription Factor,TF)的Unigenes(图8)；其中MYB、MYB-related、AP2-EREBP、bHLH以及mTERF类转录因子为多基因家族，分别占总预测量的16%、12.03%、7.19%、6.72%和6.05%，而ULT、LFY、DBP所占比例较少，均为0.07%。同时，还预测出80条与非生物胁迫有关的TF，包括NAC(67条)和HSF(13条)；72条与生物胁迫的抗性密切相关的TF，包括WRKY转录因子超家族(58条)和bZIP类转录因子(14条)以及29条可能参与调节花发育的MADS框基因。这些TF的发现表明长梗杜鹃转录调控的复杂性，也为后续长梗杜鹃逆境相关基因的表达调控以及性状的利用和改良等研究提供了新的思路。

此外，还检测出57 927个单核苷酸多态性(Single Nucleotide Polymorphsims,SNP)位点(图9)，且在Unigenes序列上呈现不均匀分布，转换类型高于颠换。其中，转换类型(A-G和C-T)有36 171个，占总变异数的62.44%；颠换类型(A-C、A-T、C-G和G-T)有21 756个，占37.56%。在所有类型中，C-T变异的SNP数目最多(50.02%)，与杜玮南等[15]所得结论一致，SNP多发生在C和T之间。长梗杜鹃SNP分析，结合Unigenes的表达量及表型，可在mRNA层面的基因型差异区别于其它物种，对了解其亲属关系和进化关系具有重要的生物学意义，也为长梗杜鹃SNP标记的后续分子实验验证及其在自然种群中的应用奠定了基础。

图8 长梗杜鹃转录组Unigene的转录因子家族分类图Fig.8 Transcription Factor family classification of Unigenes of transcriptome for R.longipedicellatum

图9 长梗杜鹃转录组Unigene的SNP变异类型分布图Fig.9 SNP variant type distribution of Unigenes of transcriptome for R.longipedicellatum

3 讨论

随着高通量测序技术的快速发展以及测序平台的升级改进，通过转录组测序获取大量EST、发掘功能基因成为一种非常高效、可靠且高输出的重要手段，既降低了成本，又保证了测序的通量和所得序列的长度与质量，适合缺乏基因组信息的非模式物种开展生物信息学相关研究。虽然杜鹃属植物种类超过1 000种，但基因和基因组学研究相对滞后，遗传信息较为缺乏。与云锦杜鹃(R.fortunei)根系(4.66 Gb)转录组相比[16]，本研究获得的信息量较饱和，共产生6.73 Gb的数据量，包含74 092条Unigenes，为后续Unigenes的功能注释和分类、CDS预测、有利基因克隆及功能研究等方面奠定了坚实的基础。De novo组装后获得的Unigenes平均长度、N50值均较大，分别为938 nt和1 616 nt，远高于其它杜鹃属植物如酒红杜鹃(R.catawbiense)[17]、沼泽杜鹃(R.viscosum)[18～19]、云锦杜鹃[16]等，说明本研究所得序列中长片段较多，组装效果较好；碱基Q20和Q30分别为98.22%、95.20%。当Q30值在80%以上就认为测序质量非常可靠[20]，而碱基Q20与碱基识别的错误率呈对数相关，其表示每100个序列碱基中仅有一个出错的概率[21]。长梗杜鹃转录组Unigenes在1 kB以上的有23 879条，且clusters有51 505条，singletons有22 587条，所得clusters数目为单独Unigenes的2倍多，进一步说明所得序列的高效性。

将长梗杜鹃Unigenes序列比对到NR、GO、COG、KEGG等7大功能数据库，共有45538条获得注释，占总体的61.46%，且所得Unigenes的表达量FPKM值集中在0.2～28.98，表明本研究检测到长梗杜鹃低表达水平基因的比例很高。在NR数据库成功注释的Unigenes中，25.73%的Unigenes注释为与葡萄同源的基因，比例远远高于其它被注释的515个物种,其可能是由于葡萄具有参考基因组且与长梗杜鹃在生活史和进化关系上较为接近的缘故。而与其比对上的17个同属物种(包括栽培品种)仅120条Unigenes，其中映山红(R.simsii)、台北杜鹃(R.kanehirai)、锦绣杜鹃(R.pulchrum)、酒红杜鹃以及台湾杜鹃(R.formosanum)相对较多，分别有54(0.14%)、14(0.04%)、13(0.03%)、10(0.03%)和7(0.02%)条，其余12个种与之比对上的基因数所占比例均≤0.01%，这可能是现有杜鹃属植物基因组信息较为缺乏或是所得Unigenes多为长梗杜鹃特有的新功能基因而未能获得相应匹配。总的未被注释Unigenes有28 554条，占总体的38.54%，这些片段可能是因为较短或是本身为非编码序列而未能在7大功能数据库中获得同源性比对[22]。在28 554个无匹配序列中，长度在300～1 000 nt的Unigenes有27 398个，占96%，且随着序列长度增加，无同源序列的Unigenes比例明显降低，与贾新平[20]所得结论一致。转录组中Unigenes的序列越长，获得注释信息的可能性就越高。

利用COG数据库对所有Unigenes序列进行比对，获得了大量长梗杜鹃某一时期生长发育过程中表达的基因，16 099条序列共26 605个功能注释信息被划分为25类，一般功能、转录和翻译的基因表达丰度较高。通过GO数据库比对，根据其参与的生物学过程、构成的细胞组分和实现的分子功能分为56个小类别，以细胞活动、代谢活动和催化活性较为丰富，表明长梗杜鹃具有较强的代谢能力。根据KEGG数据库分析其可能参与或涉及的代谢通路，共30 401个Unigenes注释到6个代谢通路大类、21类代谢途径，其中定位到代谢相关通路的基因数最多，约占总注释量的3/5，进一步证实长梗杜鹃这一时期具有较强的代谢活动。更为值得一提的是，KEGG的代谢通路分析发掘出176条与人类疾病相关的Unigenes，包括内分泌及代谢疾病(167条)和抗药性(9条)。长梗杜鹃KEGG代谢途径的分析为进一步探索其抗性机理、挖掘花期调控相关功能基因，加快在引种驯化和杂交育种等方面的研究开发提供了分子依据。

转录因子在转录水平上与DNA发生相互作用的多样性决定了其生物学功能的多样性，本研究根据组装结果预测出60个基因家族共1 488个编码转录因子的Unigenes，这些家族在长梗杜鹃生长发育过程中发挥着不同的作用。同拟南芥、棉花、草莓和番茄等大多数植物一样[23]，MYB基因家族在长梗杜鹃中普遍存在，是其中数量最大、功能最多样的转录家族之一，可能与其代谢调控、激素的信号传导、细胞周期调控有关。同时，107个植物特有转录因子AP2-EREBP和100个植物第二大类转录因子bHLH的发现，也为长梗杜鹃生长发育和应答外界环境信号过程、进化和基因间的相互作用以及抗逆调控网络的研究奠定了基础。EST文库是筛选SNP的重要来源，本研究通过转录组测序分析，共检测出57 927个SNP多态位点，包括36 171个转换类型，21 756个颠换类型，将为今后开展长梗杜鹃遗传图谱构建、重要性状基因定位、遗传多样性分析、群体选择分析、品种鉴定以及分子标记辅助选择育种等研究提供了重要参考。

4 结论

本研究采用Illumina HiSeq 4000最新高通量测序平台建立了新物种长梗杜鹃转录组数据库，获得了74 092条与该物种生长发育相关的Unigenes，通过生物信息学方法对其进行功能注释、代谢相关通路、CDS预测、TF编码能力预测、SNP检测以及表达量计算等分析。综合COG、GO和KEGG三大功能数据库的注释结果，表明长梗杜鹃在该时期转录、复制和翻译等基因表达丰度较高，具有较强的代谢活动和遗传信息处理能力。所得转录组信息作为杜鹃属植物基因组序列的重要补充，进一步丰富了该属植物的基因表达信息平台，为长梗杜鹃乃至同属植物基因克隆及功能分析、SSR、SNP分子标记的规模化开发奠定了数据基础。同时，随着新一代高通量测序技术的广泛应用，将长梗杜鹃Unigenes与其它同属植物或近缘物种进行比较分析，可为进一步研究其进化和有利基因的筛查提供更为重要的信息。

1.方瑞征,闵天禄.杜鹃属植物区系的研究[J].云南植物研究,1995,17(4):359-379.

Fang R Z,Min T L.The floristic study on the genusRhododendron[J].Acta Botanica Yunnanica,1995,17(4):359-379.

2.中国科学院中国植物志编辑委员会.中国植物志[M].北京:科学出版社,1999:1-3.

Delectis Florae Reipublicae Popularis Sinicae,Agendae Academiae Sinicae Edita.Flora reipublicae popularis sinicae[M].Beijing:Science Press,1999:1-3.

3.张长芹,高连明,薛润光,等.中国杜鹃花的保育现状和展望[J].广西科学,2004,11(4):354-359,362.

Zhang C Q,Gao L M,Xue R G,et al.A general review of the research and conservation statue of ChineseRhododendron[J].Guangxi Sciences,2004,11(4):354-359,362.

4.张长芹,罗吉凤.杜鹃花新品种‘金踯躅’和‘紫艳’[J].园艺学报,2002,29(5):502.

Zhang C Q,Luo J F.NewRhododendronvarieties-‘Jinzhizhu’ and ‘Ziyan’[J].Acta Horticulturae Sinica,2002,29(5):502.

5.张长芹,罗吉凤,冯宝均.杜鹃花新品种‘朝晖’和‘红晕’[J].园艺学报,2002,29(3):296.

Zhang C Q,Luo J F,Feng B J.NewRhododendronhybrid-‘Zhaohui’ and ‘Hongyun’[J].Acta Horticulturae Sinica,2002,29(3):296.

6.兰熙,张乐华,张金政,等.杜鹃花属植物育种研究进展[J].园艺学报,2012,39(9):1829-1838.

Lan X,Zhang L H,Zhang J Z,et al.Research progress ofRhododendronbreeding[J].Acta Horticulturae Sinica,2012,39(9):1829-1838.

7.Altschul S F,Gish W,Miller W,et al.Basic local alignment search tool[J].Journal of Molecular Biology,1990,215(3):403-410.

8.Conesa A,Götz S,García-Gómez J M,et al.Blast2GO:a universal tool for annotation,visualization and analysis in functional genomics research[J].Bioinformatics,2005,21(18):3674-3676.

9.Quevillon E,Silventoinen V,Pillai S,et al.InterProScan:protein domains identifier[J].Nucleic Acids Research,2005,33(S2):W116-W120.

10.Iseli C,Jongeneel C V,Bucher P.ESTScan:a program for detecting,evaluating,and reconstructing potential coding regions in EST sequences[C].//Proceedings of the Seventh International Conference on Intelligent Systems for Molecular Biology.Menlo Park:AAAI Press,1999:138-148.

11.Rice P,Longden I,Bleasby A.EMBOSS:the European molecular biology open software suite[J].Trends in Genetics,2000,16(6):276-277.

12.Mistry J,Finn R D,Eddy S R,et al.Challenges in homology search:HMMER3 and convergent evolution of coiled-coil regions[J].Nucleic Acids Research,2013,41(12):e121.

13.Kim D,Langmead B,Salzberg S L.HISAT:a fast spliced aligner with low memory requirements[J].Nature Methods,2015,12(4):357-360.

14.Mckenna A,Hanna M,Banks E,et al.The genome analysis toolkit:a MapReduce framework for analyzing next-generation DNA sequencing data[J].Genome Research,2010,20(9):1297-1303.

15.杜玮南,孙红霞,方福德.单核苷酸多态性的研究进展[J].中国医学科学院学报,2000,22(4):392-394.

Du W N,Sun H X,Fang F D.The research development of single nucleotide polymorphism[J].Acta Academiae Medicinae Sinicae,2000,22(4):392-394.

16.Wei X Y,Chen J J,Zhang C Y,et al.Differential gene expression inRhododendronfortuneiroots colonized by an ericoid mycorrhizal fungus and increased nitrogen absorption and plant growth[J].Frontiers in Plant Science,2016,7:1594.

17.Wei H,Dhanaraj A L,Rowland L J,et al.Comparative analysis of expressed sequence tags from cold-acclimated and non-acclimated leaves ofRhododendroncatawbienseMichx[J].Planta,2005,221(3):406-416.

18.Susko A Q,Rinehart T,Bradeen J,et al.EST-SSR development and validation in deciduous azalea breeding populations from transcriptome sequencing[C].//Plant and Animal Genomic XXIII Conference.San Diego,CA:University of Minnesota,2015.

19.Susko A Q.Phenotypic and genetic variation for rhizosphere acidification,a candidate trait for pH adaptability,in deciduous azalea(Rhododendronsect.Pentanthera)[D].Minnesota:University of Minnesota,2016.

20.贾新平,孙晓波,邓衍明,等.鸟巢蕨转录组高通量测序及分析[J].园艺学报,2014,41(11):2329-2341.

Jia X P,Sun X B,Deng Y M,et al.Sequencing and analysis of the transcriptome ofAspleniumnidus[J].Acta Horticulturae Sinica,2014,41(11):2329-2341.

21.Ewing B,Hillier L,Wendl M C,et al.Base-calling of automated sequencer traces usingPhred.Ⅰ.accuracy assessment[J].Genome Research,1998,8(3):175-185.

22.蔡年辉,邓丽丽,许玉兰,等.基于高通量测序的云南松转录组分析[J].植物研究,2016,36(1):75-83.

Cai N H,Deng L L,Xu Y L,et al.Transcriptome analysis forPinusyunnanensisbased on high throughput sequencing[J].Bulletin of Botanical Research,2016,36(1):75-83.

23.牛义岭,姜秀明,许向阳.植物转录因子MYB基因家族的研究进展[J].分子植物育种,2016,14(8):2050-2059.

Niu Y L,Jiang X M,Xu X Y.Research advances on transcription factor MYB gene family in plant[J].Molecular Plant Breeding,2016,14(8):2050-2059.

The Cultivation Project of Technology Innovation Talent of Yunnan Province

introduction：LI Tai-Qiang(1993—),male,graduate student,mainly engaged in the study of conservation biology inRhododendron.

date:2017-03-29

TranscriptomeAnalysisforRhododendronlongipedicellatum(PlantSpecieswithExtremelySmallPopulations)BasedonHighThroughputSequencing

LI Tai-Qiang LIU Xiong-Fang WAN You-Ming LI Zheng-Hong QI Guo-HaiLI Yu-Ying LIU Xiu-Xian HE Rui MA Yan MA Hong*

(The Research Institute of Resource Insects,Chinese Academy of Forestry,Kunming 650233)

To strengthen the research of resources evaluation, protection and identification of endemic and endangered species ofRhododendronlongipedicellatum, and to provide a helpful reference for genetic breeding and improvement of its agronomic traits, the transcriptome was sequenced by using Illumina Hiseq 4 000, in total, 74 092 Unigenes with an average length of 938 nt, N50of 1616 nt, Q20of 98.22%, Q30of 95.20% and GC content of 43.24% were obtained by de novo assembly and cluster with filtered data, and there were 23 879 Unigenes with more than 1 kB. Then, the Unigenes were annotated by 7 functional databases, and finally, 39 876(NR: 53.82%)，38 065(NT: 51.38%), 27 384(Swissprot: 36.96%), 16 099(COG: 21.73%), 30 401(KEGG: 41.03%), 17 518(GO: 23.64%), and 29 676(Interpro: 40.05%) Unigenes were annotated. The Unigenes were roughly divided into three functional categories(i.e. biological processes, cellular components and molecular function) and 56 sub-categories according to GO function. Most of the genes performed biological processes. KEGG functional annotation analysis showed that Unigenes could be grouped into 6 categories, 32 metabolic pathways. 176 Unigenes relating to human diseases were detected, including endocrine and metabolic diseases(167) and antimicrobial resistance(9). The 39 418 CDS were detected with functional annotation results, and after 3 194 CDS also were predicted by ESTScan with the remaining Unigenes. The 1 488 Transcription Factor(TF) coding Unigenes were predicted and 57 927 SNP polymorphic loci were detected. The analysis of the transcriptome could lay a foundation for further study of functional gene discovery and utilization, resistance mechanism analysis, classification and evolution of genetic resources, molecular marker development and molecular assisted breeding ofR.longipedicellatumand other congeneric species.

Rhododendronlongipedicellatum;transcriptome;Unigene;gene annotation;CDS;Transcription Factor;SNP

“云南省技术创新人才”培养对象项目(2016HB007)

李太强(1993—)，男，硕士研究生，主要从事杜鹃属植物保护生物学研究。

* 通信作者：E-mail:hortscience@163.com

2017-03-29

* Corresponding author：E-mail:hortscience@163.com

S685.21

A

10.7525/j.issn.1673-5102.2017.06.004