森林茶园古茶树大理茶叶绿体基因组密码子偏好性及系统发育研究

佟岩，黄荟,2，王雨华*

1. 中国科学院昆明植物研究所，资源植物与生物技术重点实验室，云南省野生资源植物研发重点实验室，云南 昆明 650201；2. 怀化学院生物与食品工程学院，湖南 怀化 418000

古茶树多指百年左右、树体高大的自然生长或栽培的茶树，主要为茶组植物的野生型、半野生型或栽培型等各类茶树资源，云南是世界上古茶树和野生大茶树资源最为丰富的地方，至今仍保留有较大面积的古茶林[1]。古茶树是重要的茶树种质资源基因库，在研究茶树的起源、进化、生物多样性、品种改良等方面都具有重要的价值[2]。裴盛基等[3]从民族植物学视角提出，以产出生态茶叶为目标，按传统方式管理，符合绿色、有机种植标准，面积在0.5 hm2以上，伴生有高于5 m 以上至少1 种形成林冠的其他野生或栽培的木本植物，覆盖度在 10%以上的栽培山茶科山茶属茶组植物群落称为森林茶园，一般古茶林都能达到森林茶园的标准，称为古茶林森林茶园。古茶林森林茶园是我国二三千年的茶叶栽培史中，劳动人民通过各种传统农业耕作模式遗留下来的，既保留了天然林的主要特征和功能，又能在人为调控作用下获得最大利用价值的茶园模式，具有重要的生态、经济和文化意义[3-4]。茶组植物大理茶（Camellia taliensis）是云南古茶树资源中分布广泛、面积大、适应性强的一种茶树资源，也是古茶林森林茶园的重要物种，其群落为物种多样性丰富的原生自然植被，调查发现白莺山古茶园、镇源千家寨等古茶林森林茶园的主要建群物种即为大理茶，其是典型的森林茶园古茶树物种[5-6]。前人基于分子标记研究发现，大叶茶可能是由栽培的大理茶与小叶茶进行杂交后再驯化而来[7-8]。此外，大理茶也是迄今发现与栽培茶化学成分最为相近，以及富含咖啡碱和茶多酚的近缘茶组植物，自古以来就被当地少数民族利用作为茶饮，有名的“大理感通茶”“银生茶”“野生乔木茶”等均来源于大理茶[5,9]。

大理茶分布于云南西南部山区，横断山脉澜沧江中段至伊洛瓦底江流域，海拔1 500～2 400 m的南亚热带山地的中山地带，野生大理茶为国家二级保护树种，群落物种多样性极为丰富[10]。茶叶品质的季节性变化及地区性差异主要是由于光、温度等环境因子对茶树多酚类、咖啡碱及茶氨酸等主要成分合成与代谢的影响导致[11]。在茶树种植中多用遮阴覆盖等农业措施调节茶树的光合作用，从而提升茶品质[11-12]。叶绿体是植物体中进行光合作用的重要器官，拥有独立于细胞核的半自主遗传系统[13]。叶绿体基因组DNA（Chloroplast genome DNA，cpDNA）具有分子量小、单倍性、单亲遗传的特点，其基因含量和基因组高度保守，基因类型与位置分布较为保守、遗传特性稳定，有些进化速率较快的基因如matK（Maturase kinase）等兼具较高的基因同源性和大量的变异位点，多态性高，在物种鉴定、系统发育关系和物种起源等方面应用广泛[13-14]。matK基因位于叶绿体trnK基因的内含子中，编码一种参与叶绿体Ⅱ型内含子剪切的成熟酶K，是叶绿体行使正常功能的关键基因之一[15]。聂传朋等[16]对茶树DNA 条形码引物进行了初步筛选，认为matK、rcbL可作为茶树品种间分类的备选序列。目前，matK作为植物DNA 条形码的标准序列之一，不仅能够用于植物种及变种之间的鉴别，还能够有效地对属间及属以上级别进行鉴定，对近缘种也有较好的分辨力。

