基于高通量测序的羊踯躅叶绿体基因组及SSR序列分析

吴林世，廖菊阳，刘 艳，李巧云，王 玲，张 娟，宋 胤，胡俊东

（1.湖南省植物园，湖南 长沙 410116；2.国家林业草原局 杜鹃工程技术研究中心，湖南 长沙 410116）

羊踯躅Rhododendron molle为杜鹃花科Ericaceae杜鹃花属Rhododendron的落叶灌木，是湖南省新优杜鹃属乡土植物。因其花色金黄，又名黄杜鹃；因其花有毒，亦名闹羊花[1]。羊踯躅树形优美，花繁色艳，是园林露地栽培的优良杜鹃花品种[2]。其不仅有良好的观赏价值，还具有一定的药用价值，可以用于祛风除湿，在治疗跌打损伤、偏正头痛、顽癣等方面也有一定的功效[3]。羊踯躅是一种微毒植物，可用于散瘀止痛，医学上也用作镇静剂[4-5]。目前，有关羊踯躅的研究报道主要集中于化学成分鉴定、药理作用、亲缘关系、繁育等方面[6-8]，有关羊踯躅叶绿体基因组的研究报道较为鲜见，仅见其全基因测序结果的相关报道[9]。叶绿体、线粒体、细胞核并称为三大遗传系统，叶绿体含DNA 聚合酶和RNA 聚合酶，还具备独立、完整的蛋白质合成系统，能自我完成RNA 的复制和转录。为了进一步补充完善羊踯躅基因组的序列特征，为羊踯躅新优园艺品种的鉴定、不同地理分布的羊踯躅亲缘关系分析提供参考，并为杜鹃属植物的系统进化、遗传多样性分析、基因工程等提供参考，本研究中对羊踯躅叶绿体基因组进行测序组装，根据已发表的杜鹃属植物基因组学的研究结果，确定羊踯躅在系统进化过程中的地位。

1 材料与方法

1.1 试验材料

羊踯躅叶片采自湖南省植物园（113°01′30″E，28°06′40″N[10]）。在8月中旬羊踯躅叶片生长旺盛时期，采集目标植株上健康无病虫害的叶片作为供试材料，送往南京集思慧远生物科技有限公司检测。

1.2 试验方法

1.2.1 羊踯躅叶绿体基因组测序

样本经硅胶干燥处理后送至实验室，检测前保存于-20 ℃低温冰箱中。严格执行Illumina 公司提供的叶绿体基因测序组装操作标准，得到测序文库，经过质检合格后，将文库采用Illumina Novaseq 平台进行测序。将经测序得到的原始数据（raw data）进行过滤，筛除接头序列和低质量读长（reads），获得高质量的数据（clean data）。按照羊踯躅叶绿体基因组序列将所得的高质量数据进行序列组装，获得叶绿体序列组装结果。

1.2.2 羊踯躅叶绿体基因组装

1）叶绿体DNA 测序序列提取。使用bowtie2 v2.2.4（http://bowtiebio.sourceforge.net/bowtie2/index.shtml）的very-sensitive-local 模式，将所得序列与南京集思慧远生物科技有限公司自建的数据库进行比对，比对一致的测序序列即目标叶绿体DNA 序列。

2）叶绿体基因序列组装。采用SPAdes v3.10.1（http://cab.spbu.ru/software/spades/）软件组装叶绿体基因组，k-mer 分别使用55、87、121，组装不依赖参考基因组[11]。组装过程中使用SPAdes、SSPACE v2.0、Gapfiller v2.1.1 等软件得到完整的pseudo genome 序列，再次进行基因组校正，参考叶绿体结构将校正后的pseudo genome 进行坐标重排，得到目标基因组序列（环状基因组序列）[12]。

1.2.3 羊踯躅叶绿体基因进化树构建

为分析羊踯躅在物种系统进化中的位置，从NCBI 中下载已知的19个杜鹃花科植物叶绿体基因组和2个外群基因组。19个杜鹃花科植物分别是白珠树属Gaultheria的4个物种，越橘属Vaccinium的5个物种，杜鹃属的10个物种，其中羊踯躅和溪畔杜鹃R.rivulare（该物种叶绿体基因组序列未上传NCBI）是自测物种；外群植物选择的是桤叶树科Clethraceae 桤叶树属Clethra的2个物种。将下载的叶绿体基因序列与羊踯躅叶绿体基因序列进行比对分析，构建杜鹃属植物的进化树。

2 结果与分析

2.1 羊踯躅叶绿体基因组的基本特征

