基于转录组数据分析药用大黄的密码子使用偏好性

高梦琦，邹建珍，霍小位，李依民，张明英，高 静，张 岗*，刘蒙蒙*

1.河北大学中医学院，河北 保定 071002

2.河北大学药学院，河北 保定 071002

3.陕西中医药大学药学院 陕西省秦岭中草药应用开发工程技术研究中心，陕西 西安 712046

药用大黄Rheum officinaleBaill.属于蓼科大黄属多年生草本植物，为《中国药典》2020年版收录的大黄3种基原植物之一，具有泻下攻积、清热泻火、凉血解毒、逐瘀通经、利湿退黄等功效[1]。现代中药化学研究表明大黄中主要含有蒽醌类、酚类等物质，其中蒽醌类化合物是主要的活性物质，具有抗炎、抗肿瘤作用[2-3]。目前药用大黄的资源主要来源于人工种植，随着过度采挖，野生资源已经濒临灭绝。因此，开展蒽醌类化合物生物合成机制研究，探索蒽醌类化合物生物合成路径，是解决药用大黄野生资源濒危的思路之一。

遗传密码子是生物体DNA与蛋白之间信息传递的信息载体。大多数氨基酸不仅由一种密码子编码，而是可以同时由2～3种密码子编码，生物对于不同的密码子的使用具有一定偏好性[4]。研究发现，密码子使用偏好性受到GC含量（鸟嘌呤和胞嘧啶含量）[5]、基因表达水平[6]、环境压力[7]、突变压力[8]、基因长度[9]、蛋白结构[10]及tRNA丰度[11]的影响。随着高通量测序技术的发展，药用植物的密码子使用偏好性研究有了显著的成果[12-13]，但是目前还没有关于药用大黄密码子偏好性的研究。生物体选择最优密码子可以体高基因的翻译效率及蛋白的表达量。因此，研究密码子的使用偏好性，有利于基因工程研究中密码子优化、最适外源宿主选择，并对生物的系统进化研究具有重要意义。

本研究利用药用大黄的转录组数据对其编码序列进行密码子偏好性分析，揭示药用大黄基因密码子使用偏好性的特征，为药用大黄基因基因工程研究中进行密码子优化及选择最佳的外源载体提供理论基础。

1 数据

药用大黄转录组（SRR10838870）来源于本课题组前期选取药用大黄种苗根、根茎、样品，利用Illumina HiSeq 2500 测序平台进行转录组测序并拼接获得Uingenes序列，利用BLASTx及ESTScan软件对上述Uingenes序列进行编码序列（coding sequence，CDS）分析。使用perl程序筛选编码序列中长度大于300 bp并且同时具有起始密码子及终止密码的序列，用于药用大黄密码子使用偏好性分析，最后共得到4 733条CDS进行分析。

2 方法

2.1 编码序列中碱基含量及中性绘图分析

使用Codon W软件分析药用大黄编码序列中的有效密码子数（effective number of codons，ENc）、密码子适应指数（codon adaptation index，CAI）及同义密码子相对使用度（relative synonymous codon usage，RSCU），计算序列中的GC含量，利用编写的perl程序计算GC1、GC2及GC3的含量，即密码子第1、2、3位的GC含量。GC12为GC1、GC2的平均值。

中性绘图分析主要用于分析GC12与GC3之间的相关性，以GC12为纵坐标，以GC3为横坐标制作散点图，并进行线性拟合回归分析。可用于评价自然选择压力及突变对密码子使用偏好性造成的影响。

2.2 GC3关联及奇偶偏好性分析

ENc-GC3关联分析主要用于分析密码子使用偏好性与碱基组成之间的关联。以ENc为纵坐标、GC3s为横坐标绘制二维散点图，并进行非线性拟合回归分析。如果基因的点分布在拟合曲线附近，则表明密码子的使用偏好性主要是受到突变作用的影响；如果基因点分布在拟合曲线下方，则表明密码子的使用偏好性主要是受到选择压力的影响。奇偶偏好性分析以A3/(A3＋T3)为纵坐标，以G3/(G3＋C3)为横坐标制作二维散点图，以（0.5，0.5）坐标为中心点，从这个中心点制作的失量表示了奇偶偏倚的程度和方向[14]。