遗传信息从mRNA 翻译成蛋白质时，由于密码子的简并性而对同义密码子的使用频率出现差异的现象称为密码子使用偏好性（Codon usage bias，CUB）[17]。不同物种、不同基因的密码子偏好性不同，是物种长期进化过程中受到自然选择、突变压力和随机遗传漂变等多种因素影响而形成的，密码子偏好性的存在可降低外源基因在宿主中的表达量，从而维持物种的稳定遗传[18-19]。本研究通过对大理茶叶绿体基因组的密码子使用与基因的表达情况、CDS 区的碱基组成等密码子使用偏好特征及其偏好性的影响因素进行系统分析，筛选最优密码子，为提高大理茶外源基因的表达转化效率提供基础。

1 材料与方法

1.1 大理茶叶绿体基因组序列的获取

大理茶（C. taliensis）叶绿体基因组参照Yang 等[20]报道，基因组序列（NC_022264.1）及注释文件均从NCBI 数据库中下载，获得的大理茶叶绿体基因组序列长度为156 974 bp，共包含 89 条编码序列（Coding sequence，CDS）。移除基因编码序列中长度小于300 bp的序列与基因重复序列，最终获得54 条符合条件的大理茶CDS 序列用于后续分析。

同时从 NCBI 数据库中下载小叶茶（C.sinensisvar.sinensis，NC_020019.1）、大叶茶云抗 10 号（C. sinensisvar.assamicacv.Yunkang10，MH019307.1）、白毛茶（C. sinensisvar.pubilimba，NC_038198.1）、德宏茶（C.sinensisvar.dehungensis，KJ806279.1）、大苞茶（C. grandibracteata，NC_024659.1）、秃房茶（C. gymnogyna，NC_039626.1）、大厂茶（C.tachangensiscv.Xingyi6，NC_056149.1）、毛叶茶（C. ptilophylla，NC_038198.1）、突肋茶（C. costata，NC_061904.1）9 个茶组植物，以及木荷（Schima superba，NC_035545.1）和紫茎（Stewartia sinensis，NC_035640.1）2 个外类群山茶科木荷属、紫茎属植物的叶绿体基因组序列。采用上述方法分别获取这11 个物种的符合条件的CDS 序列用于后续茶组植物系统发育分析。

1.2 密码子偏好性相关参数计算

用CodonW v1.4.2 软件对大理茶叶绿体基因组54 条CDS 序列的编码氨基酸长度（Laa）、有效密码子数（Effective number of codon，ENC）、相对同义密码子使用度（Relative synonymous codon usage，RSCU）、密码子适应指数（Codon adaptation index，CAI）、同义密码子中GC 含量和密码子第3 位上A、T、C、G 的含量（分别以A3s、T3s、C3s、G3s表示）进行分析。再使用CUSP（http://imed.med.ucm.es/EMBOSS）在线软件计算 GC 总含量（GCall）及密码子第1、第2、第3 位碱基组成中的GC 含量（分别以GC1、GC2和GC3表示），去除编码蛋氨酸（Met）、色氨酸（Trp）单密码子和3 个终止密码子（UAA、UAG、UGA）。用R 语言的ggplot2 软件包基于筛选到的59 个密码子的RSCU 值绘制堆积柱形图（Stacked barplot）展示密码子偏向性。利用SPSS 18.0 软件的Pearson 相关系数法对大理茶叶绿体基因组密码子偏好性相关参数进行相关性分析。

1.3 中性绘图分析

计算大理茶叶绿体各CDS 序列中密码子GC1和GC2的平均值GC12，以GC3和GC12为横纵坐标绘图散点图，分析GC3和GC12的相关性。若GC3与GC12显著相关，则表明3 个位置上碱基组成差异不显著，突变是密码子使用的主要影响因素；若相关性不显著，则表明密码子使用的偏好性受自然选择影响较大[21]。

1.4 ENC-plot 分析

