叶绿体是绿色植物生命活动的代谢枢纽，为植物进行光合作用和生物合成的半自主性细胞器，具有独立的遗传系统——叶绿体基因组^[1]。关于叶绿体的起源有很多假说，目前广为接受的是内共生起源学说，即叶绿体起源于大约15亿年前被异养宿主细胞（真核原生生物）吞并且行光能自养的蓝细菌^[2]。由于早期内共生过程中蓝细菌经历了一系列复杂演化，发生基因丢失和基因转移事件，致使稳定下来的叶绿体基因组长度大大减小，但叶绿体基因组还是保留了少量原核特征，其中包括将基因组织成类似细菌操纵子的多顺反子转录单元，即基因簇^[3]。基因簇与细菌操纵子表面上看都是由功能相关的、共表达的、具有物理连接的基因群组成，但不同的是基因簇可能来源于物种演化过程中的基因复制、新功能化以及基因组的动态重组，而并非由细菌的水平基因转移所引起^[4]。陆地植物的叶绿体基因组通常具有结构稳定、基因顺序和基因含量高度保守的特征，其中一个主要限制因素就是它被组织成了多顺反子转录的基因簇^[5]。

种子植物的叶绿体基因簇很少被破坏，但有报道表明，在被子植物的一些植物类群如牻牛儿苗科^[6]、豆科^[7]和桔梗科^[8]中存在基因簇断裂的现象。在桔梗科植物中，两个叶绿体基因簇被破坏，clpP-5'-rps12-rpl20（clpP基因簇）和rps2-atpI-atpH-atpF-atpA（rps2基因簇）。rps2基因簇编码功能不相关的多肽，第一个基因（rps2基因）编码核糖体蛋白，其他4个基因编码ATP合酶亚单位^[9]。由于翻译系统和光合系统在蛋白质组分的需求数量上有很大的差异，这两种基因必须经过复杂的调控机制才能实现共同转录^[10]。该基因簇一旦发生断裂，重新分布的片段就需要新的启动子来驱动基因转录，然而如何在新的位置进行转录尚不清楚，基因簇发生断裂后对叶绿体进化的影响也知之甚少^[11]。

桔梗科为被子植物中高度进化的一个分支，约100属，2300 ～ 2380种；广布于全球，在药用和观赏方面有着重要的经济价值。桔梗科有广义和狭义之分，狭义桔梗科只包括具有传统意义上桔梗科基本特征的风铃草亚科。广义桔梗科包含5个亚科，即风铃草亚科、半边莲亚科、线枝草亚科、苣莲亚科和腔柱草亚科^[12]。尽管存在基因成簇分布的进化限制来让叶绿体基因组保持稳定，但在桔梗科的叶绿体基因组中还是发生了多次片段倒位事件，这种大规模的结构重排导致基因位置发生变化，也促使一些基因簇发生断裂。有研究发现，叶绿体基因组中的核苷酸替换速率在不同的谱系中、在反向重复区和单拷贝区之间以及在不同功能基因组之间均有所不同^{[13, 14]}。基因簇断裂可能会影响基因的转录进程，同时，成簇排列的基因又作为独立的编码蛋白基因，其进化速率可能会因位置的变动而发生改变。本研究对桔梗科典型的单种属（桔梗属）植物—桔梗（Platycodon grandiflorus (Jacq.) A. DC）的叶绿体基因组进行了测序，共选取50种广义桔梗科植物和2种睡菜科植物（附表1 ¹），提取其共有蛋白编码序列（Coding sequences，CDS）构建系统发育关系，对rps2基因簇在桔梗科中的分布情况进行统计；并在系统发育树的背景下，计算rps2基因簇中各基因的进化速率和选择压力，并检测其经历的正选择作用。

1.   材料与方法

1.1   取样及测序

桔梗新鲜叶片采自于河南省南阳县乔端镇自然生长的桔梗植株（33°34' N，112°6' E），采集后将叶片用锡纸包裹，液氮冷冻后于−80℃保存备用。

