毛竹泛基因组-文献精读52

Haplotype-based pangenomes reveal genetic variations and climate adaptations in moso bamboo populations

基于单倍型的泛基因组揭示了毛竹种群中的遗传变异和气候适应性

摘要

毛竹（Phyllostachys edulis）是东亚地区一种在生态和经济上都具有重要价值的森林物种，在碳汇和气候变化缓解方面发挥着关键作用。然而，日益加剧的气候变化威胁着毛竹的生存。本研究中，我们为16个具有代表性的毛竹品种构建了高质量的基于单倍型的泛基因组组装，并将这些组装与427个先前重测序的品种进行了整合。对基于单倍型的泛基因组的表征揭示了广泛的遗传变异，这些变异主要发生在单倍型之间，而非品种内部。许多具有等位基因特异性表达模式的基因与气候响应相关。结合时空气候数据，研究发现超过1050个变异与温度和降水等关键气候因素相关。与气候相关的变异使得在未来排放情景下可以预测中国北部和西部地区的遗传风险增加，凸显了气温上升带来的威胁。我们的综合单倍型泛基因组揭示了毛竹的局部气候适应机制，并为应对气候压力日益加剧的这一重要竹类提供了关键的基因组资源。更广泛地说，本研究展示了长读长测序技术在解析气候敏感物种适应性状中的强大能力，推动了支持保护工作的进化知识发展。

引言

毛竹（Phyllostachys edulis）是东亚广泛分布的一种生态和经济上极其重要的森林物种，在碳汇和气候变化缓解中发挥着关键作用1。2020年，中国竹产业的市场价值达到415.8亿美元，预计到2035年将超过1386.3亿美元2,3，这突显了竹产业在全球市场中的日益经济重要性及其潜在影响。由于其生长迅速且具有高碳储量，一公顷毛竹林在60年内的碳储量可达到杉木林的两倍4。然而，日益加剧的气候变化严重威胁着植物生物多样性和生存，这将阻碍通过碳汇减少温室气体的进程，并影响《巴黎协定》的实施5,6,7。对于像毛竹这种栖息地限制的物种来说，遗传变异的持续积累对于产生能够适应气候变化的个体至关重要8,9。此前的一项研究表明，毛竹种群经历了两次瓶颈事件，导致其遗传多样性仅为0.02710。这种多样性的缺乏，加上毛竹对气候的固有敏感性11,12，凸显了在气候变化加剧背景下其生存面临的风险。因此，对毛竹气候风险和适应能力的定量评估将有助于制定保护政策和育种工作；但这类研究仍然稀缺。

近年来，通过基因组方法评估植物种群对气候变化的脆弱性取得了重大进展13。多项关于森林树种的研究证明了整合单核苷酸多态性（SNP）和环境数据，检测其内在关联，能够有效识别气候适应变异并预测气候响应14,15。基因组偏移代表了因气候迅速变化导致当前基因型-气候关系的破坏16,17。基因组偏移方法的优势在于它能够通过基因组数据预测种群对气候变化的响应和脆弱性，是常见园林实验的替代方法18,19。正向和反向基因组偏移还可以在迁移情景下预测种群的不适应性13。然而，SNP可能无法完全捕捉与气候相关的基因组信息，因为插入缺失（InDels）和结构变异（SVs）等其他类型的遗传变异通常被忽略14。此外，由于测序分辨率的限制，系统和准确检测SVs和等位基因特异性表达（ASE）具有挑战性，尽管它们在植物的物种形成和适应性进化中发挥重要作用20。这些检测遗传变异的局限性阻碍了与气候适应相关位点的识别。因此，单倍型基因组测序能够以更高的分辨率定位遗传变异，加速这些复杂基因组特征的剖析，并改善对遗传多样性、基因型-表型关联和环境适应的理解21。

近期研究也显示，泛基因组分析通过整合多个品种的全新基因组组装，可以全面阐明遗传多样性22。然而，泛基因组分析仍依赖于短读长测序，与长读长测序技术相比，在系统捕捉SVs方面的能力有限。尽管泛基因组分析变得越来越普及，但对非模式物种23（如毛竹）中的长读长单倍型和SV的研究仍然稀少。因此，使用长读长测序构建高质量的基于单倍型的泛基因组将使得与毛竹适应性相关的遗传变异的鉴定更加全面。

在本研究中，我们采用PacBio HiFi和Hi-C测序策略，为毛竹构建了基于单倍型的泛基因组，以阐明该物种广泛分布和气候适应的遗传基础。通过整合16个具有代表性的毛竹品种（RMA）的综合基因组数据集，我们高分辨率表征了全基因组的遗传变异和ASE。此外，通过利用基于图的泛基因组和高分辨率时空气候数据，我们鉴定了与局部气候适应相关的遗传位点，并量化了中国毛竹种群的气候不适应风险。我们的研究回答了以下关键问题：1）毛竹的单倍型水平基因组多样性的范围和模式是什么？2）ASE如何促进毛竹的适应性韧性？3）哪些基因组变异支撑了毛竹种群的局部气候适应？4）预测的气候变化情景将如何影响毛竹种群的气候不适应风险？回答这些问题为制定基于证据的保护和育种策略提供了关键见解，以保护这一在快速全球气候变化面前生态和经济上至关重要的物种。

结果

16个代表性毛竹品种的基因组测序和变异识别

我们根据毛竹在中国的地理分布，选择了16个具有代表性的毛竹品种（RMA），以捕捉广泛的遗传多样性（补充图1）。所有RMA的基因组均进行了全新测序，共生成1.03 Tb的PacBio HiFi读取数据，平均每个RMA的测序深度为33.74×（补充表1）。这些数据使得我们构建了16个高质量组装（32个单倍型组装），平均contig N50长度为57.0 Mb（表1和补充图2）。最终组装的平均质量值（QV）为64.26，k-mer完整性为98.20%（补充图3）。我们观察到所有组装的平均交换错误率为5.44%（表1和补充表2）。此外，我们利用Hi-C测序对三个RMA（CY、HB和HZP）进行了染色体水平的基因组组装，平均测序深度为139.24×（补充表3）。Hi-C相互作用用于锚定、排序、定向和修正contig错误连接，使得单倍型组装中99.07%的序列锚定到伪染色体（表1）。