以GC3为横坐标，ENC 为纵坐标，绘制ENC-plot 散点图，同时绘制标准曲线（标准曲线绘制公式：ENC期望=2+GC3s+29/[GC3s2+(1-GC3s)2])[22]。通过基因散点分布的位置判断影响密码子偏好性的原因，若基因分布在标准曲线上或附近，则密码子偏好性主要受突变影响；若基因位于标准曲线下方较远的位置，则密码子偏好性受到自然选择和其他因素的影响较大[23]。

1.5 PR2-plot 分析

以G3/(G3+C3)为横坐标，A3/(A3+T3)值为纵坐标，进行PR2-plot 绘图，中心点为A=T且C=G 时的值，表示密码子无使用偏好性，只受到突变的影响，其余点和中心点的距离则表示偏倚的程度和方向[24]。

1.6 最优密码子分析

根据 ENC 大小进行排序，选取前后各10%的基因构建高偏好性和低偏好性库。计算2 个库各密码子的RSCU 值和ΔRSCU 值，去除蛋氨酸、色氨酸及终止密码子UAA、UAG、UGA，将同时满足RSCU＞1 和ΔRSCU＞0.08的密码子确定为最优密码子[25]。

1.7 系统发育分析

利用10 个茶组植物和2 个山茶科外类群物种的叶绿体基因组CDS 和matK基因分别进行系统进化分析。用CLUSTALW 软件对所有序列进行多序列的比对，利用MEGA-X 软件的系统邻接法（Nerghbor-joining，NJ），设置1 000 次的Bootstrap 检验各分支的置信度，其他参数均为默认值，对12 个物种构建系统进化树，用EvolView（http://evolgenius.info/#/）软件对进化树进行美化作图。

2 结果与分析

2.1 大理茶叶绿体密码子相关偏性指数分析

大理茶叶绿体基因组全长156 974 bp，去除小于300 bp 的CDS 及重复序列后，共筛选出54 条蛋白编码基因，全长64 809 bp，占基因组全长的41%，基因编码氨基酸序列长度范围为100～2 282，平均长度为399。对这些基因的密码子组成和偏好性进行统计分析，结果如表1 所示，蛋白编码序列平均GC 含量为37.68%，其中GC1＞GC2＞GC3，显示GC 在密码子3 个位置上的分布并不均匀，GC 含量均小于50%，表明大理茶叶绿体基因组密码子偏好使用A/U 碱基，且偏好以A/U 碱基结尾。

表1 大理茶叶绿体基因密码子的GC 含量及ENC 值Table 1 GC contents and ENC values in chloroplast genome of C. taliensis

有效密码子数用于揭示密码子的偏好程度，ENC 的取值范围为21～60，ENC 值的大小与密码子偏好程度呈现负相关，当ENC 值小于45 时，意味着密码子存在偏好性。大理茶叶绿体基因组的ENC 值介于35.64～56.67，平均值为46.1，ENC 值大于45 的基因有36条，表明大理茶叶绿体大部分蛋白编码基因序列密码子的使用偏好性较弱。

大理茶叶绿体基因组密码子各参数之间的相关性分析结果显示（表2），GC1与GC2呈极显著相关（P＜0.01），但GC1与GC3及GC2与GC3之间均无相关性，表明密码子第1、2 位碱基组成情况高度相似，而第3 位上的碱基组成随机性比较大，与1、2 位碱基组成有较大差异。ENC 与GC3呈显著正相关（P＜0.05），表明密码子第3 位上碱基组成的变化对密码子的使用偏好性影响最大。

表2 密码子主要参数的相关性分析Table 2 Correlation analysis between the indexes of codon use

2.2 整体密码子偏好分析

RSCU 指某一特定的密码子在其同义密码子中出现的概率，可直观反映密码子的使用情况。结果显示（图1），大理茶叶绿体基因组中共有33 个高频密码子，其中以A、U、C、G 结尾的分别有14、16、0、3 个。此外，大理茶叶绿体基因组中还出现了偏好性极强的密码子（RSCU 值≥2），为编码亮氨酸（Leu）的UUA。说明大理茶叶绿体基因组偏好A/U结尾的密码子。

