Hortic Res│南京林大:整合GWAS、TWAS与WGCNA揭示银杏次生细胞壁与木材性状发育的遗传调控机制

导读

了解银杏木材性状和次生细胞壁（SCW）发育的遗传和调控机制对于改良木材品质至关重要。本研究整合全基因组关联分析（GWAS）、全转录组关联分析（TWAS）和加权基因共表达网络分析（WGCNA），鉴定了与木材性状和SCW发育相关的关键基因。SCW中的纤维素生物合成由GbCesA4、GbCesA7和GbCesA8/8B编码的CesA4-CesA7-CesA8复合物催化。这些CesA基因与TUBA/TUBB和EG形成共表达网络，表明纤维素合成、细胞骨架引导和细胞壁重塑之间存在协调作用。此外，GbCesA8B的功能缺失仅导致纤维素含量轻微降低，支持GbCesA8A和GbCesA8B之间可能存在功能冗余。对于半纤维素生物合成，GbCSLA9A/9B和IRX9/IRX14分别是甘露聚糖/葡甘露聚糖和木聚糖合成的主要贡献者，并与UXS、UXE、IRX7、GXMT和URG7形成共表达网络，涵盖核苷酸糖供应、运输以及聚合物延伸和修饰等环节。此外，MYB46可能通过激活CSLA9转录来调控SCW中的甘露聚糖/葡甘露聚糖生物合成。研究者对于木质素生物合成，TWAS鉴定出多个参与苯丙氨酸生物合成、苯丙烷代谢和木质素单体聚合的基因，包括ADT/PDT、PAL和PER，以及可能正向调控木质素积累的MYB91和几个bHLH基因。此外，还鉴定出几个可能参与SCW发育的转录因子，包括作为假定正调控因子的GATA9，作为潜在负调控因子的WRKY12和HB15，以及可能促进管胞扩张的ELF6。本研究发现为银杏木材性状和SCW发育的遗传调控提供了有价值的见解。

1 银杏基因组中SNP和InDel变异模式

本研究对290个银杏样本进行了平均深度为7.6×的全基因组重测序，具有高测序质量（Q20>98%，Q30>94%）和正常的GC含量（34%-38%）（表S3）。经过质控后，使用BWA软件将干净读数比对到参考基因组，平均比对率为99.34%（表S3）。随后，使用GATK4.5进行变异检测。经过严格过滤，共鉴定出76,128,881个高质量SNP和1,974,424个InDel（图1A和B）。

SNP和InDel在银杏基因组中的分布相对均匀（图1D和E）。此外，几条染色体显示出多个SNP或InDel密度较高的区域（图1D和E）。突变类型分析显示，转换突变（T:A>C:G和C:G>T:A）的频率高于颠换突变（T:A>G:C,C:G>G:C,T:A>A:T,和C:G>A:T），转换/颠换比（TS/TV）约为1.59（图1C）。值得注意的是，C:G>T:A转换是最普遍的突变类型，远远超过其他类别（图1C）。

注释结果显示，大多数SNP（94%）和InDel（88.96%）位于基因间区（图1A和B）。在外显子区域鉴定出的89139个SNP中，55032个（61.74%）是非同义SNP，可能通过改变氨基酸影响蛋白质功能（图1A）。此外，在3166个外显子InDel中，2025个（63.96%）是移码InDel，而921个（29.09%）是非移码InDel（图1B）。此外，还鉴定出1658个SNP和218个InDel导致翻译起始/终止位点的获得或丢失，从而导致蛋白质翻译失败、提前终止或延伸延长（图1A和B）。这类变异尤其有害，因为它们通常导致蛋白质功能的完全丧失。

图1.银杏的遗传变异模式。（A，B）银杏基因组中SNP和InDels的注释结果。（C）银杏基因组中SNP突变类型的分布。（D，E）银杏基因组中染色体上SNP和Indel的密度。

2 12个木材性状的全基因组关联分析

高密度银杏基因组变异图谱的建立为开展GWAS提供了丰富而精确的遗传标记。为了精准定位与银杏木材性状变异相关的关键基因，研究者使用GEMMA进行GWAS，并进一步参考EMMAX获得的GWAS结果，以提高关联位点的可靠性。结果显示两个软件包的结果高度一致（表S4，图S1-S6）。除了EMMAX对某些性状（如管胞长宽比（TLW）、管胞宽度腔径比（TWL）和管胞壁厚（TWT））检测到的InDel关联结果异常增加外，绝大多数位点都被两个软件包共同鉴定到。由于GEMMA的关联结果比EMMAX更准确可靠，研究者在后续分析中主要关注GEMMA的结果。

研究者使用GEMMA，共鉴定出1650个SNP和196个InDel与银杏木材性状相关（表S4.1）。在所研究的各种性状中，纤维素含量（CC）显示的相关信号最多，有1343个SNP（81.39%）和56个InDel（28.57%），其次是TWL和纤维素结晶度（CCD），分别有98个（73个SNP和25个InDel）和94个（68个SNP和26个InDel）变异位点（表S4.1）。此外，4个多效SNP（chr10:362509447,chr10:365643063,chr10:365709595,chr10:365902016）和3个多效InDel（chr10:367769028,chr10:36855634,chr10:370533958）显示出与木材基本密度（WBD）和木材气干密度（WOD）的双重性状关联（表S4.1）。类似地，2个多效SNP（chr2:532772394和chr12:279586611）被鉴定与TWL和TWT相关（表S4.1）。

EMMAX的结果也显示，纤维素含量是关联信号最多的性状，共检测到1282个SNP和54个InDel，所有这些都与GEMMA鉴定的纤维素含量关联重叠（表S4，图S5和S6）。基于连锁不平衡原则，研究者通过将关联位点上下游各延伸100kb来划分关联区域，鉴定出180个连锁基因（表S5）。其中，120个基因在外显子区域含有非同义SNP或移码/非移码InDel，而167个基因在上游区域含有SNP或InDel，总共有172个基因具有此类变异（表S5）。

