2023年8月1日，Horticulture Research在线发表了中国科学院华南植物园、热带作物生物育种全国实验室-中国农业科学院深圳农业基因组研究所和广西民族大学合作完成的题为“Telomere-to-telomere haplotype-resolved reference genome reveals subgenome divergence and disease resistance in triploid Cavendish banana”的研究论文。

香蕉是世界上最重要的作物之一。卡文迪许（Cavendish）型香蕉，源自单种源的Musa acuminata（AAA），约占全球香蕉产量的一半，对人类社会具有重要意义。然而，直到现在，香蕉栽培品种仍没有高质量的单体型基因组参考序列。在这里，我们报道了'Baxijiao'（巴西蕉）的端到端（T2T）和单体型基因组参考序列，包括三个单倍体组装。这三个单倍体组装的大小分别估计为477.16 Mb、477.18 Mb和469.57 Mb。尽管具有单种源起源，但这三个单倍体组装显示出很大的差异，序列一致性水平较低。在染色体1、4和7之间发现了几个大的相互易位。检测到了一些可能影响水果品质和香气的基因家族扩张，例如属于蔗糖/二糖/寡糖降解过程、蔗糖代谢过程、淀粉代谢过程和芳香化合物生物合成过程的基因家族。此外，观察到与雄蕊和花粉发育相关的基因家族扩张，这可能与卡文迪许品种的无核果和不育性有关。最后，在巴西蕉中发现的抗性基因较M. acuminata较少，尤其是在染色体3和10的基因簇中，为研究香蕉抗病性的分子分析提供了潜在目标。因此，这个T2T单体型基因组参考序列将成为生物学研究、分子育种和香蕉遗传改良的宝贵遗传资源。

我们通过整合PacBio HiFi序列、Oxford Nanopore序列和Hi-C序列生成了巴西蕉基因组。总共，我们分别获得了102.06 Gb（总基因组大小的约72倍覆盖）的PacBio HiFi reads，N50长度超过15 kb，在ONT reads中获得了48.00 Gb（约34倍覆盖），N50为80 kb，并且获得了156.14 Gb（约110倍覆盖）的Hi-C数据。根据GenomeScope 2.0的分析结果（基于k-mer为21），我们估计巴西蕉的单倍型基因组大小约为457.31 Mb，杂和率为2.855%。为了实现T2T和单倍型分离的基因组，我们采用了一种组装方案，该方案可以将HiFi和Hi-C reads结合起来组装相位基因组，并利用ONT reads填补缺口。基于HiFi和Hi-C reads，我们使用hifiasm构建了一个相位二倍体组装，生成581.77 Mb和959.85 Mb的主要contig，N50分别为14.95 Mb和16.15 Mb。在M. acuminata ssp. malaccensis的T2T基因组的辅助下，我们进一步对contig进行了排序、定向和分组。然后，利用Hi-C reads将一些短contig锚定并删除，同时将较大的contig分割成两组。我们使用NextDenovo工具构建了ONT组装，并用该组装填补了HiFi组装参考中的缺口。在填补了缺口后，生成了一个几乎完整且具有相位的参考基因组，命名为BXJ，包含33个伪染色体（每个单倍体组装有11个），总长度为1.42 Gb，单倍体组装大小为469.57-477.18 Mb（详见图1A和表1。这三个单倍体组装分别命名为BXJ1、BXJ2和BXJ3。每个单倍体组装中最长的染色体是第8号染色体，分别为BXJ1、BXJ2和BXJ3的大小分别为53.00 Mb、54.63 Mb和52.65 Mb。最短的染色体在BXJ1和其他两个单倍体组装之间略有差异；具体而言，BXJ1中最短的染色体是第2号染色体（33.52 Mb），而BXJ2（33.79 Mb）和BXJ3（32.90 Mb）中最短的染色体都是第11号。我们在BXJ基因组中使用100 kb间隔计算了基因密度、转座子密度、Copia密度、Gypsy密度、GC含量和潜在的着丝粒串联重复序列密度(详见图1A)。这是第一个用于三倍体香蕉的T2T和单倍型分离的参考基因组。