Accession	Assembly length (Gbp)	Contig N50 (Mbp)	Sequences anchored to chromosomes (%)	Completeness (BUSCO, %)*	Switch errors (%)	Annotated genes	LAI*	TEs (%)*
AJ	1.97/1.88**	72.2/46.9	98.2/99.5	98.7/96.4	6.43	55,190/52,271	16.0/15.2	65.4/65.3
CS	1.97/1.91	58.9/46.8	99.0/99.7	98.7/98.3	6.89	53,812/51,062	16.6/15.9	64.9/65.0
CY	1.98/1.96	69.3/66.9	98.9/99.7	98.8/98.7	2.36	55,903/55,465	16.0/17.1	64.5/64.1
DA	1.93/1.90	46.0/34.3	99.2/99.8	98.1/98.6	9.06	54,118/52,947	16.0/16.4	65.2/64.9
HB	1.98/1.97	58.4/54.9	98.9/99.5	99.1/98.9	3.58	54,656/55,112	16.3/16.7	64.7/64.9
HS	1.99/1.94	74.1/61.8	98.8/99.6	99.4/99.2	3.39	58,130/56,413	17.1/16.6	60.7/61.3
HZP	1.99/1.97	64.3/61.8	97.8/99.4	99.1/99.0	3.56	55,781/56,456	16.2/16.1	64.5/64.5
JZ	1.98/1.94	62.1/45.1	98.7/99.4	99.2/98.8	4.52	53,426/54,871	16.2/16.0	64.6/65.1
LY	2.00/1.93	59.3/41.0	98.3/99.6	99.4/99.0	4.40	53,416/51,582	17.0/16.2	64.8/64.6
RH	1.98/1.93	78.7/50.2	98.7/99.2	99.3/99.1	4.78	54,220/52,994	16.8/16.6	66.2/65.6
WYS	1.98/1.92	61.0/32.4	98.6/99.6	99.0/98.0	6.12	55,658/54,087	16.2/16.6	65.9/64.7
XA	1.97/1.94	63.3/60.3	98.7/99.6	98.8/99.2	5.73	55,207/52,029	16.8/17.0	65.4/64.5
XN	1.99/1.93	57.9/33.8	98.5/99.5	98.9/99.0	5.54	53,867/53,984	16.0/15.2	64.7/65.1
YA	1.99/1.95	73.8/59.7	98.1/99.5	98.8/99.4	6.49	55,876/54,896	16.4/16.8	64.9/65.1
YF	1.92/1.88	48.4/54.8	98.6/99.7	96.7/96.8	8.40	53,988/52,498	16.6/16.4	64.8/64.7
YX	2.00/1.90	79.0/46.8	98.4/99.6	99.3/98.6	5.82	55,533/53,514	17.0/16.0	64.3/63.4

为了鉴定高质量的蛋白质编码基因模型，我们对每个RMA的3或4个组织的转录组进行了测序，从186个样本中生成了1.34 Tb的RNA-seq数据（补充数据1）。通过组合方法（见方法部分），我们在所有32个单倍型组装中预测了平均54,343个蛋白质编码基因模型（表1），其中平均97.9%的基因被赋予了推定功能（补充数据2）。在这些基因模型中，平均存在92,506个双等位基因（每个RMA对应46,253对等位基因），而在单倍型组装中，平均有8090个基因以单等位形式存在（补充表4）。基因组组装和注释经过了广泛的质量评估，包括连续性、碱基准确性、结构准确性和相关性分析，结果表明单倍型组装和基因模型质量高，适合下游分析（补充数据3和补充图4）。我们还在16个组装中鉴定了完整的长末端重复序列（LTRs），并对LTR家族进行了聚类（补充图5-6）。

为了全面表征16个RMA中的遗传变异，我们选择了CYhap1作为参考基因组，因为其质量优于其他品种，并成功解析了其单倍型。因此，我们鉴定了2,123,198个SNPs和294,108个InDels（图1a），统称为本研究中的短变异。我们发现，超过一半（65.9%）的短变异存在于所有16个品种中或仅在一个品种中（补充图7）。此外，我们使用基于两种策略的五种SV检测工具，检测到26,987个SVs，包括12,417个插入（INSs）、14,494个缺失（DELs）和76个倒位（INVs）（图1a、补充数据4和补充图8）。同样，63.8%的SVs在所有16个品种中或仅在一个品种中共享（补充图7）。随后，我们利用这些SVs构建了一个基于图的泛基因组。此外，我们鉴定了与LTRs相交的SVs，这些SVs的形成可能与LTRs的存在相关（补充数据5）。

a 图：圆形图展示了毛竹泛基因组，其中同心环从外到内（从 i 到 v）依次显示了基因、转座元件（TEs）、单核苷酸多态性（SNPs）、小型插入和缺失（InDels）、以及结构变异（SVs）的密度。

b 图：显示了SVs（红色）和短变异（SNPs和InDels，蓝色）的数量，按品种间（深色）和单倍型间（浅色）进行分类。品种间变异在参考单倍型之间不存在，但在其他品种中存在。单倍型间变异存在于参考单倍型之间。x轴代表品种，y轴显示SVs/短变异的数量。

c 图：显示了单倍型间和品种间SVs/短变异的频率分布（x轴）与数量（y轴）。

d 图：展示了每100 kb基因组窗口中SVs和短变异的基因组分布之间的相关性。单倍型间变异（红色）在全基因组分布中的相关性更强（双侧皮尔森相关性检验，P值 < 0.001），相比之下，品种间变异（蓝色，双侧皮尔森相关性检验，P值 < 0.001）相关性较弱。

毛竹中的大多数遗传变异发生在单倍型之间，而不是品种内部

通过比较每个RMA与参考基因组单倍型之间的遗传差异，我们发现平均97.0%的杂合变异也存在于参考基因组的单倍型之间。因此，这些变异被分类为存在于每个品种的单倍型之间（在本研究中称为“单倍型间”），而不是品种之间（称为“品种间”）。单倍型间短变异和SV的数量分别平均是品种间短变异和SV的10.4倍和5.3倍（图1b和补充表5）。总短变异和品种间短变异的平均杂合度水平分别为98.6%和71.2%（补充数据6），这表明传统的变异鉴定方法高估了毛竹的杂合度。平均而言，我们每923.9 bp鉴定出一个短变异，每23,105.9 bp鉴定出一个品种间短变异（补充表6），这表明毛竹的实际遗传多样性低于之前的报告10。