利用已发表的RNA-seq数据，研究者进一步分析了这172个基因的表达模式，以缩小候选基因范围。排除低表达基因（平均FPKM<0.5）后，保留116个候选基因进行进一步研究（表S5）。对于木材气干密度，鉴定出15个候选基因（表S5）。其中，evm.TU.chr3.1073(TAA1)编码L-色氨酸-丙酮酸氨基转移酶，催化氨基从L-色氨酸转移到丙酮酸。此外，evm.TU.chr10.1203(XYL4)，编码木聚糖1,4-β-木糖苷酶，和evm.TU.chr7.1369(MAN)，编码甘露聚糖内切-1,4-β-甘露糖苷酶，被鉴定为与WBD相关的潜在基因（表S5）。

在与管胞形态性状（包括管胞长度（TL）、管胞宽度（TW）、TLW、TWT和TWL）相关的48个候选基因中，鉴定出几个具有调控潜力的功能基因或转录因子，包括evm.TU.chr6.1080(EIN3)和evm.TU.chr6.1069(ABI5)。此外，还鉴定出evm.TU.chr1.2638(CCR)，编码肉桂酰辅酶A还原酶，这是SCW形成中木质素前体转化的关键酶（表S5）。这些基因作为管胞形态发生和调控的潜在候选者。

3 纤维素合成关键基因的鉴定

纤维素是植物细胞壁中最主要的结构多糖，为植物细胞提供结构支撑和保护至关重要。如上所述，纤维素含量是关联信号最多的性状（表S4.1）。值得注意的是，在1号染色体（chr1）的997,196,864至998,048,937bp范围内观察到关联信号的显著富集，表明该区域可能包含银杏纤维素生物合成的关键遗传位点。在该区域内，一个潜在的候选基因evm.TU.chr1.2825被精确定位并注释为纤维素合酶（CesA）基因（表S5）。系统发育分析显示，evm.TU.chr1.2825与AtCesA8具有高度同源性，暗示这两个基因之间功能相似或进化保守（图2E）。

突变注释显示，与纤维素含量相关的峰值InDel（chr1:997452697）是一个移码插入突变（图2B和C）。该突变在evm.TU.chr1.2825的第一个外显子的第9个碱基后插入“CA”，并导致在插入后第5个氨基酸处提前终止翻译，从而导致基因功能完全丧失（图2D）。此外，在该基因中观察到两个非同义SNP，分别导致第86和447位的氨基酸替换（图2C）。根据PROVEAN预测，这两个突变分别被评估为中性（0.713）和有害（-2.583）（图2D）。因此，研究者推测第86位氨基酸的突变不会显著影响蛋白质功能，而第447位的突变可能产生实质性影响。基于InDelchr1:997452697的基因分型显示，携带该InDel的样本纤维素含量显著低于不携带的样本（图2F）。

图2.evm.TU.chr1.2825的功能表征。(A)基于SNP和纤维素含量的GWAS结果曼哈顿图。(B)基于InDel和纤维素含量的GWAS结果曼哈顿图。(C)evm.TU.chr1.2825的基因结构（蓝色矩形：外显子，黑色线条：内含子）。(D)evm.TU.chr1.2825内突变位点对蛋白质功能的影响（PROVEAN评分低于-2.5为有害突变，高于为中性突变）。(E)银杏和拟南芥中CesA家族基因的系统发育分析。(F)基于evm.TU.chr1.2825基因型的纤维素含量小提琴图。统计分析采用Kruskal−Wallis检验进行;*、**、***和****分别代表P<0.05、P<0.01、P<0.001和P<0.0001。

4 半纤维素合成的关键基因鉴定

半纤维素是植物细胞壁的重要组分，通过与纤维素形成氢键和共价键来增强强度和抗拉性。在半纤维素含量（HCC）的GWAS中，研究者共鉴定出15个变异位点，并筛选出5个候选基因（表S4和S5）。其中一个候选基因evm.TU.chr9.476位于9号染色体（chr9），被注释为编码β-甘露聚糖合酶（CSLA），该酶负责以GDP-Man为底物延伸β-1,4-甘露聚糖链，产生半纤维素中的甘露聚糖骨架（表S5）。系统发育分析表明，该基因与拟南芥中的AtCSLA9具有同源性（图3E）。

研究者进一步分析显示，在evm.TU.chr9.476的外显子中存在一个非移码缺失（chr9:94004237）和两个非同义SNP（chr9:93999325和chr9:93999567）（图3C）。具体来说，非移码缺失导致第41位的亮氨酸丢失，而两个非同义SNP分别导致第582位的色氨酸（W）被精氨酸（R）替换，以及第537位的丙氨酸（A）被苏氨酸（T）替换（图3C）。研究者使用InterPro进行功能结构域分析显示，第41位的亮氨酸位于TMHMM预测的跨膜螺旋区（TMhelix），而第582位的色氨酸和第537位的丙氨酸位于Phobius预测的细胞质结构域。研究者使用PROVEAN对变异进行进一步评估显示，只有第41位亮氨酸的缺失被预测为有害的（-9.189）（图3D）。

因此，研究者基于InDelchr9:94004237进行了基因分型。结果显示，具有纯合TCAG基因型（保留第41位亮氨酸）的个体半纤维素含量显著高于具有杂合缺失基因型（TCAG/T）和纯合缺失基因型（T）的个体（图3F）。此外，具有T基因型的个体半纤维素含量低于具有TCAG/T基因型的个体（图3F）。

