伯乐树叶绿体密码子役佣偏好分析

李佳灵，尹为治，方正，黄良鸿

（1.海南热带雨林国家公园管理局五指山分局，海南 五指山 572299；2.海南大学生态与环境学院，海南海口 570228）

“名不正则言不顺，言不顺则事不成”。本文所用基本术语均源自英文，学界用词尚未统一，造成较大混乱。为便于论述，特定义解释如下。

密码子（codon）：DNA 或mRNA 上三个相邻的碱基。共64 种，其中61 种编码氨基酸，其余3种为终止密码子。

同义密码子（synonymous codon），简称同义码：编码同一种氨基酸的多个密码子。生物体中使用的氨基酸有20 种，其中2 种各拥有1 种密码子，其余18 种各拥有2～6 种同义密码子。

密码子役佣（codon usage），简称码役：密码子被基因序列役使征用以编码氨基酸，多被译为“密码子使用”。事实上密码子是被使用，不是密码子使用什么。因此，“役佣”一词较“使用”更能精确表达“usage”原意。

同义密码子役佣偏好（synonymous codon usage bias），简称码役偏好。同义密码子非均衡役佣现象，基因偏好役用同义密码子中的某种或某几种。

同义密码子氨基酸族（synonymous family），简称码族：拥有相同数目同义密码子的氨基酸家族。一码族含2 个成员，各自拥有1 种密码子；二码族含9 个成员，各自拥有2 种同义密码子；三、四、六码族分别含1、5、3 个成员，各自分别拥有3、4、6 种同义密码子。

最优密码子（optimal codon）：在基因组多数基因中役用频率超过平均水平且在高表达基因中役用频率较高的一组同义密码子。最优密码子的筛选标准随研究目的、对象而异。

现役密码子，简称现役码：出现在基因序列，役佣于基因的密码子。

在编氨基酸：被基因编码，其密码子役佣于基因的氨基酸。

码役偏好与基因组和蛋白质相关联，生物学意义巨大[1,2]。目前，所有被研究过的物种都有着不同程度的码役偏好[3]。码役偏好受自然选择、碱基突变、基因漂变等因素共同作用[4]，还受到基因组大小[5]、tRNA 丰度[6]等因素影响。了解不同物种的码役模式，分析其影响因素，可以推断未知基因的表达，预测未知基因的功能[1]。

伯乐树Bretschneidera sinensis 系第三纪孑遗植物，单种成科，我国特有种，1999 年被列为国家Ⅰ级重点保护植物[7]。伯乐树叶绿体已经完成测序[8]，尚未有其码役等方面的研究报道。在已有的研究基础上，通过对伯乐树叶绿体基因码役偏好的分析，揭示其主要影响因素，确定最优密码子，本研究旨在为伯乐树叶绿体基因组学提供参考。

1 材料与方法

1.1 数据来源

从NCBI 数据库（https://www.ncbi.nlm.nih.gov/）中下载伯乐树叶绿体全基因组序列信息，登录号为NC_037753.1，共下载87 条基因编码序列，排除重复或较短的序列，获得53 个基因作为分析样本。样本基因均以ATG 为起始密码子，以TAG、TGA 或TAA 为终止密码子，且编码区长度大于300 bp。统计分析使用SPSS 22 及R 语言完成。

1.2 研究方法

1.2.1 码役偏好统计参量计算

运用CodonW1.4.2 和德泰生物的序列操作工具箱（http://www.detaibio.com/sms2/index.html），在线统计样本序列编码区长度（Sequence Length，缩写表示为SL），使用下述公式计算样本序列码役偏好统计参量。

1）GC 含量

式中：G、C 分别为样本基因现役码同位点碱基G、C 的出现次数；GC 为现役码同位点碱基G+C 含量；i 为碱基位点序号，取值1、2、3；N 为样本基因现役码总数，GC12为GC1和GC2的平均值，GCall为样本基因总体GC 含量。GC 含量反映方向性突变压力的强弱，GC3与码役偏好关系密切。

2）碱基偏倚度[9]

式中：PA、PG分别为样本基因现役码第3 位碱基A、G 的偏倚度，A3、T3、G3、C3分别为现役码第3 位碱基A、T、G、C 的含量，“|4”表示可计算4 种同义码的四、六码族的在编氨基酸，对于六码族的在编氨基酸，排除第3 位碱基相同的2 种现役码，仅取4 种进行计算。

根据奇偶规则（Parity Rule 2，PR2），DNA 互补链之间如果不存在选择性偏倚或突变，则碱基含量A 与T 相等，G 和C 相等。为避免密码子第3 位AT 和GC 之间突变不平衡，仅对四、六码族的在编氨基酸做现役码第3 位碱基分析[9]。

3）有效密码子数[1]（Effective Number of Codon，ENC）