在所有RMA中，对于短变异和SV，68.5%和68.8%的单倍型间变异在所有品种中检测到，频率较高，而只有0.1%和0.5%的单倍型间变异仅在一个品种中检测到，频率较低。相比之下，大多数品种间变异（46.3%的短变异和49.7%的SV）呈现低频，而只有0.2%和0.9%分别为高频（图1c）。此外，我们比较了三个相位良好的基因组单倍型，并计算了SNP的数量，结果显示来自不同品种的相同单倍型之间差异相对较小，而不同单倍型之间的变异较大（补充表7）。因此，结果表明品种中包含的单倍型间变异远超品种间变异。标准差显示，单倍型间变异的基因组分布也比品种间变异更加集中（补充表8和补充图9）。我们还观察到，与品种间变异相比，单倍型间变异的全基因组短变异和SV分布之间存在强正相关（Pearson’s r = 0.79），而品种间变异的正相关性较弱（Pearson’s r = 0.28）（图1d）。就SV长度而言，长度小于100 kb的SV中，单倍型间SV占比较多，而长度大于100 kb的SV中，95%为品种间SV（补充图10）。

通过等位基因比较构建和表征基于单倍型的毛竹泛基因组

为了更好地表征单倍型并对等位基因进行比较分析，我们基于蛋白质序列相似性和基因间距离对等位基因进行了比对和鉴定（见方法部分）。因此，通过鉴定和比较RMA中的等位基因，我们构建了基于单倍型的毛竹泛基因组。总计1738,962个基因被分为来自32个单倍型组装的74,758个基因集。基因集根据其在品种中的存在情况被分类为核心基因集（存在于所有16个品种中）、软核心基因集（存在于12–15个品种中）、可有可无基因集（存在于2–11个品种中）或私有基因集（仅存在于一个品种中）。这四类的比例分别为53.90%、16.94%、28.06%和1.10%（图2b）。此外，基因集根据等位基因组成被分为三组：双等位基因集（所有品种中检测到等位基因对）、单等位基因集（所有品种中仅检测到每对等位基因中的一个）和可变等位基因集（在部分品种中检测到等位基因对）。补充图11提供了分类示意图。在基因频率和等位基因组成类别组合的12组中，核心基因集占双等位基因集的最大比例（92.1%），而在单等位基因集中占最低比例（0.3%）（图2c）。

a 图：随着分析中纳入的毛竹品种数量（x轴）的增加，泛基因组中的基因集数量（蓝色）和核心基因集数量（红色）增加。误差棒表示平均值±标准差，n为模拟次数。

b 图：泛基因组的组成。柱状图显示每个品种中基因集的数量（y轴），按频率分类（x轴）。饼图显示每个组成类别的基因集比例：核心、软核心、可有可无和私有基因集。左侧块显示每个基因组中唯一基因集的数量（底部）和未聚类基因的总和（阴影区域）。

c 图：基因集在不同组中的分布，根据基因频率和等位基因组成进行划分。y轴表示根据品种中的基因频率划分的四组，x轴表示根据等位基因组成分类的三个基因集组。所有基因集被分为12组（核心-双等位基因（19,270），核心-可变等位基因（20,858），核心-单等位基因（169），软核心-双等位基因（561），软核心-可变等位基因（11,772），软核心-单等位基因（330），可有可无-双等位基因（267），可有可无-可变等位基因（17,010），可有可无-单等位基因（3,702），私有-双等位基因（819），私有-可变等位基因（5998），私有-单等位基因（44,104））。每组的面积与基因集数量成比例。

d、e 和 f 图：对比12个基因集组（x轴）中基因长度（d）、表达水平（e）和组织特异性指数（Tau）（f）。y轴显示基因长度（以碱基对为单位）、log2(TPM + 1)表达值和Tau特异性指数。箱线图显示中位数（中心线）、四分位范围（箱体）和1.5倍四分位范围（须）。n为基因集数量。统计显著性P值见补充表9–11。

对12个基因集组的特征分析显示了基因结构、表达模式和功能特征的差异。在相同单倍型类别中，核心基因集的基因长度、cDNA长度、编码DNA序列（CDS）数量和CDS大小均大于私有基因集，而双等位基因集的这些指标大于单等位基因集（Wilcoxon符号秩检验，P值<0.001）（图2d，补充图12-14和补充表9）。平均基因表达水平（每百万碱基转录数，TPM）从核心基因集到私有基因集逐渐下降，并从双等位基因集到单等位基因集逐渐下降，私有基因集除外（图2d和补充表10）。组织特异性指数（Tau）在核心基因集中显著低于私有基因集（Wilcoxon符号秩检验，P值<0.001），单等位基因集的Tau值大于双等位基因集，私有基因集除外（Wilcoxon符号秩检验，P值<0.05）（图2e和补充表11）。这些结果表明，核心基因集中有更多基因被表达且其表达水平更高，而私有基因集中的基因表现出更高的组织特异性。类似的模式也出现在双等位基因集和单等位基因集中。

此外，我们重点关注了核心-单等位基因集，这些基因存在于所有品种中但仅存在于一个单倍型组装中。在47个具有已知功能且在至少一个RNA-seq品种中TPM大于1的核心-单等位基因集中，鉴定出了27个与环境适应相关的基因集（补充数据7）。这些基因集的功能包括抗逆性（如基因集GS0035370，编码醛酮还原酶1（AKR1））24、抗病性（如基因集GS0058418，编码抗病蛋白RPM1）25，以及DNA损伤修复（如基因集GS0062031，编码（6-4）光解酶的动态）26。这些结果表明，一些单倍型特异性基因可能参与了毛竹的环境适应。

等位基因特异性基因表达在毛竹环境适应中的作用

等位基因特异性表达（ASE）分析揭示了16,317个基因在16对单倍型组装中表现出ASE（补充数据8）。这些基因分布在泛基因组中的8730个基因集中，其中比例最大的是核心-可变基因集（34.1%），最少的是可有可无-双等位基因集（0.3%）（图3a）。此外，等位基因特异性表达基因集（ASEG）的分布显示出高度的品种特异性，其中81.8%（7139个）ASEG仅在1-2个品种中检测到，而只有9个ASEG在所有品种中共享（图3b）。有趣的是，来自同一基因集的ASEG在不同组织间也表现出可变的等位基因特异性表达模式。例如，在不同单倍型组装中检测到的3149个ASEG中，有72%表现出在品种间变异的组织特异性表达模式，而其余28%在品种间表现出一致的组织表达（补充数据9）。