图3.evm.TU.chr9.476的功能表征。(A)基于SNP和半纤维素含量的GWAS结果曼哈顿图。(B)基于InDel和半纤维素含量的GWAS结果曼哈顿图。(C)evm.TU.chr9.476的基因结构（蓝色矩形：外显子，黑色线条：内含子）。(D)evm.TU.chr9.476内突变位点对蛋白质功能的影响（PROVEAN评分低于-2.5为有害突变，高于为中性突变）。(E)银杏和拟南芥中CSLA家族基因的系统发育分析。(F)基于evm.TU.chr9.476基因型的半纤维素含量小提琴图。统计分析采用了Kruskal-Wallis检验;*、**、***和****分别代表P<0.05、P<0.01、P<0.001和P<0.0001。

图4.管胞长度的TWARE分析。（A）管胞长度TWAS曼哈顿图。图中标注部分与管胞长度显著相关的基因。（B）与TWAS中管胞长度显著相关的基因的GO富集。（C）植物型细胞壁生物生成（GO：0009832）和植物型细胞壁组装（GO：0071668）基因与管胞长度的相关分析。SCC 表示斯皮尔曼相关系数。

5 12个木材性状的全转录组关联分析

在本研究中，研究者从290个样本中剔除低表达基因（平均FPKM<0.5）后，保留了18203个基因，并使用GEMMA软件中实现的线性混合模型（LMM）进行TWAS（表S6，图S7）。在12个木材性状的TWAS中，微纤丝角（MFA）未鉴定出显著（FDR<0.01）关联基因，而对于CCD、TW和TWT，检测到的显著关联基因数量过多——分别有8638、6739和1931个基因（表S6，图S7）。因此，未对这些性状进行进一步的深入分析。对于其余八个木材性状，共鉴定出2333个基因表现出显著关联（表S6，图S7）。值得注意的是，木质素含量显示的显著关联基因最多（1335个），其次是管胞长度（471个基因）（表S6，图S7）。

研究者对与管胞长度显著关联的基因进行GO富集分析显示，除了与核糖体形成和翻译相关的术语外，其余基因主要富集在与细胞壁发育相关的过程，包括植物型细胞壁生物发生（GO:0009832）和植物型细胞壁组装（GO:0071668）（图4B）。这表明这些关联基因可能与细胞壁的合成和组装有关。由于管胞的发育和结构形成依赖于SCW组分的沉积和有序组装，研究者进一步提出这些关联基因可能参与木材管胞的形成。在这两个术语富集的11个基因中，8个与管胞长度呈显著正相关，而3个与管胞长度呈显著负相关（图4A和C）。然而，这些基因缺乏详细的功能注释；因此，它们在银杏细胞壁沉积和组装中的具体作用需要进一步研究。

通过TWAS，研究者鉴定了多个与木质素含量显著正相关的结构基因和转录因子基因，包括编码阿罗酸/预苯酸脱水酶（ADT/PDT）的evm.TU.chr3.317和evm.TU.chr4.1696；编码苯丙氨酸解氨酶（PAL）的evm.TU.chr12.1408；编码过氧化物酶（PER）的evm.TU.chr2.1128；以及编码MYB91的evm.TU.chr8.382和evm.TU.chr9.514（图5A-G）。此外，在显著关联的基因中，所有bHLH家族成员均显示与木质素含量显著正相关，包括evm.TU.chr3.1718,evm.TU.chr3.1923,evm.TU.chr3.2311,evm.TU.chr3.2340,evm.TU.chr4.104,evm.TU.chr4.1572,evm.TU.chr4.2052,和evm.TU.chr9.1343（图5H）。

图5.木质素含量的TWAS分析。（A）木质素含量TWAS曼哈顿图。图中标注部分与木质素含量显著相关的基因。（B-G）evm.TU.chr3.317(ADT/PDT)、evm.TU.chr4.1696(ADT/PDT)、evm.TU.chr12.1408(PAL)、evm.TU.chr2.1128(PER)、evm.TU.chr8.382(MYB91)和evm.TU.chr9.514(MYB91)与木质素含量的相关性分析。(H)八个bHLH基因与木质素含量的相关性分析。

对于纤维素含量，研究者共鉴定出363个与该性状显著相关的基因（图6A）。研究者对这些基因进行GO富集分析显示，主要富集在与多糖和碳水化合物的代谢及生物合成相关的过程（图6B）。在纤维素生物合成过程（GO:0030244）、纤维素代谢过程（GO:0030243）和植物型次生细胞壁生物发生（GO:0009834）等术语中也观察到显著富集（图6B）。这些发现强调所鉴定的基因可能参与纤维素生物合成和SCW发育。在富集于纤维素生物合成过程（GO:0030244）的基因中，6个中的4个被注释为纤维素合酶（CesA）基因，并显示出与纤维素含量的显著正相关（图6A,C-F）。研究者将这些基因与拟南芥中的CesA基因进行系统发育分析显示，evm.TU.chr1.2824和evm.TU.chr1.2825与AtCesA8同源，因此分别命名为GbCesA8A和GbCesA8B（图2E）。同样，evm.TU.chr10.692和evm.TU.chr12.1759分别与AtCesA4和AtCesA7同源，并命名为GbCesA4和GbCesA7（图2E）。所有四个CesA基因在银杏木质部中均表现出高表达水平，并显示出高度相关的表达模式（图6G）。此外，研究者鉴定了一个基因evm.TU.chr7.1702，编码WRKY12，其与纤维素含量呈负相关（图6H）。

在与半纤维素含量显著相关的139个基因中，鉴定出2个基因evm.TU.chr3.293和evm.TU.chr9.476——位于不同染色体上；两者都编码β-甘露聚糖合酶（CSLA），这是一种负责甘露聚糖生物合成的酶（图7A）。表达分析显示这些基因与半纤维素含量之间存在强正相关（图7D）。系统发育分析进一步表明，这两个基因都与拟南芥中的AtCSLA9同源，因此分别命名为GbCSLA9A和GbCSLA9B（图3E）。此外，一个转录因子基因evm.TU.chr5.1668(MYB46)显示出与半纤维素含量以及GbCSLA9A和GbCSLA9B的显著正相关（图7B-D）。

