须弥葛叶绿体基因组密码子使用偏好性分析

万美暄，黄显俊，李 雪，徐 立，*

(1.大理大学基础医学院，云南大理 671000；2.大理大学药学院，云南省高校滇西道地药材资源开发重点实验室，云南大理 671000)

叶绿体(chloroplast，Cp)是半自主性细胞器，在植物光合作用、代谢反应、脂肪酸、氨基酸、维生素等生物合成中至关重要[1]。叶绿体DNA一般以双链环状存在于高等植物叶绿体中，长度在120～180 kb，具有典型的四分体结构[2]，且由于其母系遗传和适度的进化速度，在序列和结构上高度保守[3]。因此，Cp基因组测序备受关注，NCBI数据库也提供了越来越多的叶绿体基因组，这为种族遗传学和系统发育方面的研究提供了方便。

密码子(condon)是信使RNA(mRNA)翻译的基本信息单位，61个密码子编码20种氨基酸，除色氨酸和蛋氨酸，其他氨基酸均有1种以上的密码子编码[4]。然而，在蛋白质合成中同义密码子(Synonymous condon)不一定以相同的频率被使用，这种特定密码子使用频率高于其他同义密码子的现象被称为“同义密码子使用偏好性(codon usage bias，CUB)”[5-6]。CUB广泛存在于生物体中，反映了不同物种或品种之间基因的起源、进化趋势和突变模式[7]。有研究表明，CUB与基因表达水平之间存在相关性[8]。此外，密码子使用模式可能会影响 mRNA 生物合成、蛋白质折叠、翻译延伸率及其他下游表达的生物学功能[9]。因此，对CUB分析可以更好地了解物种的遗传规律、蛋白质的机制及预测新基因。

须弥葛(Haymondiawallihii)，俗称“马鹿花”，为豆科蝶形花亚科(Papilionaceae)须弥葛属(Haymondia)植物[10]。在我国主要分布于广州、四川、广西、福建、云南等省区及云南省内的思茅、红河、临沧、保山、德宏 、大理 、楚雄、怒江等州市[11]。须弥葛花较大，生长速度快，可做观赏绿化，在大理白族自治州常被采集花朵，腌制野菜[12]。有研究发现，须弥葛中含有大豆素、大豆苷、葛根素、豆甾醇等化合物[12]，具有抗心血管疾病的生理活性，还能缓解高血压患者的头痛，解决平滑肌痉挛，镇咳及抑制癌细胞生长等药性[13]。笔者基于须弥葛植物叶绿体基因组高通量测序结果，挑选出长度>300 bp的52条蛋白编码序列，对其进行整合分析，探讨须弥葛密码子使用偏性的影响因素，并最终确定最优密码子，以期为须弥葛的进化关系研究、品种改良及叶绿体基因组的应用提供参考。

1 材料与方法

1.1 试验材料的获取须弥葛植株的新鲜幼嫩叶片于2019年12月采自云南省丽江永胜县( 25°53′73.23″E，100°25′48.40″N)，经云南省高校滇西道地药材资源开发重点实验室的夏从龙教授鉴定为豆科须弥葛属植物须弥葛(Haymndiawallichii)。用无菌水将采集到的健康须弥葛叶片冲洗多次，装入有变色硅胶的自封袋中脱水干燥，委托上海生物科技有限公司的Illumina No-vaSeq 6000平台测序，最后通过组装和注释获得须弥葛叶绿体全基因组序列号(NCBI数据库，登录号为NC052910.1)。

1.2 叶绿体基因组蛋白质编码序列的获取从NCBI数据库下载须弥葛叶绿体基因组蛋白质编码序列(CDS)后，为加强数据可靠性及减少数据冗杂，剔除重复基因和长度<300 bp的基因[14]。最后将符合条件的52条CDS整合成一个fasta文件，用于后续分析。

1.3 密码子使用偏好性分析将整理的fasta文件通过CodonW1.4.2软件处理[15]，获取到须弥葛叶绿体基因组蛋白质编码基因序列的相对同义密码子使用度(relative synonymous codon usage，RSCU)、编码氨基酸的密码子数量、有效密码子数(effective number of codon，ENC)、最优密码子使用频率(frequency of optical codons，Fop)、密码子适应指数(codon adaptation index，CAI)、密码子偏性指数(codon bias index，CBI)。然后，将52条CDS依次提取整理成52个fasta文件，使用在线程序CUSP[16]计算每个基因不同位置的密码子GC含量(GCall、GC1、GC2、 GC3)，将所得数据使用SPSS 24.0和Microsoft Excel软件进行统计及Person相关性分析。