a 图：毛竹泛基因组中等位基因特异性表达基因集（ASEGs）的分布，按品种（内圈）和等位基因（外圈）分类。

b 图：ASEGs在所有品种中的频率分布（y轴）和品种（x轴）。

c 图：在毛竹泛基因组中检测到的ASEGs在四种组织（叶、茎、根茎和根）中的组织特异性分布。

d 图：使用单侧超几何分布检验对在四种组织中均检测到的ASEGs及仅在某一种组织中检测到的ASEGs进行基因本体（GO）富集分析。圆圈颜色表示统计显著性（P值），大小表示ASEGs的数量。富集因子为注释到给定GO术语的ASEGs与该术语中总基因数的比率。

e 图：一个一致的ASEG（基因集GS0010347）的ASE模式示例。该ASEG在单倍型1中表现出高表达（红色），而在单倍型2中表现出低表达（蓝色），在茎和根茎中均如此。

f 图：GS0010347中的一个6398 bp的SV的示意图，该SV可能通过改变单倍型间的蛋白质序列来引发ASE。基因区域为绿色，CDS为蓝色，5’/3’UTR为红色。

g 图：一个不一致的ASEG（基因集GS0027844）的组织特异性模式示例。在五个品种（AJ、CY、RH、XA和XN）中，ASEGs在叶中表现出单倍型1的表达高于单倍型2。相反，在根茎中，ASEGs在单倍型1中的表达低于单倍型2，相对于叶的模式。

h 图：GS0027844基因集中两个单倍型CDS之间的一个DEL和多个碱基替换的示意图。这些序列差异（红色）改变了两个单倍型之间的蛋白质序列。

ASE组织特异性分析显示，最高比例的ASEGs（39.3%）在所有组织中表达，其次是仅在叶（13.4%）、根（6.1%）、茎（5.1%）和根茎（4.9%）中表达的ASEGs（图3c）。 这些ASEGs的功能富集分析显示其与各种刺激和防御反应有关（图3d）。在所有组织中表达的ASEGs主要富集于与蛋白质生物合成和修饰相关的过程，而那些具有组织特异性表达的ASEGs则参与了发育过程，如叶中的蜡生物合成过程和茎中的细胞壁生物发生（图3d）。这些结果表明，毛竹的ASEGs可能在环境适应中起到关键作用，同时也有助于组织特异性的发育过程。

后续分析发现5156个一致的ASEGs（在所有组织中对一个等位基因表现出一致的表达模式），来自3183个基因集（补充数据10）。例如，基因集GS0010347在所有品种的茎和根茎中表现出一致的ASEGs，是单倍型组装中普遍存在的ASEG（图3e）。GS0010347的CDS与一个6398 bp的单倍型间DEL重叠（图3f），可能导致了ASE的产生。该基因集编码含黄素单加氧酶1（FMO1），实验表明其参与建立系统获得性抗性27。这些结果表明，某些单倍型间变异可能与一致的ASE事件相关。此外，我们鉴定了68个不一致的ASEGs（在不同组织中对等位基因表现出高表达切换），来自47个基因集（补充数据11）。例如，基因集GS0027844在5个品种（AJ、CY、RH、XA和XN）的叶中表现出单倍型1的ASE，但在根茎中表现出相反的模式，即单倍型2的高表达（图3g）。ASE可能是由CDS中的6 bp DEL和碱基替换引起的（图3h）。该基因集编码含有NAD依赖性表异构酶/脱水酶结构域的蛋白质。

鉴定与生物气候变量相关的基因组变异

为了识别可能与气候适应相关的位点，我们基于从图泛基因组中鉴定的SVs和427个先前重测序的毛竹品种中检测到的小变异，进行了基因型-环境关联（GEA）分析10。我们应用潜在因素混合模型2（LFMM2）和冗余分析（RDA），使用WorldClim的19个生物气候（BIO）变量，包括11个温度相关变量（BIO1–BIO11）和8个降水相关变量（BIO12–BIO19），来识别与气候相关的遗传变异（补充数据12）。首先，我们通过ADMIXTURE分析SNPs确认了种群结构，发现K = 1；这些结果与之前的研究一致10（补充图15）。LFMM2初步检测到96,638个SNPs、7456个InDels和449个SVs显著与生物气候变量相关（FDR校正P值<0.05）（补充数据13-14）。基于变量相关性和梯度森林（GF）排名（见方法部分）（补充图16-17），我们选择了六个生物气候变量用于RDA进一步筛选变异。这些变量包括年均温（BIO1）、日均温差（BIO2）、最热月最高温（BIO5）、最湿季平均温度（BIO8）、最干月降水量（BIO14）和降水季节性（BIO15）。在保留通过LFMM2和RDA两种方法识别的变异后（补充图18和补充数据14-15），我们鉴定了与生物气候变量相关的1050个适应性变异（958个SNPs、90个InDels和2个SVs）（补充数据16），代表了毛竹气候适应的核心基因组变异。此外，与123个与降水相关的变异相比，有996个变异与温度相关。RDA显示，气候效应解释了适应性变异中35%的遗传变异，这远高于地理因素的贡献（补充表12）。

为了验证气候关联性并检验这些变异的潜在适应作用，我们重点关注了与温度和降水相关的几个顶级候选变异。例如，我们在染色体19上鉴定了一个与BIO5相关的SNP，其中最显著的是chr19_24871064_SNP（LFMM2 FDR校正P值=0.019）（补充图19）。纯合形式显示出较高的BIO5值（图4b）。在较高BIO5地区（如XN、LY），G等位基因频率较高，而在较低BIO5地区（如XA、HZP），G等位基因频率较低（图4c），相关性为0.62（补充图20）。该变异位于CY_hic_hap1_01Gene018743的内含子中（图4a），而其他一些相关SNP位于该基因的下游，该基因编码与热应激相关的CCHC型锌指蛋白28。我们还检测到染色体18上的一个与BIO5相关的SV（chr18_28562210_SV）（LFMM2 FDR校正P值<0.001）。拥有此SV的品种表现出较高的BIO5值（补充图21）。距离该SV最近的基因是CY_hic_hap1_01Gene010919，编码脱落酸受体PYL8，在调节应激反应中发挥重要作用29。