图6.纤维素含量的TWARE分析。（A）基于纤维素含量的TWAS曼哈顿图。图中标注部分与纤维素含量显著相关的基因。(B)TWAS中与纤维素含量显著相关基因的GO富集分析。（C-F）evm.TU.chr10.692(GbCesA4)、evm.TU.chr12.1759(GbCesA7)、evm.TU.chr1.2824(GbCesA8A)和evm.TU.chr1.2825(GbCesA8B)与纤维素含量的相关性分析。(G)四个CesA基因之间的相关性分析。(H)evm.TU.chr7.1702(WRKY12)与纤维素含量的相关性分析。

图7.半纤维素含量的TWAS分析。（A）基于半纤维素含量的TWAS曼哈顿图。图中标注部分与半纤维素含量显著相关的基因。（B，C）evm.TU.chr5.1668(MYB46)与evm.TU.chr3.293(GbCSLA9A)和evm.TU.chr9.476(GbCSLA9B)的相关性分析。(D)evm.TU.chr5.1668(MYB46)、evm.TU.chr3.293(GbCSLA9A)和evm.TU.chr9.476(GbCSLA9B)与半纤维素含量的相关性分析。

6 加权基因共表达网络分析及枢纽基因鉴定

研究者通过WGCNA共鉴定出14个基因表达模块，包含13465个基因（图8A）。另有4738个基因未被分配到任何已鉴定的模块中，因此被归类到灰色模块（未分类模块）（图8A）。在这14个模块中，暗蓝绿色（r=0.24,P=5e-05）、灰60（r=0.23,P=7e-05）和午夜蓝（r=0.21,P=3e-04）模块与木质素含量呈显著正相关，意味着这三个模块在一定程度上与木质素合成和积累过程相关（图8A）。此外，灰60、皇家蓝和暗灰色模块与各种管胞形态性状（如管胞长度（TL）、管胞宽度（TW）和管胞壁厚（TWT））显示出显著相关性，表明它们可能与管胞发育有潜在联系（图8A）。值得注意的是，黄色模块与纤维素含量（r=0.23,P=9e-05）、半纤维素含量（r=0.19,P=9e-04）和纤维素结晶度（r=0.39,P=3e-12）呈显著正相关，表明该模块可能与银杏SCW形成和木材结构特征的调控密切相关（图8A）。

为了进一步阐明黄色模块内基因的生物学功能，研究者选择TOM>0.1的基因构建共表达网络并进行GO富集分析（图8B和C）。分析显示，与SCW生物合成及相关过程密切相关的术语显著富集，包括纤维素的合成和代谢（GO:0030244,GO:0030243）、半纤维素（GO:0010410）、β-葡聚糖（GO:0051274,GO:0051273）、木聚糖（GO:0045491,GO:0045492）、细胞壁多糖（GO:0010383,GO:0070592）和细胞壁大分子（GO:0044036,GO:0044038）（图8B）。这些结果表明黄色模块可能在银杏SCW组分的合成和重塑中发挥重要作用。

值得注意的是，共表达网络突出了4个连通性最高的枢纽基因，即evm.TU.chr10.692(GbCesA4),evm.TU.chr12.1759(GbCesA7),evm.TU.chr1.2825(GbCesA8A),和evm.TU.chr1.2824(GbCesA8B)（图8D）。所有这些基因都编码纤维素合酶（CesA），并且与TWAS鉴定的四个CesA基因完全一致（图6A）。此外，网络中还包括几个与纤维素合成或重塑密切相关的基因，如编码α-微管蛋白的TUBA(evm.TU.chr2.1590)；编码β-微管蛋白的TUBB(evm.TU.chr9.972,evm.TU.chr10.2089,和evm.TU.chr4.2321)；以及编码内切葡聚糖酶的EG(evm.TU.chr4.1432和evm.TU.chr5.993)（图8D）。

在黄色模块内，鉴定出一个以甘露聚糖生物合成相关基因evm.TU.chr3.293(GbCSLA9A)和evm.TU.chr9.476(GbCSLA9B)；木聚糖生物合成相关基因evm.TU.chr11.1607（推定的β-1,4-木糖基转移酶，IRX14）和evm.TU.chr5.417（推定的β-1,4-木糖基转移酶，IRX9）；以及木聚糖修饰相关基因evm.TU.chr6.1385（可能的葡糖醛酸木聚糖葡糖醛酸基转移酶，IRX7）为中心的局部共表达网络（图8E）。该网络还包括多个关键功能基因，如木聚糖修饰相关基因evm.TU.chr1.1515(葡糖醛酸木聚糖4-O-甲基转移酶，GXMT)和evm.TU.chr3.1397(GXMT)，以及与UDP-糖前体生成和转运相关的基因，包括evm.TU.chr8.1928(UDP-葡糖醛酸脱羧酶，UXS)和evm.TU.chr2.1702(UXS)，evm.TU.chr9.7(UDP-阿拉伯糖4-差向异构酶，UXE)，以及evm.TU.chr10.1387(UDP-鼠李糖/UDP-半乳糖转运蛋白，URG7)（图8E）。此外，IRX14和IRX9与GbCSLA9A和GbCSLA9B显示出强正相关（图8F）。

在黄色模块中，研究者还鉴定了几个转录因子基因，包括evm.TU.chr1.637(GATA9)和evm.TU.chr10.213(HB15)，它们分别与纤维素和半纤维素含量呈显著正相关和负相关（图8G）。同时，发现evm.TU.chr11.910(ELF6)与管胞长度和管胞宽度呈显著正相关（图8H和I）。