1.4 PR2 plot绘图分析PR2 plot图是以C=G和A=U为矢量发出中心点，用A3/(A3+T3)为纵坐标，G3/(G3+C3)为横坐标绘制而成，表明影响密码子偏好性的主要因素[17]。

1.5 中性绘图分析将GC3作为横坐标，GG1和GC2的平均值(GC12)为纵坐标，作二维散点图，再次绘制回归曲线，然后对GC3和GC12进行相关性分析。

1.6 ENC-plot绘图分析首先利用Excel绘制GC3和ENC 的二维散点图，并根据公式ENC= 2+GC3+29/[GC32+(1-GC3)2]绘制标准曲线，以此判断基因点在图中分布情况。然后利用ENC公式求得ENC期望，再通过 ENC比值=(ENC期望-ENC实际)/ENC期望，算出ENC比值频数分布，从而确定影响密码子偏好性的主要因素。

1.7 最优密码子的筛选以52条CDS分别得到的ENC值为基础，用Excel将 ENC值依次排序，选出ENC值最高的5个基因和最低的5个基因建立高低表达基因库，然后将这些基因序列整理成不同的fasta文件，运行Codon W 1.4.2软件得出2个库中编码各氨基酸密码子的RSCU值，计算△RSCU值(RSCU高表达，-RSCU低表达)。以高表达组RSCU>1，且△RSCU≥0.08的密码子为最优密码子[18]。

1.8 统计分析使用SPSS 24.0和Excel软件，进行Pearson相关分析，显著性水平为P<0.05或P<0.01。

2 结果与分析

2.1 须弥葛叶绿体基因组密码子偏性分析

2.1.1相对同义密码子使用度。须弥葛叶绿体基因组中的52条CDS序列的相对同义密码子使用度(relative synonymous codon usage，RSCU)结果显示(表1)，31个高频密码子(RSCU>1)中以U结尾的有16个，以A结尾的有14个，以G结尾的有1个，表明须弥葛叶绿体基因组密码子更倾向A/U结尾。所有高频密码子RSCU值均在1～2，除亮氨酸(Leu)密码子UUA>2，提示须弥葛叶绿体基因组密码子使用偏好性较弱。

表1 须弥葛叶绿体基因组中各氨基酸密码子的RSCU值

2.1.2基因组密码子碱基组成及偏性。利用软件CodonW1.4.2和CUSP分析符合条件的52条CDS序列(表2)，结果表明，52条CDS序列的总GC含量(GCall)为38.10%(范围为29.27%～46.56%)，第1位(GC1)为46.89%(范围为31.99%～59.31%)，第2位(GC2)为39.48%(范围为27.62%～54.68%)，第3位(GC3)为27.87%范围(为21.15%～36.88%)。GC2与GCall相近，GC1和GC3与GCall相差较大，表现为GC1>GC2>GC3。进一步说明须弥葛叶绿体基因组密码子偏好A/U结尾的密码子。

表2 密码子的碱基组成及偏性相关指标

52条CDS的CAI(codon adaptation index)为0.166(范围为0.110～0.294)，CBI(codon bias index)为-0.104(范围为-0.212～0.175)，FOP(frequency of optical codons)为0.349(范围为0.269～0.520)，ENC为47.54(范围为38.66～59.33)。以上指标进一步提示须弥葛叶绿体基因组密码子的使用偏性和基因表达均较弱。

ENC值与各位置GC含量的相关性分析显示(表3)，GCall与GC1、GC2和GC3含量的相关性均达到极显著相关水平(P<0.01)；GC1与GC2含量的相关性达到极显著水平(P<0.01)；GC3与GC1及GC2含量的相关性达到显著水平(P<0.05)；ENC与GC1、GC2含量的相关性均不显著(P>0.05)，但与GC3含量的相关性极显著(P<0.01)，说明第3位GC含量对密码子的使用偏性有一定程度的影响。

表3 各基因参数的相关性分析

2.2 影响密码子偏性的因素分析

2.2.1PR2-plot分析。 PR2绘图分析(图1)，所有点分布在G3/(G3+C3)上的0.26～0.67和A3/(A3+T3)上的0.38～0.57，但各基因在4个象限中明显分布不均匀。若密码子使用偏好完全受碱基突变的影响，大多数点应位于0.5水平中心线或附近，但该结果却不同。因此，就碱基频率而言，G>C，T(U)>A。由此推测，自然选择和突变都是影响须弥葛叶绿体基因组密码子使用偏性的因素。