样品基因组 DNA的提取采用E.Z.N.A^® Plant DNA kit（OMEGA）试剂盒提取基因组，经片段化（300 ～ 500 bp）后进行纯化并构建测序文库。步骤依次为DNA末端修复、3' 端加A、连接测序接头、琼脂糖凝胶电泳法回收目的片段、对目的片段进行PCR扩增，最终建成测序文库。使用FastQC v 0.11.4软件对建好的文库进行质检，质检合格后采用 Illumina Novaseq^TM平台进行双端测序，测序读长为150 bp。接下来使用Trimmomatic v 0.32^[15]过滤数据得到clean reads，利用组装软件NOVOPlasty^[16] 对clean reads 进行拼接和组装。组装完成的基因组序列利用叶绿体基因组注释软件PGA^[17]进行基因注释，参考序列为毛细钟花（Codonopsis hongii Lammers）的叶绿体基因组全序列，并对所有样本的注释结果进行人工校正。通过在线网站OGDRAW v 1.3^[18] (https://chlorobox.mpimp-golm.mpg.de/OGDraw.html) 绘制叶绿体基因组图谱。最后，通过 Banklt 平台上传序列信息（https://www.ncbi.nlm.nih.gov/WebSub/）。

1.2   构建系统发育树

从NCBI（https://www.ncbi.nlm.nih.gov）数据库下载2种睡菜科植物和49种桔梗科植物的叶绿体基因组全序列，结合自行组装的桔梗序列，共计52种被子植物，涵盖桔梗科29个属和睡菜科2个属。

利用Genious prime 2019.0.4^[19]软件提取55条共有CDS序列，使用MEGA 7.0.26^[20]软件进行比对校正，并建立共有CDS数据串联数据集，以睡菜科的两个物种作为外类群构建邻接树（Neighbor-Joining，NJ），自展次数（Bootstrap）为1000次；最大似然（Maximum likelihood，ML）树使用RaxmlGUI2软件^[21]构建，设置GTRGAMMAI替换模型，1000自展次数；最大简约（Maximum parsimony，MP）树和贝叶斯（Bayesian inference，BI）树分别使用PAUP 4.0软件和MrBayes 3.2.6软件^[22] 进行构建，MrBayes 3.2.6软件中参数设置rates为invgamma，mcmc ngen为1 000 000。综合对比分析得到的4种树，结合前人发表的系统发育关系，挑选支持度高的系统发育树用于后续进化分析。

1.3   分子进化速率分析

使用HyPhy 2.2.4^[23]软件在系统发育背景下，基于最大似然法计算基因簇内5个基因（atpA、atpF、atpH、atpI和rps2）的转换率（Transition rate，trst）、颠换率（Transversion rate，trsv）、颠换率与转换率的比值（trsv/trst，ratio）、同义替换率（Synonymous substitution rate，dS）、非同义替换率（Nonsynonymous substitution rate，dN）和非同义替换率与同义替换率的比值（dN/dS，ω）。转换率、颠换率和ratio选择核苷酸分区和HKY85替换模型计算；同义替换率、非同义替换率以及ω选择密码子分区和MG94 × HKY85 × 3_4模型计算。利用SPSS 22软件的秩和检验进行显著性分析，使用GraphPad Prism 8软件绘图。

1.4   选择压力及适应性进化分析

运行PAML 4软件中codeml程序^[24]的分支模型（Branch model）、分支位点模型（Branch-site model）和位点模型（Site model），对rps2基因簇的相关基因进行进化分析。利用分支模型中最简单的单比率模型（One-ratio model，M0）、最复杂的自由比率模型（Free-ratio model，F）和位于中间的二比率模型（Two-ratio model，Model 2）进行似然比检验，选择更优模型。在二比率模型中将rps2基因簇发生断裂的物种设为前景支，未断裂的物种设为背景支运行分析。分支位点模型同样设置rps2基因簇发生断裂的物种作为前景支，检测正选择作用对部分位点的影响。采用可以降低假阳性的检验2筛选正选择位点，对MA模型和无效模型（ω值设置为1）进行差异显著性检验。位点模型中共包含4对嵌套模型：M1a（近中性，Nearly neutral）与M2a（正选择，Positive selection）模型、M0与M3（离散，Discrete）、M7（β）与M8（β和ω）模型、M8和M8a（β和ω = 1）模型。其中M1a与M2a、M7与M8模型用于正选择作用的检测，采用M1a与M2a检验结果为可能的正选择位点。鉴定正选择位点有NEB（Native empirical）和BEB（Bayes empirical bayes）两种贝叶斯方法，前者易产生假阳性，本研究主要以BEB的结果为准。