图8.WGCNA及SCW发育中枢纽基因的鉴定。(A)模块特征基因与木材性状之间的相关性（每个单元格包含相关系数和P值）。(B,C)黄色模块中TOM>0.1的基因的GO富集分析(B)和共表达网络(C)。(D)以CesA基因为中心的共表达网络。(E)以CSLA和IRX基因为中心的共表达网络。(F)evm.TU.chr3.293(GbCSLA9A)、evm.TU.chr9.476(GbCSLA9B)、evm.TU.chr5.417(IRX9)和evm.TU.chr11.1607(IRX14)与半纤维素含量的相关性分析。(G)evm.TU.chr1.637(GATA9)和evm.TU.chr10.213(HB15)与纤维素和半纤维素含量的相关性分析。(H,I)evm.TU.chr11.910(ELF6)与管胞长度和管胞宽度的相关性分析。

木材性状是决定木材特性和应用潜力的关键因素。系统地鉴定和表征控制木材性状的关键基因，阐明木材形成的分子机制，对于提高木材品质和利用效率至关重要。银杏木材具有优异的材料特性和巨大的商业潜力。然而，与针叶树和大多数被子植物相比，研究者对其木材形成机制和参与SCW发育的调控基因仍缺乏全面的了解。为填补这一知识空白，研究者整合了GWAS、TWAS和WGCNA，基于290个银杏种质的全基因组重测序和转录组测序数据，鉴定了影响银杏木材性状的关键基因，并阐明了次生木质部发育的调控机制。

1 银杏全基因组SNP突变模式和转换偏差

在本研究中，利用平均测序深度为7.6×的全基因组重测序数据（表S3），研究者共鉴定出76128881个SNP和1974424个InDel（图1A和B），随后用于银杏12个木材性状的GWAS。全基因组SNP类型的统计分析显示，转换/颠换比（TS/TV）约为1.59（图1C），表明该群体中的核苷酸替换主要是化学性质相似的碱基之间的转换（嘌呤到嘌呤或嘧啶到嘧啶）。

值得注意的是，C:G>T:A转换的频率显著高于其他突变类型（图1C），这与在多种物种中的发现一致。已有充分证据表明，CpG二核苷酸背景下的胞嘧啶（C）常被甲基化为5-甲基胞嘧啶（5mC），后者易于自发脱氨成为胸腺嘧啶（T）。在单链DNA中，5mC更易脱氨并转化为胸腺嘧啶，且产生的胸腺嘧啶在缺乏互补链的情况下不会产生可识别的G:T错配，因此不会被DNA修复机制有效检测和修复。这些分子机制可能是全基因组范围内C:G>T:A转换显著富集的原因，为本研究中观察到的突变分布模式提供了合理的生物学解释。

2 与纤维素合成相关的候选基因

在纤维素含量的GWAS中，研究者鉴定出一个候选基因GbCesA8B(evm.TU.chr1.2825)。GbCesA8B内一个导致蛋白质功能完全丧失的移码插入（chr1:997452697）与纤维素含量降低显著相关，表明GbCesA8B可能对纤维素积累具有积极作用（图2D和F）。这一假设得到了TWAS结果的进一步支持，其中GbCesA8B的表达水平与纤维素含量呈显著正相关（图6F）。值得注意的是，在TWAS中与纤维素含量正相关的93个基因中，还鉴定了另外3个CesA基因：GbCesA4(evm.TU.chr10.692)、GbCesA7(evm.TU.chr12.1759)和GbCesA8A(evm.TU.chr1.2824)（图6C-E）。

研究者进一步发现，银杏中的GbCesA4、GbCesA7和GbCesA8A/B分别与拟南芥中的AtCesA4、AtCesA7和AtCesA8高度同源（图2E）。先前的研究表明，SCW中的纤维素生物合成由CesA4-CesA7-CesA8复合物催化，其中每个亚基都是不可或缺的。该模型在毛果杨中也得到了验证，删除PtrCesA4、PtrCesA7和PtrCesA8中的任何一个基因都会导致匍匐生长和纤维素含量急剧下降。有趣的是，拟南芥中的AtCesA4、AtCesA7和AtCesA8在木质部中特异性高表达，并表现出一致的表达模式，表明这三个CesA亚基在木质部中共表达，并以固定比例形成CesA4-CesA7-CesA8复合物。

与此一致，研究者观察到银杏中的GbCesA4、GbCesA7、GbCesA8A和GbCesA8B不仅在木质部组织中高表达，而且它们的表达水平之间存在强相关性（图6G）。这表明银杏中SCW纤维素的生物合成可能具有一定程度的进化保守性，并且可能需要像拟南芥和杨树等被子植物中观察到的那样，需要CesA4、CesA7和CesA8的协同参与。WGCNA的结果进一步支持了这一推测。WGCNA鉴定的黄色模块与纤维素含量呈显著正相关，其枢纽基因恰好对应于TWAS中检测到的四个CesA基因，即GbCesA4、GbCesA7、GbCesA8A和GbCesA8B（图8D），这四个CesA枢纽基因进一步与几个与纤维素生物合成密切相关的基因形成了共表达网络，包括TUBA(evm.TU.chr2.1590)、TUBB(evm.TU.chr9.972,evm.TU.chr10.2089,和evm.TU.chr4.2321)和EG(evm.TU.chr4.1432和evm.TU.chr5.993)（图8D）。TUBA和TUBB分别编码α-微管蛋白和β-微管蛋白，它们组装成皮层微管，引导纤维素合酶复合物（CSC）的运动，并决定纤维素微纤丝沉积的方向和细胞伸长。EG编码内切葡聚糖酶，参与细胞壁内纤维素的降解和重塑。综合来看，研究者推测CesA与TUBA/TUBB之间可能存在潜在的共表达机制，这可能有助于维持CSC与皮层微管之间的协调相互作用，以确保纤维素微纤丝的有序沉积。EG在纤维素合成相关基因模块中的共表达表明，动态重塑可能与SCW形成和纤维素生物合成同时发生。统计关联和共表达网络分析的汇聚为这些基因在纤维素生物合成中的关键作用提供了有力支持。