图1 大理茶叶绿体基因组CDS 密码子RSCU 直方图Fig. 1 Analysis of relative synonymous codon usage (RSCU) in the chloroplast genome of C. taliensis

2.3 最优密码子的确定

根据ENC 值的大小，选出5 个高表达基因和5 个低表达基因建立基因库。通过计算得出ΔRSCU 值＞0.08 的密码子有28 个（表3），为高表达优越密码子。将同时满足RSCU＞1、ΔRSCU＞0.08 的密码子确定为最优密码子，最终得到大理茶叶绿体11 个最优密码子，分别为AUU、AUC、UCA、CCA、ACA、GCA、AAU、GAU、AGA、GGA 和GGG，其中2个以C、G 结尾，其余均为A、U 结尾。

表3 大理茶叶绿体基因组最优密码子确定Table 3 Preferred codon in chloroplast genome of C. taliensis

2.4 中性绘图分析

大理茶叶绿体基因组54 个基因的中性绘图分析结果显示（图2），GC3和GC12的相关系数为0.007 8，回归系数为0.116 8，表示突变压力的贡献率为11.68%，而自然选择压力的贡献率为88.32%，说明GC3与GC12之间的相关性不强，突变压力对密码子偏好性的影响不大，表明大理茶叶绿体基因组密码子同时受到自然选择和突变压力的影响，但自然选择起到了更大的作用。

图2 大理茶叶绿体基因组中性绘图分析Fig. 2 Neutrality plot analysis of C. taliensis chloroplast genome

2.5 ENC-plot 分析

ENC-plot 分析结果显示（图3），大多数大理茶叶绿体基因都偏离了基因预测值的标准曲线，表明除由碱基突变引发的差异影响该基因组密码子偏性外，影响更多来自选择效应。ENC 值相对较小的基因位于预测曲线之下，并偏离曲线，说明这种类型偏好使用特定密码子，偏好性也较大，这部分基因的密码子偏好性更多受到自然选择等的影响。为检验基因与标准曲线的距离，统计了各基因ENC 与期望值的比值，结果显示（表4），55.6%的基因ENC 比值分布在 0～0.1 范围，7.4%分布在 0.1～0.2 及0.2～0.3 范围，表明ENC 的实际值与期望值相差较大。进一步说明其整体密码子偏好性主要受自然选择的正向影响，这种影响使基因拥有更低的ENC 值，从而获得更高的表达效率。

图3 大理茶叶绿体基因组ENC-plot 绘图分析Fig. 3 ENC-plot analysis of C. taliensis chloroplast genome

表4 ENC 比值频数分布Table 4 Distribution of ENC ration

2.6 PR2-plot 分析

通过绘制PR2-plot 图进一步分析大理茶叶绿体基因组密码子偏好性的影响因素。由图4 可知，大理茶叶绿体基因组编码基因并不是均匀地分布于4 个区域中，在右上区域分布很少。如果密码子偏好性仅受基因突变的影响，在PR2-plot 绘图分析中，4 种碱基的使用频率将相等。PR2-plot 分析显示大理茶叶绿体基因组种大部分基因受到自然选择或其他因素的影响。

图4 大理茶叶绿体基因组PR2-plot 绘图分析Fig. 4 PR2 bias plot analysis of C. taliensis chloroplast genome

2.7 系统进化分析

茶组植物是异花授粉植物，种间甚至组间的杂交比较普遍，经过长期的自然杂交和人工选育，其种间的界限比较模糊，传统的形态分类基于心皮或子房室的数目，花梗、萼片、花瓣和果实的发育状况及叶片特征等，以及细胞学、化学和数量分类学等，很难得到公认的茶组植物分类系统[26-32]。目前，基于测序技术、分子生物学和分子系统学，在DNA 水平上利用合适的分子标记解决物种间系统关系及分类是有效的方法[33-41]。