将筛选出的可能存在正选择位点的氨基酸序列提交到瑞士生物信息研究所（Swiss-Model, https://swissmodel.expasy.org/），利用同源建模原理预测编码蛋白的三维结构，查看正选择位点所在的结构位置。

2.   结果与分析

2.1   桔梗的叶绿体基因组信息

桔梗的叶绿体基因组为典型的4区环状结构，基因组大小为173 118 bp，大单拷贝（Large single copy，LSC）区为79 024 bp，反向重复（Inverted repeats，IRs）区为43 127 bp，小单拷贝（Small single copy，SSC）区为7840 bp；基因组总GC含量为38.1%，其中LSC区GC含量为37.2%，IR区为39.5%，SSC区为31.0%；基因组全序列总计110个编码基因，包含77个蛋白质编码基因、29个tRNA基因和4个rRNA基因（图1），提交至NCBI的序列登录号为OP018680。

图  1  桔梗的叶绿体基因组图谱

圈内基因按照逆时针方向转录，圈外基因按照顺时针方向转录。

Figure  1.  Chloroplast genome of Platycodon grandiflorus

Genes inside circle are transcribed counter-clockwise, and genes outside circle are transcribed clockwise.

下载: 全尺寸图片 幻灯片

2.2   系统发育关系

用4种方法构建出来的系统发育树中（附图1 ^2）），贝叶斯法和最大似然法的结果一致，邻接树和最大简约树与前面两种树的差别在于半边莲亚科的部分分支，这些分支的rps2基因簇并未发生断裂。贝叶斯树除了两个分支的自展值分别为0.94和0.99外，其余分支自展值均为1，支持度最高。综合前人研究结果，最终选取基于贝叶斯树的系统发育关系进行后续分析（图2）。

图  2  所选物种构建的贝叶斯树

Figure  2.  Bayesian inference tree of sampled species

下载: 全尺寸图片 幻灯片

2.3   rps2基因簇在桔梗科中的分布信息

rps2基因簇（rps2-atpI-atpH-atpF-atpA）在所有物种中均位于LSC区。在桔梗科的50个物种中，39个物种的rps2基因簇未发生断裂，11个物种发生了断裂，断裂位置均位于atpI和atpH基因之间，且这11个物种均来自风铃草亚科，在系统发育树上聚为一支（图2）。其中，8个物种断裂后两个片段的排列方式一致但相距较远，在系统发育树上聚为一支，剩余3个物种的断裂片段分布位置各不相同，转录方向或者排列顺序发生改变，两个片段之间的距离缩短，这3个物种在系统发育树上也聚为一支（图3）。

图  3  rps2基因簇在桔梗科中的断裂片段分布

Figure  3.  Distribution of broken rps2 gene cluster in Campanulaceae

下载: 全尺寸图片 幻灯片

2.4   rps2基因簇的进化速率分析

对基因簇中5个基因的进化速率参数进行分析，根据基因簇是否断裂，将50个物种分为未断裂（39种）和断裂（11种）两组。在均值上，除了atpH基因的ω，其余未断裂组中相关基因的进化速率参数均高于断裂组，为断裂组的1.13 ~ 7.86倍。曼惠特尼U检验（Mann-Whitney U test）结果显示，具有显著差异的有： atpA基因的dS（P = 0.023）、atpH基因的trst（P = 0.040）、atpI基因的ratio和dS（P = 0.038、P = 0.026）以及rps2基因的trst、ratio、dS（P = 0.011、P = 0.047、P = 0.010）（图4）。

图  4  rps2基因簇未断裂与断裂两组的进化速率均值 ± 标准差比较

1表示rps2基因簇未断裂的物种， 2表示rps2基因簇断裂的物种。* 表示两者的进化速率参数在P < 0.05水平上差异显著。

Figure  4.  Comparison of mean ± SD of evolutionary rate parameters of rps2 gene cluster in unbroken and broken groups

1 is species with unbroken rps2 gene cluster; 2 is species with broken rps2 gene cluster. * indicates significant differences at P < 0.05 level in evolutionary rate parameters.

下载: 全尺寸图片 幻灯片