中文文献通行译effctive 为“有效”，其实该处effctive 含义为密码子被基因序列役使征用而生效以实际编码氨基酸，即被基因役用，类似于服役。因此，译effctive 为“现役”不仅更为准确，而且能够避免因“有效”含义宽泛造成的误解。相应地，术语整体译为密码子现役数，简称现役码数。为交流方便，本文仍采用“有效密码子数”这一通行词。

式中：ENCbias为样本基因现役码偏好数，码役偏好（bias）受突变和选择压力双重影响；i 为氨基酸码族号（取值2、3、4、6）；j 为在编氨基酸序号；k 为现役码序号；n 为现役码数量；m 为在编氨基酸的现役码数目；N 为码族中在编氨基酸数目，小于或等于码族成员数；p 为现役码频数；f 为在编氨基酸的现役码一致性指数；F 为码族中在编氨基酸平均一致性指数，fij=0 时视同为非在编氨基酸，不参与F 值计算。ENCbias反映码役偏离随机选择的程度，最小值为20，反映每个氨基酸只有1个现役码的极端情况；最大值为61，反映全部氨基酸均在编且每个密码子均被均衡使用的极端情况。该值越小，码役偏好越强，基因表达强度一般越高。

式中：ENCmuta为样本基因现役码理论数，码役仅受到突变（mutation）影响；GC3为密码子第3位碱基G+C 含量。

式中：RENC 为样本基因现役码数比值，ENCmuta为现役码理论数，ENCbias为现役码偏好数。

4）同义码相对使用度（Relative Synonymous Codon Usage,RSCU）

式中：RSCU 为样本基因现役码相对使用度；x 为现役码数量；m 为在编氨基酸的现役码数目，取值2～6；i 为在编氨基酸序号；j 为现役码序号。RSCU＞1，表明其码役频率相对较高，反之亦然。

1.2.2 中性绘图分析

以GC12为纵坐标、GC3为横坐标，分析GC3与GC12的相关性，若二者呈显著相关关系，说明碱基的变异模式相同，码役受突变的影响较大；反之则表示基因有较高的保守性，码役受选择压力的影响较大[9]。

1.2.3 奇偶偏好绘图分析

以PA为纵坐标，PG为横坐标绘制散点图，分析碱基使用偏好及其影响因素。

1.2.4 有效密码子数绘图分析

以样本基因ENCbias、GC3为纵、横坐标，用R语言做散点图，并与ENCmuta曲线进行比较。当样本基因点分布于ENCmuta曲线附近时码役偏好受突变影响，样本基因点多分布于ENCmuta曲线下方较远位置时码役偏好较多受选择压力的影响[10]。

1.2.5 对应性分析

根据功能对每个样本基因进行分类，基于每个基因的RSCU 进行。将53 条样本基因分布到59 维RSCU 向量空间（排除一码族，共计59 个密码子，每个密码子为1 维），通过主成分分析降维后，检测分布在主向量轴上的样本基因的相对位置，推测码役的主要影响因素。

1.2.6 最优密码子分析

基于ENCbias排序53 个样本基因，从两极各选出10%各5 个基因作为高、低表达组分别计算RSCUtop、RSCUbot，计算RSCUtop、RSCUbot之差△RSCU。筛选△RSCU＞0.08 的现役码作为高表达优越密码子。与53 个样本基因总体现役码相对使用度RSCUall进行综合比较，将△RSCU＞0.08，且RSCUall＞1 的现役码定义为最优密码子[11]。

2 结果与分析

2.1 现役码组成分析

分析表1 得出，53 个样本基因总体GC1为45.93%，GC2为37.93%，GC3为29.10%。不同位置的GC 含量较大不同，呈现从第1 位到第3 位递减的分布趋势，样本基因现役码偏向于以A/T结尾。ENCbias大小可以反映码役偏好强弱[12]。53 个样本基因ENCbias范围为37.22～56.50，其中39 个基因ENCbias＞45，说明伯乐树叶绿体基因码役偏好较弱。

表2显示，GC3与GC1、GC2表现为相关不显著，而GC1与GC2极显著相关。结果表明，伯乐树叶绿体基因现役码第1 和第2 位碱基具有相似的组成，但与第3 位碱基不同。ENCbias值与GC 含量、SL 均无显著相关关系，说明伯乐树叶绿体基因碱基组成和序列编码区长度对码役偏好并非是主要影响因素。

表1 样本基因现役码GC 含量及有效密码子数

表2 基因现役码偏好统计参量的相关性分析

2.2 中性绘图分析

图1 表明，GC12与GC3的相关系数为0.016，双尾检验不显著，直线斜率为0.207 9，GC3与GC12的相关性很弱，突变对现役码第3 位碱基组成的影响与第1、2 位不同。说明伯乐树叶绿体基因码役偏好的形成受突变影响作用较弱，其它因素尤其是选择在此过程中可能起到重要的作用。