在三个单倍体组装中，分别鉴定出了256.48 Mb（53.76%）、258.29 Mb（54.14%）和258.84 Mb（55.13%）的重复区域，与其他香蕉基因组（例如M. acuminata: 52.62% ；M. beccarii: 55.99% ）相似。其中，最丰富的重复序列是长末端重复序列（LTRs），分别占BXJ1、BXJ2和BXJ3组装的36.36%、35.44%和33.53%。从BXJ1、BXJ2和BXJ3共预测到37,185、37,241和37,178个高置信度编码蛋白基因(表1)，略多于最新的M. acuminata双单倍体基因组v4版本（下称MAv4）预测的基因数（36,979个）[7]。BXJ1、BXJ2和BXJ3中预测的基因总长度分别为142.78 Mb、138.99 Mb和133.28 Mb。这些编码蛋白基因的平均长度分别为3,839 bp、3,732 bp和3,584 bp。使用七碱基重复序列（5'端为CCCATTT，3'端为TTTAGGG）作为查询，我们总共估计出53个端粒，并构建了BXJ1、BXJ2和BXJ3的8、6和6条端粒到端粒染色体。在BXJ1中，鉴定到的端粒长度范围从6,006 bp到33,838 bp，在BXJ2中为2,345 bp到44,905 bp，在BXJ3中为4,494 bp到26,964 bp。在BXJ2中，最长的端粒序列位于Chr10的右端，而最短的端粒位于BXJ2的Chr08右端。着丝粒是染色体中保持两个子染色单体结合并参与有丝分裂和减数分裂运动的特殊核苷酸序列。使用Tandem Repeats Finder ，确定了“巴西焦”基因组中着丝粒的大致位置，其长度分别在BXJ1的924,207 bp到3,415,439 bp之间，在BXJ2的930,592 bp到3,844,332 bp之间，在BXJ3的481,633 bp到4,080,464 bp之间。对BXJ组装的质量和完整性进行了多方面的评估。首先，根据估计，每个单倍体组装的基因组大小（469.57-477.18 Mb）略大于MAv4（约468.82 Mb）；此外，由JuiceBox 可视化的Hi-C数据在所有染色体上呈现出很高的一致性，证明了其排序和定位的准确性(图1B)。其次，将HiFi和Illumina reads对齐到这三个单倍体组装上；每个染色体的平均HiFi reads覆盖度约为200×(图1C)；此外，97.95%、97.66%和96.36%的Illumina reads，以及99.84%、99.80%和99.50%的HiFi reads分别成功对齐到BXJ1、BXJ2和BXJ3(图1D)，表明这个组装对“巴西焦”基因组具有足够高的覆盖度。第三，三个单倍体组装中观察到了较高的LTR组装指数（LAI）得分（19.84、20.65和20.22）(表1)。最后，使用BUSCO [22]评估完整性，结果显示BXJ1、BXJ2和BXJ3的完整序列分别占核心真核生物基因组的97.40%、97.80%和93.80%(表1和图1E)。总体而言，这些结果展示了BXJ参考基因组的高质量、准确性和可靠性。全基因组比较和合成分析BXJ1、BXJ2和BXJ3具有一组相似的基因组特征，如相似的基因组大小和基因数，以及相似的重复序列含量(表1)。然而，合成分析揭示了单倍体组装之间核苷酸序列一致性较低的水平(图1F)，这可能归因于巴西焦亚属不同起源的M. acuminata亚基因组。与此同时，由于进化保守的编码区域，在三个单倍体组装之间观察到了高度的同源关系。基因组一致性分析显示，在三个两两比较（BXJ1 vs BXJ2，BXJ1 vs BXJ3，BXJ2 vs BXJ3）中，29131、27571和28838个转录本分别匹配了31个、37个和37个同位基因块(Fig. 1F)。使用MCScan [23]鉴定了BXJ基因组内的大量结构变异。在Chr01、Chr04和Chr07之间鉴定出了几个大的互换易位(Fig. 2A)。具体而言，观察到了Chr01和Chr04的同源关系，表明在BXJ2中涉及Chr01长臂远端的一个互换易位，以及在BXJ1和BXJ3中涉及Chr04长臂0.47-0.50 Mb区域的互换易位(图S4A，请参见在线补充材料)。此外，还鉴定出了另一个涉及BXJ1和BXJ3中Chr01长臂远端的1.88-3.12 Mb区域以及BXJ2中Chr04长臂的互换易位(图S4B，请参见在线补充材料)。此外，还观察到了Chr01和Chr07之间的同源关系，显示了涉及BXJ1中Chr01长臂远端和BXJ2中Chr07的1.04 Mb区域的互换易位(图S4C，请参见在线补充材料)。