2.3 最优密码子分析

对于最优密码子的筛选，本研究利用药用大黄编码书写的CAI值进行由高到低排列，然后取前后各5%的编码序列，分别代表高表达与低表达的基因序列。之后分析各个密码子的平均RSCU值。当密码子在高低表达基因中RSCU值的差异＞0.3，且在高表达基因中RSCU≥1、在低表达基因中的RSCU≤1，即可被确定为最优密码子[15-16]。

2.4 对应分析

对应分析是指利用多元统计分析的方法研究不同编码序列密码子的变异趋势。利用Codon W软件中对应分析功能，将所有编码序列分布到59维向量空间，进而获得密码子使用偏好性的特点。通过分析基因位置研究基因向量及基因间值的变异程度，预测使密码子使用发生偏好性的原因[17]。

3 结果

3.1 密码子GC含量组成分析

利用Codon W软件计算药用大黄编码序列中GC含量，将药用大黄编码序列作为一个整体进行计算。结果如图1所示，GC含量范围分布在28.6%～71.9%，主要分布在40%～50%，GC的平均含量为45.6%。GC1、GC2、GC3的平均含量分别为48.63%、40.89%、44.37%。对比发现，GC2的含量最小，GC1与GC3的值相对比较接近，说明药用大黄编码序列更倾向于使用鸟嘌呤（G）和胞嘧啶（C）这2种碱基。

图1 药用大黄编码序列中GC含量组成Fig.1 Distribution of GC content ofcoding sequences(CDS) in R.officinale

3.2 中性绘图分析

药用大黄编码序列的中性图分析如图2所示，其中GC12的数值范围为0.15～0.87，GC3的数值范围为0.09～0.87。GC12与GC3的相关性系数为0.155，表明这两者之间的相关性比较弱，说明选择压力在药用大黄密码子使用偏好性的选择具有一定影响。

图2 中性图分析Fig.2 Neutrality plot

3.3 ENc-plot绘图分析

ENc的取值范围为20～61，平均值为47.58，其中ENc≤35的基因有386个（8.15%），ENC＝61的基因有251个（5.30%），说明药用大黄的编码序列除了少数密码子具有较强的偏好性，其整体的密码子使用偏好性较弱。以ENc为纵坐标、GC3s为横坐标绘制散点图（图3），图中各点代表基因分布情况。图3显示大部分基因远离期望曲线，说明密码子的偏好性主要来源于选择压力的影响，少部分基因均匀的分布在期望曲线两侧，说明部分基因收到突变压力的影响。

图3 药用大黄转录组编码序列ENC-GC3s相关性分析Fig.3 Neutrality plot (ENC vs GC3)

利用(ENCexp－ENCobs)/ENCexp统计分析期望ENc与观测ENc之间的区别[18]。如图4所示，(ENCexp－ENCobs)/ENCexp的取值主要分布在0～0.2，说明期望ENc与观测ENc两者之间的存在很小差异。进而说明在药用大黄编码序列密码子使用偏好性形成过程中，突变压力不是唯一的影响因素，但却是一个重要的影响因素。

图4 ENc比值频率分布Fig.4 Frequency distribution of ENc ratio

3.4 对应性分析

对应性分析结果见图5，GC含量不同的基因被标记上不同的颜色，编码序列中GC含量大于60%的被标记为红色，GC含量在45%～60%的被标记为蓝色，GC含量小于45%的被标记为绿色。其中，GC含量大于60%的编码序列分散的分布在坐标系中，而GC含量小于60%的编码序列则在坐标系中分布比较集中。

图5 对应性分析Fig.5 Correspondence analysis

此外，选取6个编码序列相关的参数（GC12、GC3、GCall、ENC、CAI以及Axis 1）进行相关性分析，结果如表1所示，Axis 1与GC12、GC3这2个重要参数存在显著的负相关性（r＝-0.418，r＝-0.728，P＜0.01）。结合图5与表1，可以说明GC含量与GC3含量高的编码序列取向于分布在Axis 1左侧。因此，推测以G/C结尾的密码子应集中分布在Axis 1负轴而以A/U结尾的密码子集中分布在Axis 1正轴。