图1 中性绘图

2.3 奇偶偏好绘图分析

图2 显示，大部分样本基因分布于下方区域，说明在现役码第3 位T 的使用频率高于A，分布于左下方的基因略多于右下方，表明现役码第3 位碱基C 的使用高于G。由此推测，伯乐树叶绿体基因码役模式受突变和选择压力等因素的共同影响。

图2 奇偶偏好绘图

2.4 有效密码子数绘图分析

图3 显示，样本基因ENCbias值较分散，基因间码役偏好差别明显。表3 的RENC频率分析显示，大部分样本基因RENC分布在－0.05～0.10 之间，表明伯乐树叶绿体中该部分基因的码役偏好更多受突变的影响。rps14、petD、rps8、rps18 等基因远离ENCmuta曲线，位于曲线下方，表明伯乐树叶绿体中的这些基因受选择压力影响较大，码役偏好较强。

图3 有效密码子数绘图

表3 RENC频率分布

2.5 对应性分析

主成分分析计算结果显示，第一主分量轴体现9.45%的差异，第二、三、四主分量轴分别为8.32%、7.65%和6.77%。第一轴对样本基因码役偏好有较大影响。相关分析计算结果显示，第一主分量轴与ENCbias、GCall、GC3的相关系数分别为－0.208、－0.013、－0.064，且无显著相关性，因此GC 含量并非是第一主分量轴的单独影响因素。图4 显示，样本基因中遗传结构基因分布比较集中，表明伯乐树叶绿体遗传结构基因的码役模式相对一致。

2.6 最优密码子分析

图4 基于RSCU 的对应性分析

经过计算，用ΔRSCU ≥0.08 的标准筛选优越密码子，以RSCUall＞1 为高频密码子，以同属二者的现役码作为最优密码子，结果见表4。最终确定伯乐树叶绿体有16 个最优密码子，其中15个密码子以A 或U 结尾，仅1 个密码子以G 结尾。伯乐树叶绿体基因中ΔRSCU＞0.5 的密码子有7个，分别为GCU、GGU、UUG、CAA、CGU、UCU、ACU。

3 讨论

环境选择压力对密码子第3 位碱基的影响较小，第3 位碱基的变化通常不会改变氨基酸的对应关系，因此，经常把GC3作为分析密码子役佣的重要依据[13]。分析伯乐树叶绿体基因显示，现役密码子第1、2 位与第3 位碱基组成变异相关性不显著，现役密码子偏好与碱基组成、序列长度均无显著相关关系。奇偶偏好绘图分析表明，现役密码子第3 位碱基频率T＞A，C＞G，嘧啶频率高于嘌呤。该结果与陆地棉Gossypium hisutum[14]基本一致，与降香黄檀Dalbergia odorifera[15]、大花香水月季Rosaodorata var.gigantea[16]、马尾松Pinus massoniana[17]等不尽相同，这或许是物种间差异或统计方法不同所致。

表4 最优密码子有关统计参量

在有效密码子数绘图分析中，ENCmuta曲线代表了密码子偏好仅受第3 位GC 突变影响时基因的位置[14]。有效密码子数绘图分析结果显示，大部分样本基因分布于ENCmuta曲线附近，表明伯乐树叶绿体中大部分基因现役密码子偏好差异与GC3的差异有关，主要受到突变影响。结合碱基组成分析、对应性分析等结果，说明伯乐树叶绿体基因现役密码子第3 位GC 含量对密码子偏好存在一定影响，但并非主要影响因素。伯乐树叶绿体中仍有15 个基因，如rps14、petD、rps8、rps18等，分布于ENCexp曲线下方较远处，表明伯乐树叶绿体中这部分基因主要密码子偏好较强，受选择压力的影响较大。

伯乐树叶绿体基因最优密码子分析划分了3个等级，ΔRSCU ≥0.08 有6 个，ΔRSCU ≥0.3 的有3 个，ΔRSCU ≥0.5 的有7 个。共确定出16个最优密码子，其中15 个密码子以A 或T 结尾。这符合双子叶植物偏向于役用以A、T 结尾的密码子的结论[15]。

伯乐树叶绿体基因密码子偏好更多受到自然选择等因素的影响，而突变影响较弱。最终结合高频密码子和高表达优越密码子，确定GCA、GCU、UGU、GAU、GGU、AUU、AAA、UUA、UUG、CCA、CAA、CGU、UCU、ACU、GUA、GUU 等16个伯乐树叶绿体基因最优密码子。本研究为未来伯乐树叶绿体基因改造及密码子优化等提供了参考。

致谢：承蒙王俊杰先生悉心指导和详尽建议，本文修改后更加严谨，论述多有新突破。定稿之际，谨致谢忱！