在异源多倍体起源的谱系中，基因内容和同线结构不同的同源染色体组可能共存。染色体块上基因和微同线性的存在或不存在可以用于区分亚基因组并推断系统发育。我们在这里应用基因组结构数据来推断古老异源多倍体谱系——核桃科（胡桃科）的关系，通过使用七个染色体水平基因组，其中两个是新组装的。微同线性和基因内容分析产生了相同的拓扑结构，这些拓扑结构与20世纪80年代形态学分支研究中化香和黄杞的拓扑结构相同。这里和许多早期研究中基于DNA比对的拓扑学将化香与山核桃和核桃归为一类，这可能是受到过去杂交的误导。所有可用的数据都支持胡桃科的杂交起源，其起源于已灭绝或未取样的祖先，嵌套在杨梅科或其姐妹中。Rhoiptelea chiliantha是所有其他胡桃科的姐妹，包含更多的DNA修复基因，并且似乎以2.6至3.5的速率进化。比其余物种慢几倍。核桃科，胡桃科(Juglandaceae)包括8个属63个物种（图1），包括一些世界上最具商业价值的坚果生产作物，如波斯核桃、中国铁核桃（均为核桃属）、山核桃和山核桃（山核桃属））。许多树种也是有价值的用材树种。由于其经济重要性，已经组装了几个基因组来研究真菌抗性基因和与果实品质相关的基因，他们的分析揭示了家族[1-5]中存在古老的全基因组重复(WGD)，称为胡桃类WGD[1]。此外，对核桃和伊利诺伊山核桃基因组的比较表明，在这些物种中持续存在的两个亲本亚基因组之间存在重复基因的偏向分离和不对称丢失[3,5,6]，表明该家族的无源四倍体起源。细胞遗传学研究也提出了这一点，通过假设与杨梅科共享或嵌入杨梅科的祖先来解释2n=32或64的主要核桃科染色体数，杨梅科的基本数量为x=8[7.8]。尽管过去20年进行了大量的分子研究，但该科普遍认可的属之间的关系仍不清楚，至少有6种不一致的拓扑被报道[9-19]，尽管许多这些研究中的主要信号来自相同的叶绿体matK基因序列。最成问题的是神秘的东亚特有化香树的地位，它在主要分支之间“跳跃”。今天，该物种出现在越南、台湾岛、中国大陆、韩国和日本，但化香叶子和果实化石是在北美和欧洲的古新世晚期和始新世早期发现的[20,21]。化香在核桃科中因其圆锥状果实而不同寻常（Manning[22]；我们的图1），但美洲对虾和栲树的化石叶子与现代黄杞科[20]具有独特的叶子结构特征，而现存的化香和黄杞还具有薄坚果没有空隙的壁[22]、雄花和雌花通常结合在一起的顶生圆锥花序（Stone[23]：图7.1）以及类似的花粉（Stone[23]：图7.3）。至少四个现代形态数据矩阵[9,14,19,24]无法明确定位化香属，即使包括化石类群（讨论）。然而，第一个此类分析[20]推断出化香/黄杞分支（另见曼彻斯特[25]；我们的图1a）。其他基因组和群体遗传学研究强调了核桃科最多样化的属Juglans[26,27]中古代和正在进行的杂交的作用。如果跨物种基因流（基因渗入）影响了大部分相关基因组，则基于DNA序列比对的系统发育推断将生成反映网状历史而不是分叉历史的物种树[28,29]。在这种广泛渗入的情况下，微同线性似乎保留了系统发育信号[30]。分子钟假说在分子进化、物种形成事件计时和历史人口学研究中发挥着核心作用。然而，不同谱系的替代率可能有很大差异，并且基因重复可以显示出不同的时间和功能进化模式[31]，与世代时间[32]、其他生活史特征（例如，Lanfear等人[33]）、不同的种群规模[34-38]或可能与DNA有关修复系统[34,39,40]。由于比较不同类群中参与DNA修复的基因数量存在技术难度，最后一种可能性迄今为止受到的关注最少。在这里，我们利用同线性信息和基因存在/缺失来重建核桃科的系统发育历史，利用其两个共存的同源亚基因组和染色体水平基因组组装。我们的基于基因内容和微同线性的系统发育方法建立在Zhao等人[30]的工作基础上。和Pettetal.[41]，但与他们的研究不同，他们的研究重点是异源四倍体进化枝。解决核桃家族的系统发育对于通过了解该家族性状进化的方向来更好地利用可用的基因组非常重要，因为正确的根拓扑对于生物地理学和分子钟研究至关重要。我们还推断了DNA替代率、过去的种群规模和DNA修复基因的相对数量，以测试替代率可能随着DNA修复基因家族的扩展而变化的可能性。在其他近期研究的背景下，我们的结果强调了核心真双子叶植物木本分支中多倍体后进化急剧放缓的模式，呼吁对此类分支的分子钟测年更加谨慎。Rhoipteleachiliantha的总组装基因组包含408.19Mbp，重叠群N50和支架N50值分别为6.97Mbp和24.34Mbp。在整个基因组中，91.71%被排序并定向为16条伪染色体（补充图1）。R.chiliantha基因组由约39.35%的重复序列和32,505个预测的蛋白质编码基因组成（补充表1）。Engelhardiaroxburghiana组装大小为884.78Mbp，整个基因组的97.73%有序并定向为16条假染色体（补充图1）。E.roxburghiana基因组由约57.11%的重复序列和30,590个带注释的蛋白质编码基因组成（补充表1）。组装的完整性评估显示，BUSCO胚胎植物_odb9数据库中96.1％和93.7％的通用单拷贝直向同源物分别存在于R.chiliantha和E.roxburghiana中，表明这两个基因组组装相当完整（补充表2）。R.chiliantha和E.roxburghiana的伪染色体的圆形图和点图在同源染色体之间产生了丰富的同线性信息（图2a和补充图2、3）。共线块（同线块）的数量为399个（山核桃）、403个（山核桃）、324个（核桃）、444个（山核桃）、338个（小果山核桃）、384个（化香）和463（R.chiliantha），分别。为了确定每个物种中最近和更古老的Ks峰(WGD)，我们绘制了物种内每个同线性块中包含的共线基因对的中位Ks值的分布。在Ks分布中，我们在C.illinoinensis、E.roxburghiana、J.mandshurica、J.microcarpa、J.regia、P.strobilacea和R.chiliantha的每个基因组中检测到两个多倍化事件（图2b）。R.chiliantha的Ks第一个峰值为0.17，而其他6个物种的Ks为0.36-0.48。所有物种的第二个Ks峰都很窄，指向共享的古老γ-WGT。四重同义第三密码子颠换率（4DTv）的分布进一步支持了这些结果（补充图4）。图1 核桃科基于形态学和分子的系统发育。系统发育树是根据Wing和Hickey对30个叶子、果实、花序和花粉特征进行的分支分析修改而来的，这些特征为现存和灭绝的类群（这里省略了后者），其中Rhoiptelea、Myricaceae和Fagaceae为外类群。有关最近四项形态分支分析的结果，请参阅正文和补充图17。Cyclocaryapaliurus(Batal.)Iljinsk。在此被视为继22,98和Alfaropsisroxburghiana(Wall.)Iljinsk之后的Pterocaryapaliurus。被视为Engelhardiaroxburghiana（比较补充图16）。水果照片并未按比例绘制，而米亚比例尺是基于各属最古老的化石出现，但没有进行正式分析。b右侧的系统发育树改编自Mu等人15，他们使用核RAD-Seq和全叶绿体比对数据。图2 全基因组重复检测。a一个Rhoipteleachiliantha同源染色体的圆形图。MCScanX中使用了相互最佳命中同源基因对。中心线连接染色体上的共线块，其上显示基因密度。不同的颜色代表源自最近共享WGD事件的不同同源染色体对。b为七个核桃科物种的每个共线块的中值Ks的Ks分布。源数据作为源数据文件提供。为了确定R.chiliantha是否与其他核桃科物种共享最近的WGD事件，我们使用了基因组同线性分析[42,43]和Pfeil等人的“多基因家族”系统发育树方法[44]。使用J.regia作为胡桃科sensustricto的代表。在基因组同线性分析中(图3a)，R.chiliantha和J.regia的染色体形成八个四联体，每个四联体由一对来自R.chiliantha和J.regia的同源染色体组成。点的颜色清楚地区分了种间同源染色体之间的直系同源(黄色：较小的Ks)和旁系同源(绿色：较大的Ks)关系，表明它们之间具有共享的WGD。Rhoipteleachiliantha和胡桃科的其余物种显示出类似的模式(补充图5-9)。利用Pfeil等人的系统发育树方法，超过92%的基因树(bootstrap值≥80)支持拓扑结构(((R.chiliantha,J.regia)(R.chiliantha,J.regia)),Quercuslobata),(((R.chiliantha,J.regia),R.chiliantha),Q.lobata)and(((R.chiliantha,J.regia),J.regia),Q.lobata)(图.3b和补充数据1)；核桃科的其余物种也表现出类似的模式（补充数据1）。这两个结果提供了明确的证据，表明R.chiliantha和其余的核桃科共享相同的WGD事件，该事件一定存在于它们最近的共同祖先中。尽管R.chiliantha和其他胡桃科共有胡桃属WGD，但它们的同源基因对的Ks分布显着不同。为了量化R.chiliantha和其他核桃科物种的替代率差异，我们计算了焦点物种与其姐妹分支的共同祖先之间的Ks值（更多详细信息参见图3c及其图例）。通过这种方式计算，Ks现在反映了一个物种自其共同祖先分化以来的进化速率（图3c）。对R.chiliantha和其他物种的Ks分布进行了Wilcoxon秩和检验（Mann-Whitney检验）。结果表明，R.chiliantha与其他6个物种之间的Ks分布存在显着差异（P<0.01；图3d和补充图10）。接下来我们使用等式量化了速率。1，根据BudvaricarpusSerialis假设为85MYA，BudvaricarpusSerialis是一种与R.chiliantha[45]相似的果实结构化石，也是核桃科最古老的化石。它来自晚土伦阶到桑顿阶，这个时期跨越89.8到83.6MYA。每个站点的替代率估计为1.60×10−9,1.95×10−9,1.47×10−9,1.54×10−9,1.51×10−9,1.89×10−9,0.55×10−9分别为C.illinoinensis、E.roxburghiana、J.mandshurica、J.microcarpa、J.regia、P.strobilacea和R.chiliantha（补充表3）。与其他核桃科植物相比，Rhoipteleachiliantha的有效种群规模较小根据近中性理论，有效种群规模（Ne）可以深刻影响分子进化的速率[38,46]。与其他核桃科物种相比，在PSMC图中，R.chiliantha的有效种群规模在5.0至0.5MYA期间急剧下降（图4a）。此外，在2.0到0.1MYA的大部分时间内，R.chiliantha的Ne都小于其他物种的Ne（图4a）。我们还使用全基因组杂合性来计算每个物种的长期平均值。2，其中μ是每年的突变率，g是世代时间，假设为30年（方法）。在所有分析的物种中，Rhoipteleachiliantha的Ne值最低（8,900），其他物种的Ne值在10,700至25,700之间（补充表4）。我们比较了谱系之间的Ka/Ks，发现R.chiliantha和其他胡桃科物种之间存在显着差异（图4b和补充图11）。总而言之，这些结果表明R.chiliantha的有效种群规模较小，如果有的话，这将导致更高的替代率，而不是这里发现的较慢的替代率。图3 Rhoipteleachiliantha与其他胡桃科物种共享核桃全基因组重复(WGD)，但在多倍化事件后表现出慢得多的进化。显示R.chiliantha和J.regia的同源染色体的点图。同线区块的颜色表示其Ks中值。一个物种中的每对同源染色体与另一个物种中的一对同源染色体表现出共线关系。点颜色区分直系同源（黄色：较小的K）和旁系同源（绿色：较大的K）关系，指出两个物种共享的WGD。b三种基因家族的系统发育（见正文）。顶部、中部和底部的系统发育树分别来自1型、2型和3型基因家族。对于括号中显示的分数，分母是Bootstrap值≥80的基因树总数（2655棵树），分子表示显示相应拓扑的Bootstrap值≥80的基因树的数量。左侧的系统发育树（2464棵树）支持共享WGD，右侧的系统发育树（123棵树）支持独立WGD。(b)和(c)中的红星标记了WGD事件，红色圆圈表示物种形成事件。c相对速率测试用于估计不同物种的进化速率。Rch1和Rch2是R.chiliantha的WGD衍生同源物，Jre1和Jre2是J.regia的WGD衍生同源物；Rch1-Jre1和Rch2-Jre2是由于物种形成而形成的直系同源基因对。Ks(S1-Rch1)是自物种形成事件（标记为S1）以来Rch1中每个同义位点的预期取代数，其他Ks估计值的定义类似。d使用双尾两样本Wilcoxon秩和检验（Mann-Whitney检验）比较R.chiliantha和J.regia的分子进化速率。***P<0.001。源数据作为源数据文件提供。与其他核桃科相比，R.chiliantha中与DNA修复和重组相关的基因更多Rhoipteleachiliantha有323个被鉴定为参与DNA修复和重组的基因，而其他核桃科有271至302个此类基因（补充数据2和补充数据3）。这些数字是通过将所有研究物种的DNA修复和重组基因分配给260个直系群（基因家族）来推断的（补充数据3）。R.chiliantha和其他核桃科物种之间每个直系类群的基因数量显着不同（图4c）。我们还检查了R.chiliantha和其他每个物种之间的DNA修复和重组基因的不对称模式。共有43个直系类群，其中R.chiliantha有2个基因拷贝，J.regia最多有1个基因拷贝；相反，只有7个直系群，其中J.regia有两个基因拷贝，而R.chiliantha最多有一个拷贝。也就是说，R.chiliantha和J.regia的比较中存在43:7的强烈不对称性。类似的数字为43:17、37:14、34:10、49:9和43:11通过将R.chiliantha与E.roxburghiana、C.illinoinensis、J.mandshurica、J.microcarpa和P.strobilacea进行比较而获得。还有27个直系类群，其中R.chiliantha有两个拷贝，而其他核桃科物种只有一个。在这27个直系群中的26个中，R.chiliantha中的两个拷贝被确定为源自核桃WGD，而在一个直系群中，这两个拷贝是转座衍生的。都是54R。27个直系群中的chiliantha基因拷贝在花蕾中表达，其中45个（83%）也在叶子中表达，这表明大多数（如果不是全部）都具有真正的功能（补充数据4）。相比之下，我们发现只有5个直系群，其中R.chiliantha最多有1个基因拷贝，而其他核桃科物种至少有1个基因拷贝(补充数据3)，并且6个基因拷贝在J.chiliantha的柔荑花序、雌花和叶中表达。(补充表5）。图4 与DNA修复和重组相关的基因家族中的有效群体规模和基因数量。aPSMC估计的七种核桃科物种的有效种群规模。bR.chiliantha（红色）和J.regia（蓝色）的Ka/Ks分布。Ka/Ks被评估为焦点物种及其与其他物种的共同祖先之间每个同义位点的预期替换数。对Ka/Ks分布进行双尾双样本Wilcoxon秩和检验（Mann-Whitney检验）。***P<0.001。c与七个核桃科物种DNA修复和重组相关的基因家族中基因数量的小提琴图。水平线上的数字代表R.chiliantha与其他物种之间的双尾双样本Wilcoxon秩和检验的P值。源数据作为源数据文件提供。我们对核桃WGD中的R.chiliantha中保留了两个拷贝的基因进行了KEGG富集分析，但在其他核桃科物种中最多只保留了一个拷贝。这些基因在转录偶联修复（RPB1）和碱基切除修复（Aprataxin）等多种功能中显着富集（补充图12）。最近的两项研究表明，转录偶联修复在DNA修复中发挥的作用比之前认为的要大得多[47,48]。对于七个核桃科物种，我们根据整个染色体或种内共线块上保留的祖先基因的数量对每对同源染色体进行亚基因组分配[49-51]（补充说明1；补充数据5和补充数据6)。亚基因组分配后，我们使用三种方法来推断核桃家族中的亲本谱系和亚基因组关系，即微同线、局部基因含量和DNA序列比对（同线块的数量、矩阵长度和基因家族编号在补充表中给出）6).基于微同线性和基于基因内容的方法（通过划分同源染色体来分配亚基因组）都产生了拓扑，其中七个隐性亚基因组和杨梅形成了一个进化枝，该进化枝与七个显性亚基因组（染色体）形成的进化枝是姐妹。具有更多祖先基因的在这里被称为“显性”，那些具有较少祖先基因的被称为“隐性”）。MCScanX[52]中的不同参数设置，即A5G25、A10G25或A15G25，其中'A'是共线基因对的最小数量，'G'是相邻块之间插入基因的最大数量，没有改变拓扑结构，与我们的假设一致，即核桃科是异源多倍体起源的（图5a、b和补充图13）。当显性和隐性亚基因组通过划分种内共线块来分配时（允许我们包括一个额外的物种，Pterocaryastenoptera，它是在支架水平上组装的），A5G25参数设置产生了一个拓扑，其中子基因组形成了一个分支，是M.rubra的姐妹（补充图14a），而A10G25和A15G25设置产生的拓扑结构中隐性亚基因组和杨梅形成一个分支（补充图14b、c）。在所有分析中，无论亚基因组是通过划分同源染色体还是种内共线块来分配，化香总是与黄杞归为一组，具有较高的bootstrap支持（图5a、b、补充图13和补充图14a-d）。图5 从全基因组微同线性[Syn-MRL]、基因内容和序列比对获得的核桃科系统发育树。由同源染色体分配的亚基因组(七个核桃科物种和五个外群的显性和隐性亚基因组；显性亚基因组（D）；隐性亚基因组(R)。核桃科物种有Caryaillinoinensis(Cil);Engelhardiaroxburghiana(Ero);Juglansmandshurica(Jma);Juglansmicrocarpa(Jmi)、Juglansregia(Jre);Platycaryastrobilacea(Pst);Rhoipteleachiliantha(Rch);外群是Betulapendula(Bpe);Carpinusfangiana(Cfa);Myricarubra(Mru);Ostryopsisdavidiana(Oda);Quercuslobata(Qlo)。该面板显示了Syn-MRL在A5G25下推断的系统发育（A：调用共线块所需的锚对的最小数量，G：共线块中两个（相邻）锚对之间的最大干预基因数量），（a）来自微同线性包括显性和隐性亚基因组、仅显性亚基因组或仅隐性亚基因组，(b)来自基因存在/不存在或(c)DNA序列比对，包括两个亚基因组、仅显性亚基因组或仅隐性亚基因组。(a,b)中显示了每个节点的超快引导(UFBoot)支持(%)，(c)中显示了每个内部节点的局部后验概率。在通过序列比对方法获得的树中，通过划分同源染色体或种内共线块来分配亚基因组，核桃科的显性和隐性亚基因组形成了一个与M.rubra姐妹的进化枝（图5c和补充图14e））和黄杞和化香相互分离（图5c和补充图14e）。我们还使用GRAMPA53来确定多倍体进化枝及其亲本谱系最可能的位置，而无需先验的亚基因组分配。GRAMPA需要一棵物种树作为开始，我们用ASTRAL-Pro54推断出它。GRAMPA推断的最佳多标记树(MULtree)与我们的子基因组感知序列比对方法获得的树一致（补充图15）。从基因组结构数据中获得的单型亚洲化香属的特殊位置——无论使用的共线块大小或亚基因组如何，其保持相同（图5a、b、补充图13和补充图14a-d）——无论是在我们的研究中（图5c和补充图14e）还是在过去20年[9,10,14–16,18,19]中，都没有使用DNA比对来恢复。这意味着化香的放置在基于基因组结构的基因组中系统发育或基于DNA比对的系统发育一定是错误的。由于目前尚不清楚真正的系统发育，是否有形态学证据支持化香属的一种或另一种放置？Wing和Hickey[20]对核桃科化石和现存物种中的30个性状进行了分析，得出了化香属/黄杞属及其中美洲近缘种、Alfaroa属（有5个中美洲物种）和Oreomunnea属（有2个墨西哥物种），但两者均不包括在内这里。在东南亚，黄杞包括七种，其中两种，越南、老挝、缅甸和中国西南部的E.roxburgiana和中国东南部的E.fenzelii[55]，有时被认为是单独的Alfaropsis属（补充图16显示了核ITS系统发育整个进化枝）。虽然Wing和Hickey[20]没有公布他们的数据矩阵，但他们的文本强调了化石和现存化香属和黄杞共有的四到七叶结构特征以及花粉和圆锥花序形态。此外，两个属的坚果壁均较薄，无空隙，而核桃属、山核桃属和枫杨属的坚果壁有空隙[22]。受Wing和Hickey研究的影响，Manchester[25]的核桃科直观系统发育显示出相同的分组（我们的图1a）。还有四个最新的形态学数据矩阵可供使用，我们在本研究中重新分析了它们。Manos和Stone[9]为40个现存分类群编码了64个字符，其形态树使化香属处于黄芪分支和核桃分支之间的三分体中（Manos等人[10]：图3A）。Hermsen和Gandolfo[24]修正并扩大了Manos和Stone[9]的矩阵，为28个现存类群和6个灭绝类群编码了64个字符。他们的数据产生了（Platycarya、Engelhardia、Alfaroa、Oreomunnea）进化枝，但Rhoiptelea的位置尚未解决（我们的补充图17显示了根据他们的数据生成的Neighbour-Net），可能是因为他们未能包括最古老的Rhoiptelea化石，一个85年前的果结构，是所有胡桃科植物中最古老的化石[45]。Larson-Johnson[14]对来自37个现存分类单元和27个已灭绝分类单元的89个字符进行了编码，Zhang等人[19]。编码了47个现存类群和113个已灭绝类群的73个字符。尽管在Larson-Johnson[14]的系统发育史中发现了化香类群与黄杞属的化石（图4），但这两个矩阵都未能明确地确定化香属的位置。因此，形态特征与此处发现的化香的位置并不矛盾，并且基于Wing和Hickey[20]可能支持这一点（补充图17）。基于DNA比对的树都将黄杞和化香放置得很远（图5c和补充图14e）。对此的解释可能是核桃科之间的后多倍体基因流（如核桃属[26,27]中所证明的），这将导致不一致的系统发育信号。不需要预先指定的物种树的基因树不和谐测试使用多物种合并(MSC)模型下的基因树四重体来计算四重体计数一致性因子(qcCF)[56,57]。我们对每组四个分类单元进行了这样的分析，以评估MSC模型下的预期qcCF与数据的拟合程度。通过同源染色体分配亚基因组，我们对七个类群的5,882个基因树进行了MSCquartets分析。应用Holm-Bonferroni方法对多重检验进行调整，35项检验中的12项（约34.29%）拒绝原假设，转而支持显着的不一致。特别值得注意的是，所有四重奏，每个由黄杞属、化香属、山核桃和核桃属物种之一组成（补充数据7），都进行了重要的测试，表明它们之间可能存在古代杂交。如果亚基因组由种内共线块分配，则获得类似的结果（参见补充图18和补充数据8）。如上所述，我们缺乏中美洲Alfaroa属（5种）或Oreomunnea属（2种）和第二种亚洲黄杞属的代表性基因组，最好是该属的模式种E.spicata（补充图16）。它们的包含可能会影响化香与黄杞分支的精确附着。我们的Ks分析、基因组同线性（共线性）和系统发生学结果共同证实了核桃科共同祖先（包括Rhoiptelea）中发生了WGD事件（图3和补充数据1）。我们对核桃科基础上异源多倍体事件的假设（最初从细胞遗传学[7,8]推断）得到了系统发育的支持，该系统发育显示隐性亚基因组分支是杨梅的姐妹（图5a-b，补充图13和补充图14b-c)，与Xiao等人的结果一致[6]。他们以序列相似性为标准，发现山核桃的一个亚基因组与杨梅表现出更高的平均同一性。令人放心的是，在之前对核桃和山核桃基因组的研究[3、5、6]以及我们的研究中，已经检测到两个亚基因组之间重复基因的不成比例丢失（即，有偏分馏）（图6和补充图19-24）。隐性亚基因组进化枝是杨梅属的姐妹，这可能意味着核桃科的一个亲本是与杨梅科密切相关的物种，而另一个亲本已灭绝或未采样。然而，显性亚基因组和隐性亚基因组的基于比对的树（图5c）显示核桃科形成了一个与杨梅属姐妹的进化枝，这意味着核桃科的亲本要么嵌套在杨梅科内，要么与杨梅科有一个最近的共同祖先。杨梅科。未来的工作应该对杨梅科的额外基因组进行采样，最好是Canacomyricamonticola、Comptoniaperegrina或MyricaGale。在Rhoiptelea以及其他核桃科物种中检测到两个Ks峰，对应于核心真双子叶植物共享的核桃WGD和γWGT（图2b）。R.chiliantha对应于核桃WGD的第一个Ks峰比其他6个科内成员小2.1至2.8倍（图2b）。正如Doyle和Egan58所指出的，异源多倍体中同源物的Ks峰会高估其分子进化速率，因为亚基因组会在两个祖细胞之间的物种形成点发生分歧，而不是在杂交（和基因组加倍）事件本身。因此，来自WGD的同源物的K充其量是后代谱系替代率的粗略估计（参见Thomas等人[53]：图1）。然而，更直接的估计（图3c）得到了类似的结果：在胡桃样WGD之后，R.chiliantha的速率比其他物种慢2.6至3.5倍。其他证据支持R.chiliantha的替代率极慢。BUSCO评估结果（补充图25）显示，在R.chiliantha中，embryophyta_odb10数据库中24%的通用单拷贝直向同源物是完整的，并且是重复的，比其他核桃科物种多3-4倍（5.1-6.8%））；这些重复的BUSCO组中的大多数都源自核桃WGD（补充数据9）。基于沿每条Q.lobata染色体的100个基因滑动窗口，Rhoiptelea的保留直系同源基因的百分比高于其他六种核桃科物种的显性亚基因组和隐性亚基因组（图6和补充图19-24）。我们还发现R.chiliantha保留了更多的WGD基因重复（补充图26和补充表7），表明Rhoiptelea基因重复中积累的有害突变较少。然而，我们不能排除其他假设，例如WGD59后基因的缺失差异，需要进一步研究。物种的种群规模会影响选择和漂移的强度，从而导致种群中近中性突变的固定概率不同。在小群体中，弱有害突变可能会保留为有效的中性等位基因，因此通过负选择暴露和从基因库中删除的弱有害突变会减少，这应该会导致分子进化的高速率[34-36]。然而，我们的PSMC、Ka/Ks比值和杂合性结果（图4和补充图11）表明，无论以何种方式推断，R.chiliantha的Ne都较小，因此不太可能是其Ne解释了其较慢的替代率。相反，与其他核桃科植物相比，R.chiliantha中参与DNA修复和重组的基因数量更多（图4c，这是七个物种中与DNA修复和重组相关的基因家族中基因数量的小提琴图，补充材料数据3)以及Rhoiptelea中富含转录偶联修复(RPB1)和碱基切除修复(Aprataxin)的核桃WGD基因表明，更有效的DNA修复是核桃WGD后R.chiliantha替代率降低的主要原因。在图2b中，对应于核桃WGD的Ks峰值范围为0.17至0.48，而对应于γWGT的Ks值范围为1.70至2.01。假设核桃WGD发生在大约85MYA45，简单的计算告诉我们，在γ-WGT之后但在核桃WGD之前的时间段（即从~120到~85MYA），核桃科的进化率为7.7到22比胡桃样WGD后（即从85MYA至今的时间段）高出几倍。换句话说，分子在核桃科WGD进化后，核桃科的生长速度明显减慢。Smith和Donoghue[32]首先提出了核心真双子叶植物木质分支的急剧放缓，并已在多项基因组研究中进行了量化。例如，水杨酸WGD（柳科）的Ks峰值出现在60-65MYA[60,61]范围内，范围为0.34至0.56，而柳树所属的核心真双子叶植物的γ-WGT的Ks峰值约为2.5[31]，表明放缓了5.0倍。在蔷薇科亚科Pomoideae（包含苹果和其他果树）中，PomoideaeWGD的Ks峰值在48-50MYA[62,63]范围内为0.27至0.39，而γ-WGT的Ks峰值在2.5[31]左右，表明放缓了4.5倍。忽视这种急剧的速率放缓会导致对进化事件的不可靠估计[64]，例如核心真双子叶植物共享的γ-WGT与核桃类WGD[6]的时间的校准。总之，化香与黄杞的位置对于胡桃仁壁的进化以及该科的生物地理学具有重要意义。微同线性和基于基因内容的系统发育方法包含迄今为止被低估的系统发育信号。该信号对于研究异源多倍体生物的系统发育可能特别有价值，因为异源多倍体同时涉及杂交和物种形成。对于此类谱系，基因组结构允许区分共存的同源染色体并单独使用它们来推断生物体的祖先。更重要的是，任何能够改变基因顺序和/或内容的大型基因组结构变异可能很少在物种之间水平传播，从而使得基于基因顺序或内容的树推断对基因渗入更加稳健。千叶菊（中国云南省西稠县，北纬23°22′42.73′′N，104°47′17.21′′E）和黄杞（中国贵州省望谟县，25°15′30.62′N，105°5740.7′′E）的鲜叶被收集用于提取基因组DNA并测序。各物种的永久凭证已存入北京师范大学植物标本馆（张BNU20180707-4、曹BNU20200818-1）。在IlluminaHiSeqXTen测序平台上分别对Rhoipteleachiliantha和E.roxburghiana进行44.28和103Gb150bp双端测序测序，然后使用GenomeScope2.0[65]通过Jellyfishv2.3[66]构建的k-mer计数直方图评估其基因组聚体长度为17。Rhoiptelea和E.roxburghiana的基因组大小分别为~410Mbp和~880Mbp（补充图1）。对于R.chiliantha，总共使用43Gb(~105×)PacBio单分子长读数和103Gb(~108×)插入片段大小为350bp的Illumina短读数进行初始组装和后续校正。此外，来自Hi-C文库的总共40.86Gb(~100×)原始数据被用于染色体规模的基因组组装。对于E.roxburghiana，在组装和后续校正中总共使用了71Gb（~75×）PacBio单分子长读数和103Gb（~108×）插入大小为350bp的Illumina短读数。来自Hi-C文库的总共305Gb(~347×)原始数据用于染色体规模的基因组组装。每个基因组组装的详细信息在补充说明2中给出。为了代表核桃科的主要谱系（如图1所示），我们从公共数据库下载了Carya[67]、Juglans[3]和Platycarya[27]的五个染色体水平基因组。我们还下载了杨梅科（Myricaceae）和其他除了γ-WGT外没有经历进一步多倍体事件的山毛榉的基因组，作为外群[68-72]用于后续分析。杨梅科有3属50种，是胡桃科的姐妹科。为了注释R.chiliantha和E.roxburghiana的基因组，使用了基于同源性的推理、从头开始预测和RNA测序转录(RNA-seq)的组合（详细信息参见补充说明3）。最终的基因集使用京都基因和基因组百科全书(KEGG)自动注释服务器(https://www./tools/kaas/)73进行功能注释，以执行KEGGOrthology(KO)注释。 我们使用来自七个物种（C.illinoinensis、E.roxburghiana、J.mandshurica、J.microcarpa、J.regia、P.strobilacea和R.chiliantha）的蛋白质序列和注释信息来检测并确认WGD。BLASTP[74]（e<10−10和top5匹配）用于搜索每个物种内蛋白质序列的所有潜在同源基因对，并通过MCScanX[52]在默认设置下识别共线基因对。我们使用TBtools[75]分别可视化R.chiliantha和E.roxburghiana16条染色体的共线块。利用DupGen_finder[31]中7个物种的种内共线基因对信息以及7个物种与Quercuslobate[69]之间的种间共线基因对信息,我们识别了5种基因复制模式,包括全基因组复制(WGD)、串联重复（TD）、近端重复（PD）、转置重复（TRD）和分散重复（DSD）。对于WGD的所有共线基因对（锚定基因），我们使用KaKs_Calculator2.0[77]中的GammaMYN算法[76]估计每个同义位点的同义替换(Ks)、每个非同义位点的非同义替换(Ka)和Ka/Ks。使用Perl脚本计算共线基因对之间的四重同义第三密码子颠倒率(4DTv)(https://github.com/JinfengChen/Scripts/blob/master/FFgenome/03.evolution/distance_kaks_4dtv/bin/calculate_4DTV_Correction。pl）。为了确保物种内每个共线块中Ks和4DTv的独立性，我们遵循Schnable等人的方法[49]。并估计每个共线块的中值Ks和4DTv值。Ks值等于0或大于5[31]以及4DTv距离等于0或大于1被排除。R包ggplot2[78]用于绘制共线基因对的Ks和4DTv分布。 每个采样的核桃科物种的WGD事件的推断和替代率的估计根据共线基因对的Ks分布，我们发现R.chiliantha的第一个峰值位于0.17，远小于其他核桃科的0.36至0.48（图2b）。为了解释这种差异，我们测试了两个假设：(1)R.chiliantha独立经历了不同的WGD，或(2)R.chiliantha共享相同的核桃WGD，但在WGD后进化速度要慢得多。在后一种假设下，可以使用为相对速率测试开发的程序来估计每个物种的进化速率[79]。为了测试WGD事件是否共享，我们使用了基因组间同线性分析[42,43]（图3a）和Pfeil等人的“多基因家族”方法[44]。（图3b）。对于基因组间同线性分析[42,43]，我们使用BLASTP（e<10−10，top5匹配）和MCScanX来确定R.chiliantha和其余核桃科物种之间的基因组共线块。使用KaKs_Calculator2.0中的Gamma-MYN算法计算共线基因的Ks后，我们使用WGDI[80]可视化R.chiliantha和其余核桃科物种之间的共线块的中位Ks值。如果R.chiliantha和其他核桃科共有一个WGD事件，那么它应该先于它们彼此分歧；物种之间的同源染色体应表现出不同的分歧，WGD共线块的中位Ks值应大于R.chiliantha和其余核桃科物种之间的直系同源共线块的Ks值（图3a）。然而，如果R.chiliantha和其他物种经历了独立的WGD，这些WGD事件将在它们彼此分歧之后发生，因此一个物种中的每条染色体应与其他物种中的任何同源染色体距离相等，从而产生相似的中位数WGD和正交共线块的Ks值。对于“多基因家族”方法[44]，我们使用了R.chiliantha或其他胡桃科物种（C.illinoinensis、E.roxburghiana、J.mandshurica、J.microcarpa、J.regia、P..strobilacea)和至少一份Quercuslobata(参见补充说明4和补充图27中的详细信息)。在这组基因家族中，我们区分了三种类型：（i）R.chiliantha和其他核桃科都有两个基因拷贝的家族；(ii)R.chiliantha有两个拷贝而其他物种只有一个的科；(iii)R.chiliantha只有一个拷贝而其他物种有两个拷贝的科。如果R.chiliantha和其他核桃科(SP)物种共享WGD事件，则基因谱系树的拓扑结构应为(Q.lobata,((SP,R.chiliantha),(SP,R.chiliantha)))或(Q.lobata,((SP,R.chiliantha),SP))或(Q.lobata,((SP,R.chiliantha),R.chiliantha))（图3b中的左侧部分）。相比之下，如果R.chiliantha和其他核桃科(SP)物种的WGD独立发生，则拓扑结构应为(Q.lobata,((SP,SP),(R.chiliantha,R.chiliantha)))或(Q.lobata,((SP,SP),R.chiliantha))或(Q.lobata,((R.chiliantha,R.chiliantha),SP))（图3b中的右侧部分）。在确定R.chiliantha和其他物种共享WGD事件后，我们使用相对速率测试[79]来估计每个物种的进化速率（Ks）（图3c，补充图27）。平均进化率的计算公式为： 其中r是每年每个同义位点的核苷酸取代率，Ks是共同祖先和焦点物种之间每个同义位点的预期取代数，t是分歧时间[64]。
为了确定核桃科的七个物种是否经历了种群收缩/扩张，我们使用成对顺序马尔可夫合并（PSMC；v.0.6.5-r67）方法[81]从一个二倍体个体的全基因组信息中推断出Ne的变化。除E.roxburghiana和R.chiliantha外，本研究中的全基因组测序数据均来自于其他地方使用的已发表的核桃科基因组（补充表1）。我们使用BWAv.0.7.15的BWA-mem算法[82]将七个个体的修剪读数与他们各自的参考基因组进行比对，并使用SENTIEONDNAseq软件包v.201808.08[83]对比对序列进行排序、删除PCR重复项并重新比对插入缺失每个人。然后，我们使用获得的二进制比对图（BAM）文件和SAMTOOLSv.1.2[84]来执行变体调用。使用质量调节器-C设置为50、最小映射质量为20来过滤SNP。在为PSMC准备输入文件时，我们将最小深度设置为平均基因组深度的三分之一，将最大深度设置为平均值的两倍vcfutils.plvcf2fq(https://github.com/lh3/samtools/blob/master/bcftools/vcfutils.pl)管道中的基因组深度。替代率设定为本研究估计的值，世代时间假设为30年[85]。我们使用ANGSD[86]计算每个物种的一个个体的杂合度。每个个体的BAM文件作为ANGSD的输入，参数设置为“-C50-minQ20-minMapQ30”。全基因组杂合度用于估计每个物种的长期平均值：根据近中性理论，有效种群规模(Ne)较小的种群选择效率较低46。因此，Ne较小的物种可能会积累更多的非同义突变，导致较大的Ka/Ks。我们比较了R.chiliantha和其他核桃科物种之间Ka/Ks的分布，以了解两个物种中具有相似Ks值的基因。四种主要的DNA修复途径：碱基切除修复(BER)、核苷酸切除修复(NER)、同源重组修复(HR)和非同源末端连接(NHEJ)在真核生物中高度保守[87]。我们使用KAAS（KEGG自动注释服务器：http://www./kegg/kaas/）对7个核桃科物种的蛋白质进行注释，然后统计与DNA修复和重组相关的基因数量（KEGGOntology：ko03400）。收集了R.chiliantha的两个组织（叶和花蕾）的11个个体进行RNA测序。我们还下载了J.regia2的叶、柔荑花和雌花的双端RNA测序数据。使用Trimmomaticv0.32[88]修剪低质量碱基或读数后，我们使用Hisatv.2.1.0[89]将修剪后的读数与参考基因组对齐，并使用htseq-countv.0.11.2[90]来计算读数数量。使用R脚本计算每个基因每百万映射读数中每千碱基转录本的片段数(FPKM)。通过这种方式，我们使用RNA-seq数据来验证J.regia和R.chiliantha中表达的基因有多少与DNA修复和重组相关。为了推断核桃WGD中R.chiliantha保留2个拷贝的基因的功能，而其他核桃科物种最多只有1个拷贝，我们使用R包clusterProfile[91]进行KEGG富集分析。 为了推断核桃属WGD的亲本物种，我们使用了7个核桃科物种和5个在γ-WGT之外没有经历过进一步WGD的山毛豆目代表作为外类群，即Quercuslobate[69]、Betulapendula[70]、Carpinusfangiana[68]、Ostryopsisdavidiana[71]、和Myricarubra[72]。根据保留的祖先基因的数量[49-51]，对每对同源染色体（或种内共线块，详细信息见补充说明1）进行亚基因组分配。我们使用Q.lobata和七种核桃科中的每一种之间的相互最佳命中（RBH）来获得保留的祖先基因的列表。我们将具有更多祖先基因的染色体表示为“显性”，将具有较少祖先基因的染色体表示为“隐性”。这种分配亚基因组的方法产生与Zhu等人相同的结果[3]。根据所有基因的数量使用了略有不同的方法，Xiao等人[6]。使用进化距离来分配亚基因组。亚基因组分配后，我们在三个数据集中应用了三种系统发育推断方法（基于ASTRAL-Pro[54]比对、基于全基因组微同线性和基于局部基因内容），通过包括(i)显性和隐性来推断核桃科的系统发育。所有胡桃科和上述五个外类群的亚基因组；(ii)仅每种核桃科及其外类群的主要亚基因组；(iii)仅每个胡桃科和外类群的隐性亚基因组。选择每个基因具有最长转录本的蛋白质序列进行进一步分析。对于DNA序列比对方法，我们使用OrthoFinderv.2.4.0[92]来获取每个数据集的基因家族。为了避免太大的基因家族[93]引入的多重序列比对错误，我们将基因家族的大小保持在小于100个基因以进行下游分析。使用MAFFTv7.475[94]比对每个基因家族的蛋白质序列，然后使用PAL2NALv.14[95]转换为CDS比对。我们使用IQ-TREEv2.1.2构建基因家谱树，然后将其用作ASTRAL-Pro54的输入数据来推断物种树。补充表6提供了三个数据集中基因家族的详细信息。对于基于全基因组微同线性的系统发育推断，我们使用具有似然性的同线性矩阵表示（Syn-MRL）管道[30]从上述三个数据集中重建系统发育树，即显性和隐性亚基因组（数据集1），仅显性亚基因组（数据集2），并且仅限隐性亚基因组（数据集3）。继赵等人之后[30]，BLASTP[74]用于搜索所有潜在的种间和种内同源基因对（e<10−10和top5匹配）。使用MCScanX[52]按照A5G25的推荐参数检测种间和种内的同线性块（A：调用共线块所需的最小锚对数，G：两个（相邻）锚对之间的最大干预基因数在共线块中）。我们还探索了A10G25和A15G25参数设置，这将从分析中删除具有少于10或15个共线基因对的较小共线块。在所得的同线性网络中，节点是共线块中的共线基因，边用于连接共线基因对。使用Infomap算法v.1.6.0(https://github.com/mapequation/infomap)，我们在映射方程框架内检测到同线性簇，并通过十次试验设置两级划分模式(--clu-N10−2)。对于具有A5G25设置的数据集1，整个同线性网络总结了67,954个成对同线性块的信息，并包含270,276个节点和1,976,821个边。使用Infomap算法v.1.6.0将整个网络分配给19,849个集群。这些簇被转换成二进制存在-不存在数据矩阵，其中行和列分别代表物种和簇。继Pett等人之后，我们删除了仅包含一个（子）基因组的簇[41]。由此产生的19行×19,202列矩阵用于在Mk R FO模型下使用IQ-TREEv2.1.2推断系统发育。补充表6提供了数据集2和3中具有A10G25和A15G25设置的矩阵的详细信息。对于基于基因内容的系统发育推断[41]，我们遵循Zhao等人的方法[30]。使用OrthoFinderv.2.4.0分配三个数据集的基因家族，然后使用基因存在/不存在矩阵使用IQ-TREEv2.1.2构建系统发育树。对于每个数据集，基因存在/不存在矩阵分别有19行×17,149列、12行×17,926列和12行×16,894列。由于GRAMPA[53]可以在多倍体存在的情况下推断基因重复和丢失，并确定多倍体进化枝及其亲本谱系最可能的位置，因此我们还使用此方法来推断核桃WGD的亲本谱系，而无需执行先验亚基因组分配。首先，我们使用OrthoFinderv.2.4.0[92]获取7个核桃科物种和5个外群的基因家族，并使用IQ-TREEv2.1.2构建基因树。然后，我们使用ASTRAL-Pro从基因树推断物种树。我们将核桃科的七种物种指定为可能的多倍体谱系（H1节点），GRAMPA搜索了第二亲本谱系（H2节点）的可能位置。最后，GRAMPA返回所考虑的每个多标记树（MUL树）（以及原始单标记物种树）的协调分数。具有最小协调分数的MUL树是H2节点的最佳放置。文章标题：Genome structure-based Juglandaceae phylogenies contradict alignment-based phylogenies and substitution rates vary with DNA repair genes 作者信息：Ya-Mei Ding, Xiao-Xu Pang, Yu Cao, Wei-Ping Zhang, Susanne S. Renner*, Da-Yong Zhang* & Wei-Ning Bai* 原文链接：https:///10.1038/s41467-023-36247-z 文内图片及封面图片来源原文
|