图1 须弥葛叶绿体基因组PR2-plot分析

2.2.2中性绘图分析。以中性绘图分析GC12和GC3的关系(图2)，GC12的取值0.340 1～0.557 5，GC3的取值0.211 5～0.368 8，所有基因位于对角线以上(除Ycf2基因)。Pearson相关分析发现，GC12和GC3之间的相关性较弱(R2=0.125 8，r=0.355 0)，相关系数r值距离1较远。该结果提示自然选择是须弥葛叶绿体基因组密码子使用偏性的主要影响因素。

图2 须弥葛叶绿体基因组中性绘图

2.2.3ENC-plot分析。ENC-plot绘图分析(图3)，大部分基因散乱地分布在标准曲线下方，只有少数基因分布在曲线附近或曲线上。另外，52条CDS的ENC比值在-0.1～0.3(表4)，其中，25个基因的ENC值和预期接近，分布在-0.050～0.050；而其他27个基因的ENC值不在该范围内。若碱基突变是影响密码子的主要因素，大部分基因应在曲线上或曲线附近，且ENC比值在-0.05～0.05，但结果却相反。提示碱基突变虽是影响须弥葛密码子偏好性的因素，但受自然选择更多，这与中性绘图分析结果相一致。

图3 须弥葛叶绿体基因组ENC-plot分析

表4 ENC比值频率分布

2.3 最优密码子确定该研究将ENC值最高的5个基因(ycf3、clpP、ycf2、rpl2、rpoA)和 ENC值最低的5个基因(rpll6、aptF、petD、psbA、rps12)分别作为高低基因表达组，建立高低表达基因库(表5)，计算得到△RSCU。RSCU>1为高频密码子，△RSCU≥0.08的密码子为高表达优越密码子，选取二者共有的密码子为最优密码子，据此符合条件的有10个密码子，其中8个以A/U结尾。

表5 须弥葛叶绿体基因组最优密码子分析

3 讨论

密码子使用偏好广泛存在于多种生物体内，可能受 GC 含量、基因长度、基因功能、选择、突变等多种因素影响，是一种复杂而又不可避免的现象[19]。以往大多利用全基因组信息研究物种的密码子偏好[20]。随着近年来新一代测序技术的飞速发展，出现了大量的基因组和转录组数据，Machado等[21]通过使用深度基因组群体测序推断黑腹果蝇同义位点的选择。Chu 等[22]利用 Ribo-seq 和 RNA-seq 方法研究了密码子使用偏好如何影响拟南芥的翻译模式。Guan等[23]通过对菲牛蛭密码子使用模式的分析，发现基因进化的基本要素是突变压力和选择。该研究通过PR2绘图、中性绘图、ENC-pltot等对影响须弥葛密码子偏好性的因素进行了分析整理，并最终确定最优密码子。

由于密码子的简并性，氨基酸的类型虽不会因GC3同义突变所改变，但其仍被认为是影响氨基酸类型的一个重要因素[24]。该研究发现，须弥葛叶绿体基因组总GC含量(GCall)是38.10%，且GC1>GC2>GC3，表明密码子不同位置GC含量不同，密码子偏好使用以A/U结尾的密码子，这与文心花、中华绒螯蟹、豆蔻等植物偏好性一致[16，25-26]。RSCU是确定密码子使用偏好性的主要指标，若RSCU<1，则表明该密码子使用频率偏低；RSCU>1，则说明该密码子使用频率偏高；当RSCU=1时，该密码子使用无偏好性[27]。该研究发现，须弥葛叶绿体基因组高频密码子中RSCU>1的有18个，大部分是以U和A结尾，G和C结尾的仅2个。最终确定最优密码子是10个，且8个是以A/U结尾，这与GC含量分析结果相同。

ENC通常表示密码子偏性的强弱，其ENC>45表示密码子偏好性较弱[28]。须弥葛ENC平均值为47.54，则说明密码子使用偏好性较弱。PR2-plot、中性绘图、ENC-plot分析结果显示，自然选择和突变在须弥葛叶绿体密码子偏好形成中都起作用，但前者作用更大。这与“玉铃铛”枣、翠雀叶绿体基因组中的密码子偏性影响结果一致[29]，但与蒙古韭不同[30]。由此可知，造成植物密码子偏好的因素存在差异，不同因素在不同物种间发挥作用也不同。

该研究通过生物信息学的方法介绍了影响须弥葛叶绿体基因组密码子使用偏好和模式的因素，并筛选出最优密码子，这对须弥葛的进一步开发运用有着重要意义。