将数据进行初步整理，得到羊踯躅叶绿体基因组图谱（图1）。由图1 可见，经测序组装得到的基因组序列属于典型的四段式结构，全长为197 877 bp。其中，大单拷贝区（LSC）长度为110 189 bp，小单拷贝区（SSC）长度为26 bp，反向重复区a（IRa）和反向重复区b（IRb）长度均为43 831 bp，GC 含量为36.00%。羊踯躅叶绿体共有146个基因。其中：92个基因为蛋白编码基因，占全部基因数量的63.01%；其次是tRNA 基因，数量为46，占全部基因数量的31.51%；剩余的均为rRNA 基因，数量为8，占全部基因数量的5.48%。

图1 羊踯躅叶绿体基因组图谱Fig.1 Chloroplast genome map of R.molle

由图1 可见：IR 区域有4个rRNA 基因，分别是rrn16S，rrn23S，rrn4.5S和rrn5S；IR 区域有9个tRNA 基因，分别是rnA-UGC、trnI-CAU、trnI-GAU、trnL-UAG、trnM-CAU、trnN-GUU、trnR-ACG、trnR-UCU、trnV-GAC；IR 区域还有4个蛋白质编码基因，分别是rpl32、rps15、rps16、ycf4。这17个基因在IRa 区域和IRb 区域均存在，且除trnI-CAU基因外均为2个拷贝，trnI-CAU基因在IRa 区域和IRb 区域各出现2 次，存在2个拷贝，所以trnI-CAU基因为4个拷贝。

将测序结果进一步分析，对其中已知功能的基因进行分类，结果见表1。由表1 可知，根据功能可将其分为4 大类，分别是光合作用类、自我复制类、其他基因类以及未知功能基因类。光合作用类基因，按功能又可进一步分为6 小类，分别是光系统I 的亚基、光系统Ⅱ的亚基、细胞色素b/f 复合体的亚基、NADH-脱氢酶的亚基、二磷酸核酮糖羧化酶大亚基和ATP 合酶亚基。在光合作用功能相关的基因中：种类最多的是光系统Ⅱ的亚基基因，有psbA、psbB、psbC等15 种；其次是NADH-脱氢酶的亚基基因，有ndhA、ndhB、ndhC等11 种；种类最少的是二磷酸核酮糖羧化酶大亚基基因，仅有rbcL。自我复制类基因可分为核糖体大亚基、核糖体小亚基、RNA 聚合酶亚基、rRNA 基因和tRNA 基因共5 小类，其中tRNA 基因数量最多，有46个基因。未知功能的基因有3个，分别是lhbA、ycf3、ycf4。

表1 羊踯躅叶绿体基因组编码的基因†Table 1 Genes present in chloroplast genome of R.molle

羊踯躅叶绿体具内含子的基因长度见表2。由表2 可知，将已知基因功能进行分类后，有15个基因具有内含子。就分布位置来说，处于LSC区的有11个，处于IR 区的有4个。就内含子长度来说：自我复制功能相关trnk-uuu基因的内含子Ⅰ的长度最长，达2 499 bp，其次是未知功能基因ycf3，其内含子Ⅰ的长度为1 683 bp，核糖体小亚基rps12不具备内含子；具有内含子Ⅱ的仅rps12，其长度为538 bp。就外显子长度来说，外显子Ⅰ的平均长度为144 bp，外显子Ⅱ的平均长度为248 bp，远低于内含子Ⅰ的平均长度（966 bp），仅rps12具有外显子Ⅲ，长度为114 bp。外显子Ⅰ长度最长的是NADH-脱氢酶的亚基中的ndhB，长度为721 bp，该基因的外显子Ⅱ长度也最长，为758 bp；外显子Ⅰ长度最短的是细胞色素b/f 复合体的亚基中的petB和petD，长度仅为6 bp，外显子Ⅱ长度最短的是核糖体小亚基rps12，长度为26 bp。

表2 羊踯躅叶绿体具内含子基因的长度Table 2 Gene length of introns in chloroplasts of R.molle bp

2.2 羊踯躅叶绿体基因组的密码子偏好性

羊踯躅叶绿体基因组中密码子组成见表3。由表3 可知，羊踯躅叶绿体基因组中共有22 312个密码子，基因编码区的序列占63.01%。叶绿体基因组中密码子组成共20 种氨基酸。其中：亮氨酸数量最多，有2 456个（占11.01%），属于编码率最高的氨基酸；数量最少的是半胱氨酸，仅260个（占1.17%），属于编码率最低的氨基酸。