rps2基因簇发生断裂的11个物种根据断裂片段排布方式可分为4组，桔梗（Platycodon grandiflorus (Jacq.) A. DC）、Codonopsis minima Nakai和毛细钟花（Codonopsis hongii Lammers）的进化速率参数大多高于断裂方式一致的8个物种均值（附图2 ^3））。秩和检验结果显示，atpI的trst、ratio、ω和dS（P = 0.019、 P = 0.040、 P = 0.019、 P = 0.040），rps2的ω和dS差异显著（P = 0.019、 P = 0.019）。将桔梗、C. minima和毛细钟花3个物种与未断裂组进行均值比较，发现这3个物种的很多值也高于未断裂组，秩和检验结果显示所有参数均无显著差异（附图3 ^3））。对5个进化速率参数做秩和检验，结果显示在桔梗科物种（50个）中，除了dS无显著差异外（P = 0.678），其余参数值均存在显著差异（P_trsv = 0.000、 P_trst = 0.028、 P_ratio = 0.000、 P_ω = 0.000、 P_dN = 0.000）（附表1 ^3））。

2.5   选择压力及适应性进化分析

分支模型中，以断裂组为前景支，未断裂组为背景支，M0与Model 2的似然比检验结果显示，rps2基因在前景支和背景支之间受到的选择压力存在显著差异（ω_前景支 = 0.875，ω_背景支 = 0.431， P = 0.029）；atpA、atpF、atpH和atpI 基因无显著差异，即每个分支受到的选择压力一致（ω_atpA = 0.089、ω_atpF = 0.356、ω_atpH = 0.026、ω_atpI = 0.121）。分支位点模型中前景支与背景支的选择压力均无显著差异，未检测到正选择位点。位点模型中，atpA、atpF、atpH基因的M1a和M2a模型的似然比检验结果无显著差异（P = 1.000、P = 1.000、P = 1.000），atpI、rps2基因的似然比检验结果存在显著差异（P = 0.021、P = 0.000）。atpI基因的M2a模型中检测到1个正选择位点124H（后验概率 > 95%）；rps2基因中M2a模型中检测到7个正选择位点23D、70R、105P、106N、127N、164K、240R，其中23D、70R、127N 3个位点的后验概率超过99%（附表2、3 ^4））。

基于M1a与M2a的检验结果，将ATPI和RPS2的蛋白质序列导入到Swiss-Model中，以桔梗为例，基于同源建模原理得到蛋白质的三维结构图，atpI基因的M2a模型检测到的正选择位点124H位于S（弯曲）位置。rps2基因的M2a模型检测到的正选择位点23D和164K位于E（β-折叠），70R位于C（环或无规则卷曲）位置，105P和106N位于T（氢键转角）位置，127N位于H（α-螺旋）位置，240R位于不规则区域（图5）。

图  5  ATPI（A）和RPS2（B）蛋白质的三维结构

Figure  5.  Three-dimensional structure of ATPI (A) and RPS2 (B) protein

下载: 全尺寸图片 幻灯片

3.   讨论

3.1   桔梗科rps2基因簇的断裂及其系统发育意义

桔梗科植物叶绿体基因组中有大规模的重组变异现象发生，包括倒位、基因丢失以及IR区的收缩/扩张等^[25]。最早在疗喉草（Trachelium caeruleum L.）的叶绿体基因组中报道过大规模的重排现象，包括7 ～ 10个大倒位片段、2个基因丢失、IR区的收缩/扩张以及基因簇的断裂等^[26]。重排片段的两端往往伴随着大量重复序列和基因复制，而断裂的基因簇片段与重复序列一样，也分布在重排片段的末端^[27]。保守叶绿体基因簇在序列、基因含量和排列上的差异可能有助于系统发育分析^[28]。在桔梗科中，rps2基因簇发生断裂的位置皆位于atpI ~ atpH基因之间，且发生断裂的物种均来自风铃草亚科，在系统发育树上聚为一支；断裂片段摒弃原本紧密连接的物理位置进行重新排布，重新排布方式相同的8个物种在系统发育树上也聚为一支；断裂片段排布方式各不相同的3个物种均来自蓝钟花族，在系统发育树上同样聚为一支，说明rps2基因簇断裂可能具有系统发育意义。本研究结果验证了前人的猜想，同时也为构建桔梗科复杂的系统发育关系提供了一个新思路。另外，rps2基因簇发生断裂的物种均位于桔梗科系统发育关系的基部分支，随着物种的不断进化，断裂现象越来越少，表明rps2基因簇可能更倾向于成簇转录。