对10 种茶组植物和2 种山茶科植物共12种植物分别基于叶绿体基因组matK基因序列和蛋白编码序列构建系统发育树（图5），结果显示，2 种系统发育的拓扑结构相似，均将白毛茶、小叶茶、大苞茶、德宏茶和大叶茶云抗10 号聚为一大类群，大理茶和秃房茶聚为一支，白毛茶和小叶茶聚为一支，两个外类群物种木荷和紫茎聚为一支，突肋茶和大厂茶分别位于单独的分支。

图5 大理茶及茶组植物叶绿体系统发育树Fig. 5 Phylogenetic tree of C. taliensis chloroplast genome

3 讨论

密码子受自然选择、基因突变、蛋白质结构等多因素的共同影响逐渐形成一定的偏好性，与基因的GC 含量、长度、翻译效率、表达水平及tRNA 丰度等因素有关[25,42-45]。本研究中大理茶叶绿体基因组编码基因的 GC 含量为37.68%，密码子的GC1、GC2和GC3含量均小于50%，偏好使用A/U 结尾的密码子和碱基，与前人对普通油茶（C. oleifera）[46]、鸟王茶（C. sinensisvar.niaowangensis）[47]、香花油茶（C. osmantha）[48]的密码子使用偏性研究相一致。大部分大理茶基因ENC 值偏大，平均ENC 值为46.1，说明大理茶叶绿体基因密码子使用偏性较弱，暗示大部分大理茶叶绿体基因在序列水平上具有较低的表达潜力。

植物基因组中存在大量密码子偏好性现象是基因组进化的一个重要特征，是物种在长期的环境适应过程中形成的[49]。大理茶中11个最优密码子大部分以A/T 结尾，与大部分双子叶植物一致。ENC-plot、中性绘图和PR2-plot 等密码子使用性影响因素分析表明，大理茶叶绿体基因密码子偏性受到较小的碱基突变影响，更多受到自然选择等因素的影响。大理茶的野生、半野生和栽培型大都长期生长于生态环境稳定的森林环境中，加之云南当地少数民族生态化的栽培和管理模式，使得人工选择、基因突变等对大理茶叶绿体基因组密码子偏好性的影响较小，且密码子偏性较弱。而在人工栽培驯化的普通油茶叶绿体基因组密码子偏好性中，显示出突变和选择的共同作用[46]。

现代茶组植物的分类系统中影响较大的主要有Sealy 分类系统、张宏达分类系统、闵天禄分类系统等[26-29,50]，这些分类系统对大理茶物种的认定都比较清晰。陈亮等[30]根据形态学演化推断茶组可能有2 条演化线路，从原始茶亚属—大理茶—大厂茶—秃房茶的演化或是从原始茶亚属植物—大理茶—普洱茶、白毛茶和茶。黄荟等[33]研究显示，小叶茶与德宏茶聚为一支，大叶茶与大苞茶聚为一支，大理茶单独聚为一支，表明茶组植物基于叶绿体基因组的分类系统并未以子房室的数量进行聚类。本结果显示，5 室子房的大理茶与3 室子房的秃房茶聚为一支，与Li 等[51]基于全叶绿体基因组构建的ML 树一致。此外，根据赵东伟等[52]对大苞茶的研究发现，大苞茶没有独立的分布区，总是与大理茶和茶相伴出现，其形态性状也处于大理茶和茶之间的过度形态。本研究基于叶绿体matK基因序列和CDS的系统发育树显示大苞茶与茶的关系更近。

大理茶是茶组植物中较为原始的物种，是茶的重要野生近缘种，是具有重要经济价值和研究价值的古树茶资源。云南当地少数民族对大理茶资源的开发利用和保护历史悠久，至今仍保有完整的群落。本研究结果为大理茶重要农艺性状改良、外源基因密码子优化、提高基因表达效率、分子育种等研究提供了科学依据。此外，随着基因组测序的发展，结合基因组、细胞器基因组深入研究密码子使用规律，并结合形态、细胞、化学等研究共同揭示大理茶的起源、进化，以及大理茶如何参与大叶茶的起源进化也是后续研究的重要方向之一。