此外，研究者通过TWAS和WGCNA鉴定了多个可能参与调控纤维素生物合成的转录因子，为理解银杏SCW形成的调控网络提供了额外见解。WRKY12(evm.TU.chr7.1702)的鉴定表明它可能在银杏纤维素生物合成中发挥负调控作用（图6H）。AtWRKY12在拟南芥中的功能已被充分表征：AtWRKY12的功能缺失突变导致NST2表达激活和髓细胞中SCW的异位沉积，而AtWRKY12的同源物MIRWRKY12的过表达则恢复了AtWRKY12突变体髓细胞的SCW增厚缺陷。这些发现共同支持WRKY12作为SCW形成抑制因子的作用。与此一致，研究者的TWAS结果显示，银杏中WRKY12表达与纤维素含量呈显著负相关，表明WRKY12也可能参与银杏木质部SCW发育的负调控。

值得注意的是，在以GbCesA4、GbCesA7、GbCesA8A和GbCesA8B为枢纽基因的纤维素生物合成相关共表达网络中，研究者还鉴定了GATA9(evm.TU.chr1.637)和HB15(evm.TU.chr10.213)，它们分别与纤维素/半纤维素含量呈显著正相关和负相关（图8G）。这一观察结果与先前的研究一致。在毛果杨中，GATA9的过表达导致NAC结构域转录因子和下游SCW生物合成基因的上调，从而促进束间纤维的SCW增厚；相反，在拟南芥中，HB15的过表达抑制SCW形成并扰乱维管发育。基于这些观察，研究者认为GATA9和HB15可能分别在银杏木质部SCW沉积中起正向和负向作用。然而，它们的确切作用机制和影响程度仍有待进一步的实验验证来阐明。

如前所述，AtCesA4、AtCesA7和AtCesA8编码SCW形成过程中纤维素生物合成所需的纤维素合酶异构体。在毛果杨中，5个同源基因PtCesA4、PtCesA7A、PtCesA7B、PtCesA8A和PtCesA8B被鉴定为SCW纤维素沉积的贡献者。通过基因敲除进行的功能表征显示，PtCesA4的敲除突变体和PtCesA7a/b或PtCesA8a/b的双敲除突变体表现出匍匐生长习性和纤维素含量显著降低，而PtCesA7a、PtCesA7b、PtCesA8a或PtCesA8b的单敲除未显示明显的生长缺陷。这一观察结果表明PtCesA7A与PtCesA7B之间，以及PtCesA8A与PtCesA8B之间存在功能冗余，使得在纤维素合成过程中一个拷贝的丢失可以得到补偿。相比之下，PtCesA4与PtCesA7A/B和PtCesA8A/B功能相似但不可冗余，在纤维素合成中发挥重要作用。

整合功能注释和系统发育推断，研究者在银杏中鉴定了4个CesA同源基因，即GbCesA4、GbCesA7、GbCesA8A和GbCesA8B，它们被认为参与SCW纤维素生物合成（图2E）。研究者还观察到，携带导致基因功能完全丧失的GbCesa8B插入突变的个体仅显示出轻微的纤维素含量降低，而没有明显的生长缺陷（图2F）。考虑到在毛果杨中观察到的PtCesa8A和PtCesa8B之间的功能冗余，研究者假设银杏中的GbCesa8A和GbCesa8B也可能表现出功能冗余，GbCesa8B的丢失可能由冗余基因GbCesa8A补偿。这一假设需要通过基因敲除实验进一步验证。值得注意的是，在GbCesa4或GbCesa7中未检测到导致功能完全丧失的突变。这一观察结果表明，涉及GbCesa4和GbCesa7的纤维素合成可能缺乏替代的冗余基因，因此这些基因在进化过程中受到更强的纯化选择。

3 与半纤维素合成相关的候选基因

在半纤维素含量的GWAS中，研究者鉴定出一个候选基因GbCSLA9B(evm.TU.chr9.476)。该基因内的一个非移码缺失（chr9:94004237）导致第41位亮氨酸丢失，并与半纤维素含量降低显著相关（图3D和F）。此外，GbCSLA9B与拟南芥中的AtCSLA9显示出高序列同源性（图3E）。AtCSLA9编码一种定位于高尔基体膜的β-甘露聚糖合酶。该蛋白通常包含奇数个跨膜结构域，其催化活性位点朝向高尔基体腔，并催化半纤维素中β-1,4-甘露聚糖或葡甘露聚糖骨架的合成。

一致地，InterPro结构域预测表明GbCSLA9B也是一种具有多个跨膜结构域的膜结合蛋白，与半纤维素含量降低相关的第41位亮氨酸缺失位于跨膜结构域内。因此，研究者推测GbCSLA9B可能在银杏的甘露聚糖生物合成中也发挥重要作用。缺失突变体中观察到的半纤维素含量降低的一个可能解释是，跨膜区氨基酸的丢失可能影响GbCSLA9B的膜结合能力，从而限制其在甘露聚糖/葡甘露聚糖合成中的催化功能。