3.2   桔梗科rps2基因簇与进化速率的关系

研究发现，在细菌和动物线粒体基因组以及被子植物的叶绿体基因组中，核苷酸进化速率和基因组重排率之间存在显著的正相关^[29-31]。究其原因，有人认为在线粒体中这种相关性与DNA复制的准确性有关，在叶绿体基因组中则可能是由DNA修复机制的异常所引起^{[6, 11, 31]}。被子植物叶绿体的DNA修复中涉及4类核编码基因，它们产生抑制重复DNA序列间重组的蛋白质，通过防止不正当重组保证基因组的稳定性，一旦编码基因发生突变就会影响错配修复系统的正常运作，使得叶绿体基因组极其不稳定^[32]。那么，这种相关性很有可能是由错配修复系统受到干扰的结果。与其他物种相比，rps2基因簇发生断裂的物种经历了更多的重排事件，但其相关基因进化速率的均值降低；其中，非同义替换速率无显著差异，同义替换速率在atpA、atpI、rps2基因中显著降低。根据断裂片段的排布方式进一步探究，发现桔梗、C. minima和毛细钟花3个物种中多个基因的进化速率值均高于未断裂组，而断裂后atpA-atpF-atpH片段不发生变动，仅atpI-rps2片段进行重新排布的8个物种，其进化速率参数大多为零或接近零，这种结果可能是计算中的假阳性导致，而这也是基因簇断裂后进化速率均值降低的主要原因。同时，rps2基因簇断裂后的进化速率表现出的不一致性变化，可能与断裂片段的重新排布方式或断裂片段间的距离有关。

3.3   桔梗科rps2基因簇相关基因的适应性进化

生物进化的本质是基因进化^[33]。对叶绿体结构和功能特征的研究表明，质体基因可能不是作为单个位点进化，而是经历着不同的进化力量^[13]。由于基因功能的差异以及各个基因在染色体上的位置不同，基因之间关系复杂，不同基因的核苷酸替换速率有时相差很大，非同义替换速率在基因间的差异比同义替换速率的差异大^[33]。rps2基因簇这5个基因的非同义替换率存在显著差异，而同义替换速率并未出现差异，可能是因为rps2基因簇的5个基因分别隶属于两类功能基因。这也体现出了基因间的进化速率异质性，表明进化速率在不同的蛋白编码基因中可能是独立且高度可变的。

叶绿体ATP合酶（CF₁-CF_o ATP synthase）主要利用光合作用过程中产生的跨膜质子动力势催化ATP的合成，使得光能转化为化学能^[34]。它包含2个亚复合物—CF₁和CF_o，CF₁主要参与ATP的合成或水解，由α、β、γ、δ和ε 5种亚基构成；CF_o构成跨膜质子通道的同时也是固定CF₁的支架，由Ⅰ ~ Ⅳ 4种亚基构成^[35]。在rps2基因簇的这5个基因中，atpA、atpF、atpH和atpI分别编码CF₁的α亚基和CF_o的Ⅰ、Ⅲ和Ⅳ亚基^[36]。其中atpI基因中检测到了正选择位点，说明该基因经历了适应性进化，可能会对植物的光合作用及能量转换产生一定的影响。rps2基因是小亚基核糖体蛋白的一种，编码核糖体蛋白S2，参与核糖体前体的加工和组装，主要与蛋白S3、S5、S8交联，同16S RNA的3′结构域相互作用^[37]。该基因也是叶绿体基因组中关键的功能基因，在物种适应不同环境中起着十分重要的作用。

4.   结论

桔梗科叶绿体基因组大规模的重排导致rps2基因簇发生断裂，这种断裂可能具有一定的系统发育意义。rps2基因簇断裂后，进化速率均值降低，但在不同的断裂方式中具有差异性。同时，基因间进化速率和选择应力的差异可能与基因的功能和结构相关。本研究揭示了叶绿体基因之间广泛的进化速率异质性和基因组结构变异，为桔梗科叶绿体进化提供了新的见解，但由于本研究选取物种有限，得出的结论可能不具备普遍性，今后还应扩大物种数量，以获得更准确的结论。