a 图：与最热月最高温（BIO5）相关的位点位于染色体19上的CY_hic_hap1_01Gene018743的内含子和下游。 d 图：与最干月降水量（BIO14）相关的位点位于染色体13上的CY_hic_hap1_01Gene048977的下游。 b, e 图：不同毛竹品种基因型间BIO5（b）和BIO14（e）值的比较。箱线图显示中位数（中心线）、四分位范围（箱体）和1.5倍四分位范围（须）。BIO5（b）的样本数n = 268（0/1）和159（1/1），BIO14（e）的样本数n = 375（0/1）和50（1/1）。统计显著性通过双侧Wilcoxon秩和检验确定，P值分别为9.29e-10 < 0.001（b）和1.99e-27 < 0.001（e）。 c, f 图：与BIO5（c：chr19_24871064_SNP）和BIO14（f：chr13_64739621_SNP）相关的候选适应性SNP的等位基因频率。圆圈表示等位基因频率，红色表示替代等位基因，蓝色表示参考等位基因（c），黄色表示替代等位基因，蓝色表示参考等位基因（f）。地图颜色反映了毛竹种群分布范围内BIO5（c）和BIO14（f）值的变化。

对于与降水相关的BIO14，我们在染色体13上鉴定了一个相关的变异（补充图22）。最显著相关的变异是chr13_64739621_SNP（LFMM2 FDR校正P值<0.001）。携带此变异的品种表现出较低的BIO14值（图4e）。在BIO14较低的HZP和CS地区，G等位基因频率较高（图4f），且存在强正相关（Pearson’s r = 0.81）（补充图23）。这些变异位于CY_hic_hap1_01Gene048977的下游（图4d），该基因邻近一个编码与干旱应激相关的受体样蛋白激酶（RLK）的基因30。

预测脆弱的毛竹种群

利用识别出的与气候相关的变异和未来气候预测，我们计算了不适应风险（RONA），该值代表毛竹为适应未来条件所需的预期等位基因频率偏移。考虑了四个通用环流模型（GCMs）：澳大利亚社区气候与地球系统模拟耦合模型版本2（ACCESS-CM2）31、第二代CMCC地球系统模型（CMCC-ESM2）32、戈达德空间研究所E版本2.1耦合GISS海洋模型（GISS-E2-1-G）33和气候跨学科研究模型版本6（MIROC6）34，它们参与了世界气候研究计划耦合模型比较计划第6阶段（WCRP CMIP6），基于两个共享的社会经济路径（SSPs）。两个SSPs包括低温室气体排放情景（SSP126）和高温室气体排放情景（SSP585）。结果显示在高排放情景SSP585下的RONA值高于低排放情景SSP126下的RONA值。对于所有温度相关变量，RONA值均高于降水相关变量，并且SSP585与SSP126之间的差异更大（图5a）。对于温度，整体趋势显示未来将上升；因此，我们特别关注BIO5。在BIO5预测下，最高的RONA出现在XN，这可能是由于其目前暴露于最高温度和该地区预计的最大升温（图5b）。LY、YF和AJ地区也显示出较高的BIO5 RONA值，表明这些种群可能需要采取针对极端高温事件的保护措施。

a 图：比较2061年至2080年期间在19个生物气候变量（BIO1–BIO19，x轴）和四个气候模型（ACCESS-CM2、CMCC-ESM2、GISS-E2-1-G和MIROC6）下，不同排放情景SSP126（红色）和SSP585（蓝色）的平均不适应风险（RONA，y轴）。误差棒表示平均值加减标准误差。

b 图：在高排放情景（SSP585）和最热月最高温（BIO5）下，从2061年到2080年基于四个气候模型（ACCESS-CM2、CMCC-ESM2、GISS-E2-1-G和MIROC6）的毛竹种群的平均RONA估计值。地图颜色表示BIO5的预测气候变化，红色越深表示温度增加越显著。圆圈大小表示不同种群的RONA值。

c, d 图：梯度森林（GF）预测的2061年至2080年期间在四个气候模型下毛竹分布区域的局部基因组偏移平均值图，分别为SSP126（c）和SSP585（d）。颜色从蓝到红表示基因组偏移的增加。

e, f 图：显示SSP126（e）和SSP585（f）下局部（红色）、正向（绿色）和反向（蓝色）基因组偏移的RGB地图。颜色越亮（接近白色），表示沿各轴基因组偏移值相对较大；颜色越暗（接近黑色），表示值相对较低。下方面板是（e, f）的双变量散点图，包含1:1线。

此外，结合所有生物气候变量的GF模型预测了各区域的局部基因组偏移，代表了气候变化下的脆弱性。预测显示风险最大的是西北地区（图5c, d）。除了局部基因组偏移外，我们还计算了正向和反向基因组偏移（图5e, f 和补充图24, 25）。图5e, f中将三种基因组偏移（局部偏移、正向偏移和反向偏移）以红、绿、蓝色带在RGB色彩空间中进行组合可视化。较亮的单元格（接近白色）和较暗的单元格（接近黑色）分别表示各轴上较大的和较小的值。来自北部地区的大多数品种表现较亮（图5e, f），表明其偏移值相对较高，这暗示即使有迁移，这些品种仍然比南部地区面临更大的脆弱性。然而，在大多数北部地区，正向和反向偏移均低于局部偏移（图5e, f下方面板），这表明辅助迁移在某种程度上可以促进适应未来气候变化。与RONA结果一致，SSP585情景下的所有基因组偏移均高于SSP126情景，高基因组偏移的区域（较亮区域）也更大，这表明与化石燃料开发相关的极端气候变化（SSP585）可能比更可持续的情景（SSP126）使毛竹种群面临更大的适应挑战和潜在风险。毛竹在主要自然分布区在SSP126情景下仍处于相对安全的位置。然而，在SSP585情景下，一些主要的自然生长区域，尤其是最西部的两个自然生长区域，将面临风险（图5f）。

讨论