此外，TWAS发现GbCSLA9B和GbCSLA9A都与纤维素含量呈显著正相关，并且两者都与编码拟南芥中β-甘露聚糖合酶的AtCSLA9显示出高序列同源性（图7A和3E）。这进一步支持了研究者的推测，即GbCSLA9A和GbCSLA9B参与甘露聚糖/葡甘露聚糖生物合成。值得注意的是，GbCSLA9A和GbCSLA9B与IRX14(evm.TU.chr11.1607)和IRX9(evm.TU.chr5.417)一起形成了一个以这些基因为中心的局部共表达网络（图8E）。该网络还包括多个已知参与半纤维素生物合成的基因，包括UXS(evm.TU.chr2.1702)、UXE(evm.TU.chr9.7)、URGT(evm.TU.chr10.1387)、IRX7(evm.TU.chr6.1385)和GXMT(evm.TU.chr1.1515和evm.TU.chr3.1397)。该共表达网络所属的黄色模块与半纤维素含量呈显著正相关（图8A和E）。

在该网络中，枢纽基因IRX14和IRX9已被报道通过β-1,4-糖苷键催化木糖残基添加到木聚糖骨架上。有趣的是，在GbCSLA9A/9B与IRX14/IRX9之间观察到强相关性（图8F）。这一观察结果表明，尽管不同银杏个体SCW中半纤维素含量存在相当大的变异，但两种主要的半纤维素组分甘露聚糖和木聚糖的生物合成可能表现出协同关系，确保甘露聚糖和木聚糖的沉积比例保持在一定的范围内。

此外，该共表达网络中的其他基因参与半纤维素生物合成的多个阶段。UXS催化UDP-葡萄糖转化为UDP-木糖，为木聚糖骨架延伸提供活化的木糖供体。UXE将UDP-木糖转化为UDP-阿拉伯糖，为木聚糖的阿拉伯糖基化提供底物。IRX7参与葡糖醛酸木聚糖/O-乙酰葡糖醛酸木聚糖（GX/ACGX）生物合成所需的还原端四糖的形成，该过程对于正常的GX/ACGX合成和提供骨架延伸的关键起始结构至关重要。在拟南芥irx7突变体中，ACGX含量降低超过60%。GXMT催化GX中葡糖醛酸取代基的4-O-甲基化，GXMT1突变体中GX的4-O-甲基化程度比野生型降低约75%。URGT介导UDP-半乳糖从胞质溶胶转运到高尔基体腔，从而支持甘露聚糖骨架的半乳糖基化。综上所述，这个以GbCSLA9A/9B和IRX14/9为中心的共表达网络可能负责银杏中的半纤维素生物合成，涵盖从前体供应、核苷酸糖转运到多糖聚合和合成后修饰的多个步骤。

研究者还鉴定了一个转录因子基因MYB46(evm.TU.chr5.1668)，其表达水平不仅与半纤维素含量显著正相关，而且与GbCSLA9A/9B也显著正相关（图7）。MYB46是参与SCW形成的NAC-MYB级联网络中的次级主调控因子，并已被证明可以激活被子植物中与SCW生物合成相关的一系列结构基因和转录因子。值得注意的是，先前的研究表明，在拟南芥中，MYB46直接结合CSLA9启动子并驱动其转录活性。因此，研究者推测可能存在一个保守的MYB46-CSLA9调控模块，参与银杏木质部中甘露聚糖的合成。

4 与木质素合成相关的候选基因

关于木质素含量，研究者鉴定了两个基因evm.TU.chr3.317和evm.TU.chr4.1696，它们编码ADT/PDT，并且都与木质素水平显著相关（图5B-C）。ADT/PDT是苯丙氨酸合成的关键限速酶，为苯丙烷代谢途径提供初始底物。ADT/PDT可以直接催化L-阿罗酸脱水生成苯丙氨酸，或者先催化预苯酸脱水和脱羧生成苯丙酮酸，然后通过转氨作用转化为苯丙氨酸。由于苯丙氨酸作为苯丙烷途径的初始底物，其生物合成效率直接决定了通过该代谢途径的碳通量，最终影响木质素的潜在积累。除了前体合成外，TWAS还鉴定了苯丙烷途径的几个下游关键基因，如evm.TU.chr12.1408(PAL)和evm.TU.chr2.1128(PER)，它们都与木质素含量显著相关（图5D-E）。PAL作为该途径的第一个限速酶，将苯丙氨酸脱氨生成肉桂酸，从而为木质素单体生物合成提供底物。相比之下，PER作用于终端步骤，直接催化木质素单体的聚合和交联，这直接决定了SCW的木质素沉积。重要的是，在研究者的结果中，evm.TU.chr3.317、evm.TU.chr4.1696、evm.TU.chr12.1408和evm.TU.chr2.1128都与木质素含量呈显著正相关（图5A-E）。这些发现表明，TWAS鉴定的结构基因覆盖了木质素生物合成的多个关键步骤，包括初始底物供应、底物转化和单体聚合。值得注意的是，研究者还鉴定了两个MYB91基因（evm.TU.chr8.382和evm.TU.chr9.514），它们与木质素含量呈正相关（图5F-G）。先前的研究表明，MYB转录因子通过激活多种结构基因以增强木质素积累，在调控木质素生物合成中发挥核心作用。特别是，MYB91在梨中瞬时过表达时已被证明能显著提高木质素水平。这与研究者在银杏中的发现高度一致，表明MYB91也可能参与银杏的木质素积累。

此外，研究者鉴定了8个bHLH基因，它们与木质素含量呈显著正相关，包括evm.TU.chr3.1718(bHLH6)、evm.TU.chr3.2340(bHLH29)、evm.TU.chr3.1923(bHLH144)、evm.TU.chr4.1572(bHLH80)、evm.TU.chr9.1343(bHLH130)、evm.TU.chr3.2311、evm.TU.chr4.104和evm.TU.chr4.2052（图5H）。这一发现表明bHLH转录因子可能作为银杏木质素生物合成的正调控因子。来自其他物种的累积证据已经确定了多个bHLH转录因子在木质素调控中的关键作用。例如，在毛白杨、陆地棉、白菜、本氏烟和银白杨×腺毛杨中，bHLH家族的多个成员已被证明可以激活参与木质素生物合成的关键基因的转录，包括CCoAOMT、PER和POD。这种转录激活最终促进了木质素的积累和沉积，突出了bHLH因子在木质素生物合成中的保守调控作用。