进一步对密码子进行同义密码子相对使用度（RSCU）分析，结果见表3。由表3 可知，有32个密码子的RSCU 值大于1。其中：有16个氨基酸的密码子以U 结尾，14个氨基酸的密码子以A结尾，仅有2个氨基酸的密码子以G 结尾。可见，氨基酸的密码子以A 和U 结尾的较多，以C 和G 结尾的密码子的出现频率相对较低。同时，羊踯躅叶绿体中4个NCG 型（N 代表4 种碱基中的任一种）密码子的RSCU 值较低，UCG、GCG、CCG、ACG 的RSCU 值分别为0.477 6、0.468 4、0.468 0、0.368 8。CG 含量的降低有利于抑制mRNA的降解，从而使蛋白产量增加。4个NUA 型密码子具有较高的RSCU 值，UUA、GUA、AUA、CUA 的RSCU 值分别为2.152 2、1.400 8、0.910 2、0.786 6，说明羊踯躅可能属于DNA 甲基化程度较高的物种[13]。

表3 羊踯躅叶绿体基因组中密码子偏好性†Table 3 Relative synonymous codon usage(RSCU)in chloroplast genome of R.molle

2.3 羊踯躅叶绿体基因组的简单重复序列（SSR）

在羊踯躅叶绿体基因组中已注释的编码序列之间，存在大量的间隔序列，这里面包含了大量的重复序列，这些序列的组成影响着染色体的空间构象，也影响了编码序列的表达，因此研究重复序列对于分析基因调控网络的功能具有重要意义。羊踯躅叶绿体基因组的简单重复序列分析结果见表4。由表4 可知：在单核苷酸重复单元中以A/T 为主碱基，其数量占97.26%；在二核苷酸重复单元中，均为AT/AT，重复频率集中在5 次；三碱基的主要重复基序为AAG/CTT，其数量占简单重复序列位点重复单元总数量的11.48%，其次是AAT/ATT，占简单重复序列位点重复单元总数量的9.24%；四碱基至六碱基重复出现频率相对较低，均在0.6%以下。A/T 和AT/AT 占简单重复序列位点重复单元总数量的64.70%，可见羊踯躅叶绿体基因组富含AT，与叶绿体简单重复序列位点以AT 为主的理论吻合[14]。

表4 羊踯躅叶绿体基因组的简单重复序列Table 4 Simple sequences repeat(SSR)analysis of chloroplast of R.molle

简单重复序列的重复次数决定着重复碱基序列的长度，从而影响简单重复序列的多态性[15]。羊踯躅叶绿体组简单重复序列中从单碱基到六碱基各基元重复次数为5 ～10 的均有较多分布。单碱基重复次数为5 ～10 的较多，共179个，占简单重复序列位点重复单元总数量的50.14%；双碱基到六碱基重复次数为3 ～4 的较多，且随着重复次数的增加，其简单重复序列位点重复单元数量呈下降趋势。

2.4 羊踯躅叶绿体基因组的系统发育

基于叶绿体基因组构建的羊踯躅系统进化树如图2所示。由图2 可见，外群物种桤叶树科的城口桤叶树C.fargesii和云南桤叶树C.delavayi与杜鹃花科物种的亲缘关系较远，杜鹃花科中杜鹃属植物距离白珠树属和越橘属距离一致。进一步分析杜鹃属内物种的亲缘关系发现，映山红亚属Subgen.Tsutsusi映山红R.simsii、溪畔杜鹃与其他亚属杜鹃距离较远，杜鹃亚属Subgen.Rhododendron的照山白R.micranthum、着生杜鹃R.kawakamii、大天顶杜鹃R.datiandingense、秀雅杜鹃R.concinnum4个物种的聚类距离与映山红亚属的聚类距离相近，羊踯躅亚属Subgen.Pentanthera物种与常绿杜鹃亚属Subgen.Hymenanthes的灵宝杜鹃R.henanensesubsp.lingbaoense、马缨杜鹃R.delavayi、阔柄杜鹃R.platypodum3个物种的距离较近。进一步分析羊踯躅与常绿杜鹃亚属这3个物种的遗传距离发现，羊踯躅与灵宝杜鹃的遗传距离为6.621×10-3cM，与马缨杜鹃的遗传距离为8.411×10-3cM，与阔柄杜鹃的遗传距离为6.740×10-3cM。可见进化树的聚类结果与经典植物分类学研究结果是一致的，在植物分类学研究中羊踯躅是杜鹃属羊踯躅亚属五花药组的唯一物种，在进化树中羊踯躅单独成组，聚类分析及遗传距离分析结果表明，与羊踯躅亲缘关系最近的物种是杜鹃属常绿杜鹃亚属的灵宝杜鹃。

图2 基于叶绿体基因组构建的羊踯躅系统进化树Fig.2 Phylogenetic tree of R.molle based on chloroplast genome

3 结论与讨论