单倍型解析基因组、基于图的泛基因组和属级泛基因组的引入和应用极大地丰富了我们对物种或分类群基因组多样性的理解，为揭示重要性状的遗传基础提供了强大的工具23,35,36。在此背景下，我们对毛竹（Phyllostachys edulis）这一具有重要经济和生态价值的非木材资源进行了研究，通过采用第三代PacBio HiFi和Hi-C测序技术获得了16个代表性毛竹品种的单倍型基因组，并构建了综合泛基因组。这些基因组资源不仅更全面地捕捉了毛竹基因组的异质性，还为深入了解毛竹对多样环境条件的适应性提供了宝贵的遗传信息。然而，尽管单倍型基因组和泛基因组相较于传统的折叠基因组有显著优势，在实际应用中仍存在挑战。主流的组学分析工作流程，如转录组学和表观基因组学，仍主要依赖于将测序数据比对到线性参考基因组，而未能充分利用单倍型基因组和泛基因组中蕴含的丰富多样性信息。因此，在不断提高基因组数据的准确性和完整性的同时，必须加强这些高质量基因组资源的分析框架和应用环境，以确保其最佳利用。

为了充分利用这些基因组资源，我们整合了单倍型基因组和泛基因组，揭示了关于毛竹基因组结构的宝贵见解。我们的研究发现，毛竹的两个单倍型之间的差异大于两个毛竹品种之间的差异。鉴于毛竹长时间的无性繁殖及其67年的开花周期37,38，变异的主要来源可能是一个单倍型内罕见的体细胞突变。无性繁殖使品种中积累的变异难以传递，因为缺乏减数分裂阻止了同源染色体之间的遗传物质交换（补充图26）。我们假设不同区域的毛竹种群的共同祖先的两个单倍型之间已经存在差异，并且来自不同区域的毛竹中体细胞突变的积累未超过两个祖先单倍型之间的原始差异。这些因素导致了单倍型间变异在数量上超过不同品种之间的遗传变异。此外，我们发现传统变异检测方法可能高估了杂合度。考虑单倍型基因组时，我们发现普遍的杂合位点在参考基因组中也是杂合的，不应被视为品种间的变异位点。过滤掉这些单倍型间变异会导致检测到的杂合度下降，并表明遗传多样性低于原始估计。

由于毛竹的遗传多样性较低，对其适应多样环境条件的深入研究显得尤为重要。我们观察到核心-单等位基因集和等位基因特异性表达（ASE）现象与环境适应密切相关，并鉴定出两组与气候相关的杂合变异位点，这可能意味着单倍型在毛竹环境适应中起着重要作用。我们的研究还显示在高排放情景SSP585下，毛竹种群面临显著的适应性风险（图5），强调减排措施在缓解气候变化压力方面的重要性。尤其在西北地区，我们建议优先进行栖息地恢复（图5c），并在解决潜在竞争风险的同时，考虑对北部种群的辅助迁移（图5e–f）。值得注意的是，我们的样本只包含了中国毛竹的主要自然分布区域，而未包含一些人类迁移种群或极端种群。补充这些种群，甚至全球毛竹品种，可能有助于识别更多适应极端环境的变异。对于不涉及迁移的风险预测如RONA和局部偏移，这些样本的缺失影响较小。但对于正向和反向偏移分析，加入更多的种群可能揭示更适合毛竹种植的地区，并识别更适合迁移到极端区域的毛竹种群。然而，基因组偏移在保护规划中的应用仍处于初期阶段，其预测的实用性需通过精心设计的实验进行实证验证，如共园实验或受控环境测试，将基因组偏移预测与种群在环境变化中实现的适应性结果进行比较14,19。

方法

样本采集

为了优化对遗传和环境多样性的代表性，我们根据先前的一项系统发育研究选取了16个代表性毛竹品种（RMA），该研究确定了该物种在中国的主要自然分布区10,39（补充图1）。我们的采样策略旨在通过覆盖毛竹的所有主要栖息地，捕捉其广泛的遗传多样性，确保对种群的全面代表性。16个RMA各自来自16个不同的地区（补充表1和补充图1），共同捕捉了毛竹种群中的广泛变异。此外，我们还使用了先前研究中获得的427个重测序品种的遗传信息10（补充数据17），这些品种覆盖了中国毛竹的所有主要自然分布区，从而增强了本研究中样本的遗传代表性。在每个地区，毛竹笋于2020年4月采集用于DNA提取。同时，在每个地区从同一毛竹根茎上采集年轻的叶片、茎、根和根茎作为RNA-seq样本。DNA样本使用硅胶珠快速干燥，RNA样本置于RNA稳定液中。所有样本在干冰上保存和运输。代表性标本随后存放于国际竹藤中心。

基因组和转录组测序

高分子量基因组DNA（gDNA）提取过程中严格控制质量。使用Nanodrop 1000分光光度计和Qubit检测仪（Thermo Fisher Scientific, CA, USA）评估纯度和数量。使用SMRTbell Express Template Prep Kit 2.0（Pacific Biosciences, CA, USA）构建约15 kb插入的SMRTbell文库，并使用SageELF系统（Sage Science, Beverly, MA, USA）将其按窄范围大小选择，以提高测序准确性。在Sequel II平台（Pacific Biosciences, CA, USA）上使用2-3个SMRT Cells 8 M进行测序，每次运行30小时以最大化数据产量和质量。

为了进行染色体水平的组装，按照标准协议构建Hi-C文库，并在DNBSEQ-T7平台上采用150 bp双端策略进行测序。对于RNA-seq，每个RMA的三个生物重复的组织（包括叶、茎、根和根茎）被收集（补充数据1）。使用CTAB方法提取总RNA，随后使用NEBNext Ultra RNA Library Prep Kit构建文库，并在DNBSEQ-T7平台上进行150 bp双端测序。本项目的所有文库制备和测序均由安诺优达基因科技有限公司完成。

基因组组装

使用CCS算法v6.2.040从子读段生成HiFi读取，并使用Hifiasm v0.16.1-r37541将其组装成contigs。对于CY、HB和HZP，使用Hi-C测序读取作为Hifiasm的输入，以协助生成单倍型组装。组装质量通过BUSCO v5.4.342评估，使用Embryophyta_odb10数据集（1614个植物直系同源基因）来评估结构准确性。使用LTR_retriever v2.9.043计算LTR组装指数（LAI），突变率为8.51 × 10−8 38，毛竹基因组估算大小为1.92 Gb10。使用Meryl v1.4144为原始测序读取和组装生成Meryl数据库。然后使用Merqury v1.344通过比较组装和原始读取的k-mer谱计算质量值（QV）和k-mer完整性。Hi-C读取通过3D-DNA pipeline v17012345将scaffolds锚定到染色体。最终的接触图使用Juicerbox v1.5.246可视化。