值得注意的是，木质素和类黄酮生物合成途径竞争共同前体对香豆酰辅酶A。查尔酮合酶（CHS）是类黄酮生物合成途径中的关键酶，催化对香豆酰辅酶A和丙二酰辅酶A缩合形成查尔酮，这是类黄酮的前体。因为对香豆酰辅酶A也是木质素生物合成的关键前体，CHS活性关键地影响木质素和类黄酮生产之间的平衡。有趣的是，在苹果中，MhbHLH130通过直接结合其启动子来抑制MhCHS的表达，从而使更多的碳通量流向木质素生物合成。与该机制一致，研究者发现evm.TU.chr9.1343(bHLH130)在银杏中与木质素含量呈显著正相关。因此，研究者推测银杏中的bHLH130可能通过对香豆酰辅酶A的代谢命运调控，对木质素积累产生积极影响。

5 与木材密度和管胞发育相关的候选基因

在与WOD相关的区域中发现了一个基因evm.TU.chr3.1073(TAA1)，编码L-色氨酸-丙酮酸氨基转移酶（表S5）。TAA1催化氨基从L-色氨酸转移到丙酮酸，生成吲哚-3-丙酮酸（IPA），这是吲哚-3-乙酸（IAA）的前体。该途径是植物中IAA生物合成的主要途径。已知IAA通过ARF转录因子调控形成层细胞的维持和分化。因此，研究者提出evm.TU.chr3.1073可能通过影响IAA合成来影响形成层细胞的维持或分化，从而可能影响木材密度。

与WBD相关的基因evm.TU.chr10.1203(XYL4)和evm.TU.chr7.1369(MAN)，分别编码木聚糖1,4-β-木糖苷酶和甘露聚糖内切-1,4-β-甘露糖苷酶（表S5）。XYL4和MAN在木材发育过程中半纤维素的降解和SCW结构的重塑中起关键作用。这些结果表明银杏的木材密度可能与SCW重塑有关。

管胞形态性状似乎与evm.TU.chr6.1080(EIN3)和evm.TU.chr6.1069(ABI5)有潜在关联（表S5）。EIN3和ABI5分别编码ETH和ABA的响应因子，在相应的信号转导途径中发挥关键作用。已发现ABA和ETH通过促进拟南芥和棉花纤维中H₂O₂的产生来调节细胞伸长和SCW合成。在陆地棉纤维细胞发育的早期阶段，低浓度的H₂O₂促进细胞伸长。然而，在过渡阶段，ABA和ETH可能通过调节高水平的H₂O₂和钙离子来抑制纤维伸长并促进SCW合成。因此，研究者推测EIN3和ABI5可能通过ETH和ABA信号转导影响银杏木质部管胞细胞的发育。

此外，研究者在与纤维素/半纤维素合成相关的黄色模块中鉴定了evm.TU.chr11.910(ELF6)，其表达与管胞长度和宽度呈显著正相关（图8H-I）。先前的研究表明，ELF6缓解了拟南芥花器官中扩展蛋白基因的表观遗传沉默效应，从而促进心皮细胞的伸长。根据本研究的结果，推测ELF6可能通过类似的机制影响银杏木质部的管胞形态。

在本研究中，研究者整合了GWAS、TWAS和WGCNA，以阐明银杏木材性状变异和SCW发育的遗传基础和调控机制。多条证据，包括显著的GWAS信号、表达-性状关联、共表达网络和系统发育分析，一致支持GbCesA4、GbCesA7、GbCesA8A和GbCesA8B在银杏SCW纤维素生物合成中的协同核心作用。这些CesA基因与微管相关基因（TUBA和TUBB）和内切葡聚糖酶基因（EG）的共表达表明纤维素生物合成、细胞骨架引导和动态细胞壁重塑之间存在紧密协调。此外，GbCesA8B功能缺失突变的轻微表型下降表明，银杏中的GbCesA8A和GbCesA8B可能具有功能冗余，类似于毛果杨中的PtrCesA8A和PtrCesA8B，而GbCesA4和GbCesA7中缺乏功能缺失突变可能是由于缺乏冗余基因，使得这些基因受到更强的纯化选择。除了结构基因，TWAS和WGCNA还鉴定了几个可能参与SCW发育调控的转录因子，其中GATA9可能作为正调控因子，而WRKY12和HB15可能作为SCW形成的负调控因子。

对于半纤维素生物合成，研究者鉴定出GbCSLA9A/9B是甘露聚糖/葡甘露聚糖合成的主要贡献者，与木聚糖生物合成基因（包括IRX9和IRX14）形成一个紧密连接的共表达网络。该网络涵盖参与半纤维素生物合成多个阶段的基因，如UXS、UXE、IRX7、GXMT和URG7，从核苷酸糖的供应和运输到聚合物的延伸和修饰。这表明一个协调机制，维持SCW中主要半纤维素组分的平衡沉积。此外，MYB46与半纤维素含量和CSLA9的显著关联进一步支持了可能存在一个调控银杏中甘露聚糖合成的MYB46-CSLA9模块。

参与苯丙氨酸生物合成（ADT/PDT）、苯丙氨酸代谢（PAL）和木质素聚合（PER）的关键基因都与木质素含量显著相关，表明银杏中的木质素生物合成似乎通过前体供应和下游聚合过程进行调控。此外，包括MYB91和几个bHLH家族成员在内的多个转录因子被认为对木质素积累有积极影响，其中bHLH130可能通过调节木质素和竞争代谢途径之间的碳通量分配来影响木质素积累。

往期精彩回顾：