根据序列相似性，在单倍型分辨的BXJ基因组和MAv4之间获得了18,263个等位基因组(表S6，请参见在线补充材料)。这些等位基因均匀地分布在基因组的11条染色体上。大多数等位基因的编码序列在三个单倍体组装中与MAv4中的相似度很高（平均值分别为96.50％、96.59％和96.62％）。此外，通过计算等位基因对之间的Ka/Ks值，发现大多数等位基因经历了纯化选择(Ka/Ks <1，图S5，请参见在线补充材料)。

使用单子叶植物中‘巴西焦亚’和15个其他物种的同源基因进行系统发育分析，其中包括香蕉栽培品种的主要祖先M. acuminata（A基因组）、M. balbisiana（B基因组）、M. schizocarpa（S基因组）和M. troglodytarum（T基因组）。系统发育树表明，‘巴西焦亚’的三个单倍体组装聚集在一起，并且是M. acuminata的姐妹类群(图2B)。在选择的单子叶植物物种的612,509个基因集合中，共鉴定了38,387个直系同源基因家族。基于系统发育树，识别了潜在的扩增和收缩基因家族。共有447个由1,551个基因组成的显著扩张和收缩基因家族特异于‘巴西焦亚’。共117个由694个基因组成的基因家族在‘巴西焦亚’基因组中显著扩展。GO富集分析表明，这些扩张基因在199个生物过程中显著富集(表S7，请参见在线补充材料)。最显著的生物过程与应对胁迫和蔗糖诱导的翻译抑制有关(图2C)。有趣的是，其他一些显著富集的生物过程与碳水化合物来源和花粉发育有关，包括对二糖的反应、对蔗糖的反应、对碳水化合物的反应、花粉接受、花粉识别、花粉-柱头相互作用和花粉萌发。同时，330个由857个基因组成的基因家族显著缩小，富集了231个生物过程(表S8，请参见在线补充材料)。富集的大部分生物过程（30％）与植物生长、形态发生和激素有关，如子叶发育、果实发育、花药发育、花器官形态发生、参与繁殖的发育过程、生长素外流、激素介导的信号通路等。此外，一些富集的生物过程（14％）与对外部生物和非生物因素的应对有关，包括对共生真菌的反应、对细菌的反应、免疫应答激活信号转导、红光或远红光信号传导途径以及对蓝光的细胞响应。在单倍体组装中还鉴定了潜在的扩张和收缩基因家族(表S9，请参见在线补充材料)。有趣的是，一些显著扩展的基因家族在芳香化合物生物合成过程、蔗糖降解过程、二糖降解过程、寡糖降解过程、蔗糖代谢过程和淀粉代谢过程中富集，这些可能对果实品质很重要(图S6和表S10-S12，请参见在线补充材料)。此外，一些扩展的基因家族涉及到花药，尤其是花药壁淀粉层的发育，例如花药形态发生和花药壁淀粉层的形成。