基因组注释

转座元件（TE）注释结合了在DNA和蛋白质水平上的从头和同源性方法。使用RepeatModeler v2.0.347构建从头重复库，并使用RepeatMasker v4.1.2-p148将TEs分类，映射到该库和Repbase v21.1249。基因预测结合了同源性、基于RNA-seq和从头方法。七种植物物种（凤梨、拟南芥、大豆、大麦、烟草、水稻和玉米）的同源蛋白质从Phytozome v1150下载，并使用GeneWise v2.4.151比对到基因组。RNA-seq读取经Trimmomatic v0.3852预处理，并使用HISAT2 v2.1.053比对到基因组。AUGUSTUS v3.4.054在转录组训练下进行从头预测。GETA v2.5.5（GitHub - chenlianfu/geta）用于整合所有方法的证据，为每个品种生成高质量的基因集。使用BUSCO v5.4.3评估基因集的完整性。基因功能注释根据与以下数据库的最佳匹配使用DIAMOND55进行分配：NCBI Nr v2021082456、Swiss-Prot v2021082457、KOG v2003072158、eggNOG v5.059、InterPro v8560、Pfam-A v34.061、GO（与eggNOG和InterPro整合）和KEGG数据库62（基于KAAS v2.163）。所有数据库访问日期为2021年12月11日。

切换错误计算

切换错误通过calc_switchErr v1.021计算。简而言之，使用pbmm2 v1.5.0（GitHub - PacificBiosciences/pbmm2: A minimap2 frontend for PacBio native data formats）将PacBio HiFi读取对齐到单倍型1，并使用DeepVariant v1.4.064识别SNP。随后在Whatshap v1.165上进行SNP分相，然后使用MUMmer4 v4.0.0rc166中的nucmer对两个单倍型基因组进行比对。使用show-snps命令识别SNP，并通过将这些SNP与对齐读取中识别的SNP进行比较来计算切换错误。对于有Hi-C数据的品种（CY、HB和HZP），使用bwa v0.7.1767将Hi-C测序读取对齐到单倍型1，使用DeepVariant调用SNP，并从HiFi和Hi-C SNP中计算切换错误。

等位基因识别

在任何两个单倍型组装之间确定等位基因是基于蛋白质序列相似性和相对位置的组合。蛋白质序列相似性通过使用BLASTP v2.9.0 + 68对每个单倍型基因组进行比对来计算。通过使用Minimap2 v2.24-r112269对组装进行比对来确定来自不同基因组的两个基因之间的相对位置，并保留40 kb范围内的基因对（毛竹基因组中的平均基因距离）进行进一步分析。选择40 kb作为阈值是因为它近似于毛竹基因组的平均基因距离（约2 Gb，含约50,000个基因）。该阈值允许一定的定位灵活性，使等位基因可以相隔一个平均基因间隔。根据更接近的基因组距离和更高的序列一致性，保留最佳互惠比对作为等位基因对。满足其他单倍型的相似性阈值但不在最佳比对等位基因对中的基因被识别为单倍型特异性重复基因（补充数据18）。为提高透明度和可重复性，该步骤的脚本已公开在GitHub（GitHub - ZhaoGroupLab/moso-bamboo-pangenome）。

泛基因组构建

通过迭代将每个基因组的基因整合到聚合基因集中构建泛基因组。任何基因组都可以作为初始参考。对于每个附加基因组，其基因根据识别的等位基因映射到现有泛基因组中的等位基因集。没有映射等位基因的基因被放入新集。此过程反复进行，直至所有基因组被整合。通过基因组添加顺序排列生成泛基因组聚合曲线。

泛基因组中的基因集定义

泛基因组基因集根据其在品种中的频率进行分类。核心基因集包含所有16个品种中的基因。软核心基因集包括12至15个品种中的基因。可有可无基因集包含2至11个品种中的基因。私有基因集包括仅存在于一个品种中的基因。此外，为表征等位基因在单倍型中的存在与否，我们将所有基因集按等位基因分为三类。如果基因集中每个基因都有等位基因存在，我们将其定义为双等位基因集。相反，如果基因集中所有基因都没有等位基因存在，我们将其定义为单等位基因集。第三类介于双等位基因集和单等位基因集之间。如果基因集中有些基因有等位基因存在，而其他则没有，我们将其定义为可变等位基因集。此外，包含重复基因的单等位基因集被重新分类为可变等位基因集。因此，考虑到品种出现和等位基因存在，我们最终将所有基因集分为12类：核心-双、核心-可变、核心-单、软核心-双、软核心-可变、软核心-单、可有可无-双、可有可无-可变、可有可无-单、私有-双、私有-可变和私有-单。

等位基因特异性表达

对16个代表性毛竹品种采集的3或4个组织（叶、茎、根茎和根）的三重RNA-seq数据进行了等位基因特异性表达（ASE）分析。使用Trimmomatic v0.39预处理读取，并使用HISAT2 v2.1.0比对到相应的单倍型组装。使用StringTie v1.3.570计算每百万碱基转录数（TPM）值，并使用DESeq2 v1.34.071比较等位基因。P值<0.05且绝对log2（倍数变化）>1的等位基因被分类为具有ASE。

SNP和InDel识别

我们基于全基因组比对和读取比对策略识别SNP和InDels。对于全基因组比对策略，使用MUMmer v4.0.0rc1中的nucmer将单倍型组装对齐到CY单倍型1参考。使用show-snps命令识别SNP，并使用svim-asm v1.0.272识别InDels。对于读取比对策略，使用pbmm2 v1.5.0（GitHub - PacificBiosciences/pbmm2: A minimap2 frontend for PacBio native data formats）将PacBio HiFi读取对齐到单倍型1，并使用DeepVariant v1.4.063识别SNP和InDels。随后保留两种策略检测到的SNP和InDels。根据查询-参考比较将变异转换为参考基因组坐标。使用BCFtools v1.973的合并命令整合跨品种的SNP和InDels。

结构变异识别