表1 药用大黄编码序列对应参数的相关性分析Table 1 Correlation analysis of R.officinale gene-related parameters

3.5 偏倚性分析

偏倚性分析结果如图6所示，大部分药用大黄转录组编码序列样本点分布于中心点附近，说明其碱基中嘌呤与嘧啶出现的频率相对一致，少部分编码序列偏离中心点较远，进一步说明药用大黄密码子使用偏好性形成过程中主要受到突变压力的影响。

图6 编码序列偏倚性分析Fig.6 Parity rule 2 (PR2 ) -bias plot

3.6 基因表达水平对密码子使用偏好性的影响

药用大黄编码序列的CAI的范围为0.261～0.835，表明药用大黄的基因表达水平比较高。同时，对CAI与其他几个重要的参数进行了相关性分析，包括ENc、GC3、GC。结果如表1及图7所示，CAI与ENc及GC之间都存在显著的正相关性（r＝0.07，P＜0.01；r＝0.695，P＜0.01）。因此，基因表达水平在药用大黄密码子使用偏好性形成过程中起到关键作用，基因表达水平高并且基因中GC含量高，密码子使用偏好性程度越高。

图7 中性图分析Fig.7 Neutrality plot

3.7 最优密码子分析

药用大黄编码序列高、低表达序列的RSCU值如表2所示。从中筛选到29个密码子作为药用大黄编码序列的最优密码子，在这些最优密码子中，以U（15个）及A（13个）结尾的较多，表明药用大黄转录组编码序列的最优密码子偏好以U或A结尾。

表2 药用大黄转录组高、低表达基因的RSCU比较分析Table 2 Optimal codons of R.officinale genes based on the RSCU values

4 讨论

在生物长期的进化过程中，形成了对密码子的使用的固定模式。今年来，随着二代测序技术的发展，探索天然产物合成路径并进行异源合成，成为了研究热点。因此，探索药用大黄基因密码子使用偏好性，可为研究其密码子使用特点，筛选最佳蛋白异源表达载体及密码子优化奠定理论基础。

本研究对药用大黄转录组数据中4733条完整编码序列的密码子使用偏好性进行分析，结果显示，药用大黄转录组编码序列的GC、GC3平均含量分别为45.6%、44.73%，表明药用大黄转录组中的编码序列偏好使用富含A、U且以A/U结尾的密码子，这与小麦Triticum aestivumL.[19]、水稻Oryza sativaL.[20]、玉米Zea maysL.[21]密码子使用模式不一致，这些植物的密码子富含GC碱基，并且大多数密码子以G/C结尾，药用植物川贝母转录组中的编码序列的GC和GC3平均含量都接近于50%，并没有表现出密码子使用偏好性[15]，表明碱基组成在药用大黄响密码子偏好性形成过程中起到重要作用。

此外，碱基突变和遗传选择也影响着密码子使用偏好性形成的重要因素。本研究中，药用大黄转录组GC12与GC3含量具有显著相关性，ENc-plot绘图及偏倚性分析分析结果表明，大部分基因远离于标准曲线及坐标轴中心点，基因的期望ENc值与观测ENc值具有较大差异，说明这部分基因密码子偏好性的形成是以遗传选择为主导因素；少数的基因分布于标准曲线附近及坐标轴中心点，基因的期望ENc值与观测ENc值基本保持一致，说明这部分基因密码子形成偏好性的形成是以突变为主导因素。

在构建蛋白表达载体时，不同的基因密码子使用偏好性差异很大，因此分析药用大黄的密码子使用特点对于载体的密码子优化具有重要意义。本研究筛选得到药用大黄29个最优密码子，大多数密码子以A/U结尾，从药用大黄编码序列中的GC含量可以看出, 其偏好使用A或U碱基结尾的密码子。由此推断药用大黄基因的密码子偏好性的产生与其基因进化过程中AU到GC突变压力要比GC到AU突变压力高有关。

本研究对药用大黄转录组数据中编码序列的密码子使用特点进行了研究，探讨了影响其密码子使用偏好性形成的影响因素，并筛选出药用大黄基因的29个最优密码子，可为蒽醌类化合物的异源合成及其合成路径研究提供理论依据。

利益冲突所有作者均声明不存在利益冲突