植物基因组通常包含大量的抗病性（R）基因。植物的R蛋白主要分为两类：膜结合型模式识别受体（PRRs）和胞内核苷酸结合亮氨酸重复受体（NLRs）。PRRs可以是受体样激酶（RLKs）或受体样蛋白（RLPs）[24]。在本研究中，我们首先使用RGAugury [25]预测了‘巴西焦亚’和M. acuminata中的抗性基因类似物（RGAs）。与MAv4中的124个相比，在BXJ1、BXJ2和BXJ3中预测到了68、52和69个NBS编码基因(表2)。对于RLK和RLP型的RGAs的预测还表明，与M. acuminata相比，‘巴西焦亚’携带的R基因较少(表2)。MAv4和‘巴西焦亚’基因组中上述预测的RGAs的分布详见补充图S7（请参见在线补充材料）。此外，我们还使用NLR-annotator [26]来鉴定可能与抗病性相关的基因组区域。预测的NLR位点可能是具有完整或部分开放阅读框架的基因，也可能是被赝化的序列在研究基因组中的痕迹[26]。在BXJ1、BXJ2和BXJ3中分别检测到了123、126和140个NLR位点(表2；表S13-S15，请参见在线补充材料)，与MAv4中的128个位点进行比较[7]。NLR位点不均匀地分布在11条染色体上(图3)。一些NLR位点倾向于集群分布。在BXJ基因组中观察到了三个主要的NLR位点群：两个在Chr03中，一个在Chr10中(表3)。值得注意的是，BXJ基因组中大部分NLR位点位于基因间区域（BXJ1：37.4％；BXJ2：59.5％；BXJ3：38.6％）。与MAv4中的122个基因相比，在三个单倍体组装中只预测到了77、51和79个抗性基因，与NBS编码的RGAs的预测结果相似。特别是在Chr03长臂的末端基因簇中（BXJ1：43.02-43.80 Mb；BXJ2：41.36-42.09 Mb；BXJ3：36.85-37.62 Mb），所有NLR位点都被预测为赝化序列(图3和表3)。此外，在Chr10的基因簇中，BXJ2中的所有NLR位点被预测为赝化序列（21.93-22.73 Mb），而BXJ1（19.29-20.01 Mb）和BXJ3（20.93-21.73 Mb）中的大多数NLR位点根据基因注释仍然具有功能性(图3和表3)。总体而言，‘巴西焦亚’中鉴定到的抗性基因比M. acuminata中的要少得多。

为了确定同源亚基的表达模式，我们关注了52,593个表达基因（总共18,263个三联体中17,531个有1:1:1基因对应关系，占95.99%），并且至少在一个组织中有>0.5个每百万转录本（TPM）的同源亚基表达(表S16，参见在线补充材料)。针对这52,593个表达三联体，我们首先针对每个个体组织和综合分析，标准化了BXJ1、BXJ2和BXJ3同源亚基的相对表达。然后，在所有四个组织中进行了全局分析(图4A)。平均来说，约40％的同源亚基三联体显示出不平衡的表达模式，其中一个同源亚基相对于其他两个同源亚基具有更高（亚基优势）或更低（亚基抑制）的表达量。针对每个组织，平衡的三联体比例从果实的49.45％到叶片的55.00％不等(图4B)。亚基优势的三联体在每个组织中都以适度的频率出现（从11.95％到12.83％），而亚基抑制的三联体更为频繁（从33.06％到37.72％）。显著的是，统治组中来自同源三联体的基因比平衡组或抑制组的基因具有更高的绝对转录量(图4C)。此外，我们发现非平衡的三联体显示出比平衡的三联体稍高的Ka/Ks值(图4D)。对于非平衡三联体进行了GO富集分析，但这些基因未观察到显着的功能富集。非平衡表达基因密度在图4E中呈现，显示这些基因在染色体末端区域密度更高。