基于读取比对或全基因组比对到CY单倍型1参考的五个流程识别全基因组结构变异（SVs）：（i）Minimap2 v2.24-r1122 + cuteSV v1.0.1374；（ii）Ngmlr v0.2.775 + sniffles v1.0.1275；（iii）pbmm2 v1.5.0 + pbsv v2.8.0（GitHub - PacificBiosciences/pbsv: pbsv - PacBio structural variant (SV) calling and analysis tools）；（iv）Minimap2 + svim-asm v1.0.2；和（v）Nucmer v4.0.0rc166 + Assemblytics v1.2.176。

SVs的过滤和合并

基于上述五个流程在和跨品种间过滤和合并SVs。初始条件下，SV必须至少由两个调用器支持才能保留。通过至少一种全基因组比对方法检测到的SVs>50 kb。过滤后，SVs按以下方式合并：长度较大的DEL和INV（在大于5 bp和小于100 bp的范围内）合并。彼此相距10 bp以内且BLAST比对身份和覆盖率>90%的INS合并。

基因组区域中变异的频率识别

通过使用BEDtools v2-2.25.077中的makewindows工具将基因组分为100 kb的基因组窗口，识别SVs和短变异（SNP和InDels）的频率。使用coverage命令计算每个窗口中的突变数量，并按突变计数降序排列窗口。

基于图的泛基因组构建

使用vg v1.3878按以下步骤构建基于图的泛基因组。使用vg construct命令构建初始图。使用gbwt命令生成.gbwt文件。使用snarls命令生成.snarls文件。index命令构建距离索引文件并最小化索引。此外，使用vg giraffe79将427个重测序品种的短读取比对到基于图的泛基因组。

基于重测序读取的SNP和InDel调用

原始测序读取按我们前期研究中的相同流程进行处理以确保一致性。简而言之，使用BWA v0.7.17将过滤后的重测序读取比对到CY单倍型1参考基因组。对齐的读取（BAM文件）使用SAMtools v1.980排序，并使用GATK v4.2.081去除重复项。使用GATK的联合调用方法进行SNP和InDel调用。我们根据读取为每个品种获得基因组变异调用格式（GVCF），并基于质量直接过滤SNP，删除质量分数低于50的变异。

为确定修剪的SNP集，我们使用PLINK v1.982。随后使用ADMIXTURE v1.3.083进行多个随机种子的重复种群结构评估。种群结构分析显示K = 1（补充图15），与我们之前的发现一致10。

气候相关变异的识别

从基于图的泛基因组识别的SVs由16个代表性毛竹品种构建，先前生成的427个重测序毛竹品种中检测到的SNPs和InDels用于气候相关变异的识别。我们保留了小于0.05的次要等位基因频率、大于0.2的最大缺失数据比例、以及跨所有变异类型（包括SNPs、InDels和SVs）的最小深度为3的变异，使用VCFtools v0.1.1384。由于毛竹种群中的大多数变异发生在单倍型之间而不是品种内部，我们过滤掉次要基因型频率小于0.05的位点，留下1,467,461个SNPs、103,955个InDels和4,643个SVs用于全基因组气候相关变异的识别。我们使用在R包LEA v3.6.085中实现的潜在因子混合模型2（LFMM2）测试与19个生物气候变量（BIOs）的相关性。根据ADMIXTURE v1.3.0的结果和先前研究的结果10，我们将K设置为1进行LFMM2分析，并使用假发现率（FDR）校正进行多重检验以获得校正后的P值。我们保留校正后的P值<0.05的遗传变异作为气候相关变异。此外，我们在R包gradientForest v0.1.3486中执行梯度森林（GF）分析以排名19个BIO的重要性并检查它们之间的相关性。选择6个相关系数|r | <0.6的BIO用于使用R包vegan v2.6.487进行冗余分析（RDA）。显著变异定义为超过均值3.5标准差的RDA轴负载。通过保留每个气候变量中LFMM2和RDA检测到的变异识别气候相关适应变异。此外，我们使用RDA和部分RDA（pRDA）量化地理和环境的相对贡献。经度和纬度值表征地理的解释变量，以上RDA中使用的6个BIO表征气候的解释变量。

RONA和基因组偏移的计算

未来气候预测通过世界气候研究计划耦合模型比较计划第6阶段（WCRP CMIP6）中的四个通用环流模型（GCMs）获得：澳大利亚社区气候与地球系统模拟耦合模型版本2（ACCESS-CM2）、第二代CMCC地球系统模型版本2（CMCC-ESM2）、戈达德空间研究所E模型版本2.1耦合GISS海洋（GISS-E2-1-G）、气候跨学科研究模型版本6（MIROC6）。考虑了两个共享的社会经济路径（SSPs）——低排放情景（SSP126）和高排放情景（SSP585）——针对两个20年时期（2061-2080年为SSP126和2081-2100年为SSP585）。SSPs表示共享社会经济路径和代表浓度路径（RCPs）的组合。SSP126是SSP1-RCP2.6情景的缩写。SSP1（可持续性，走绿色道路）假设逐渐转向更可持续的世界，注重人类福祉和减少不平等。RCP2.6代表一个非常低的强迫水平的减排情景88,89。类似地，SSP585是SSP5-RCP8.5情景的缩写。SSP5（化石燃料发展，走高速公路）假设快速的经济增长由化石燃料驱动，能源需求高且减缓温室气体排放的努力有限。RCP8.5代表一个非常高的基线排放情景88,89。

使用pyRona v0.3690在SSP126和SSP585情景下计算所有采样区域的不适应风险（RONA）。随后在R中执行梯度森林（GF）分析以预测不同未来情景下的基因组偏移，使用19个BIOs。随后根据方法学13计算局部、正向和反向偏移，该步骤的脚本已公开在GitHub（GitHub - ZhaoGroupLab/moso-bamboo-pangenome）。局部偏移是测量居住种群对气候变化的脆弱性的指标，通过估算气候适应性位点的等位基因频率预测变化来适应当地气候随时间的变化。相反，正向偏移假设种群具有无限的迁移能力，通过识别最小预测偏移，如果通过基因流动，繁殖体或等位基因可以迁移到范围内的任何适宜栖息地来计算。反向偏移代表当前范围内的任何种群在未来是否可能预适应特定位置。反向偏移通过识别未来气候中假设种群和当前气候中的种群之间的最小偏移来计算。对于正向和反向偏移，我们将迁移距离设置为无限，因为毛竹长时间通过无性繁殖并主要通过移植引入。