叶绿体是光合作用的关键场所，能够独立自我完成RNA 的复制及转录。羊踯躅通过光合作用，可获得其生长过程中必需的有机物及能量[16]。选择叶绿体基因组进行测序组装，来研究物种的亲缘关系和种群进化关系，是因为其有独特优势，测序组装相对容易些[17]：与核基因组相比，其基因组较小，长度一般为120 ～160 kbp；与线粒体基因相比，其基因组结构相对保守。本试验中对观赏及药用特性俱佳的羊踯躅的叶绿体基因组进行了测序，并获得了较好的组装和注释结果。

测序组装结果表明，羊踯躅叶绿体基因组具有高度保守的四段式结构，共注释有146个基因，其中大部分基因为蛋白编码基因，数量为92，tRNA 基因为46个，rRNA 基因为8个。羊踯躅叶绿体基因密码子偏好A 和U 结尾，单核苷酸重复基序主要为A/T（97.26%）。在所选19个杜鹃属物种和2个外群物种的系统进化树中，羊踯躅单独成组，与其关系最近的物种是灵宝杜鹃。

羊踯躅叶绿体基因组高度保守的四段式结构，与已经报道的爵床科黄猄草属菜头肾Championella sarcorrhiza[18]、漆树科盐肤木属盐肤木Rhus chinensis[19]等绿色植物物种类似，均有1个大单拷贝区、1个小单拷贝区和2个反向重复区。羊踯躅叶绿体基因组的结构功能与菜头肾、盐肤木等物种的基因组也高度相似，基因表达主要集中在光合作用和自我复制的相关功能，与叶绿体的主要功能一致，即将光能转变为化学能，将CO2和水转变为糖。

密码子在基因组与蛋白质的联系中发挥着重要作用。揭示密码子偏好性的影响因素，能为基因组研究及其遗传改良提供理论基础[20]。基因的碱基组成是密码子偏好性的普遍影响因素，牛元等[21]认为影响密码子偏好性形成的主要原因是选择和碱基突变，密码子的应用模式会影响基因的表达，进一步影响到密码子偏好性。叶绿体基因组中的密码子偏好性可揭示物种基因组的进化关系。羊踯躅叶绿体基因密码子偏好A 和U 结尾，同义密码子相对使用度分析结果也充分证明了这一观点，与对蝶形花科香槐属植物永椿香槐Cladrastis yunchunii[22]、马兜铃科细辛属辽细辛Asarum heterotropoides[23]、蔷薇科李属大山樱Prunus sargentii[24]的相关研究结果一致。

叶绿体DNA 中富含碱基A/T。羊踯躅叶绿体基因组的单核苷酸重复基序中A/T 占97.26%，椭圆叶花锚Halenia elliptica叶绿体基因组中单核苷酸重复序列中A/T 占94.53%[25]，菜头肾叶绿体基因组的单核苷酸中A/T 占93.48%[18]；羊踯躅叶绿体DNA 简单重复序列中A/T 和AT/AT 占简单重复序列位点重复单元总数量的64.71%，椭圆叶花锚和菜头肾叶绿体基因组的单核苷酸重复序列中其占比依次为61.90%、61.65%[18,25]。

在所选19个杜鹃属物种和2个外群物种的系统进化树中，羊踯躅单独成组，这一结果与经典植物分类学研究结果一致。羊踯躅与常绿杜鹃亚属物种的亲缘关系较近，与其关系最近的物种是灵宝杜鹃。在通过叶绿体基因序列构建得到的系统进化树中，杜鹃花科不同属、杜鹃属不同组的物种被区分开来，这一结果可以为新优杜鹃或者野外变异杜鹃的鉴定提供参考；根据叶绿体基因组的序列，可以定位其在系统发育进化树中的位置，分析与其亲缘关系较近的杜鹃种，可为杜鹃属植物的遗传多样性研究提供参考。

在研究过程中发现，羊踯躅野外分布地生境的破坏较为严重，羊踯躅野外群落较为稀少。国内羊踯躅的引种繁育主要依赖科研院所和花卉苗木培育企业。因地理位置、环境、栽培方式和管护的不同，羊踯躅产生了不同程度的变异，但目前仅选育出1个变异品种‘金踯躅’R.molle‘Jin Zhizhu’。目前，叶绿体基因测序组装技术存在一定局限性：成本虽然有所降低，测序费用仍然偏高，难以使用该技术对经济价值相对较低的物种大规模开展研究工作；基因组序列组装虽然不断接近完成，但依然存在较多空白，而且基因组中的高重复、复杂区域依旧是基因组组装面临的重要问题；因其错误率高，需要得到大量测序数据进行纠正；基因组数据研究是其他功能研究的基础，但表型性状的遗传机理十分复杂，如何有效开展多组学研究也是今后面临的主要问题。为给羊踯躅的衍生新品种杜鹃的初期新品种鉴定提供参考，下一步将利用分子标记技术研究羊踯躅杜鹃花种质的遗传多样性，从基因组水平揭示其遗传变异程度。