全转录组共表达网络分析鉴定大豆种子蛋白质含量相关候选基因

《Agriculture》：Transcriptome-Wide Co-Expression Network Analysis Identifies Candidate Genes Associated with Seed Protein Content in Soybean

　　

大豆中与最终种子蛋白质含量相关的转录过程尚未被充分表征，特别是在与野生大豆杂交的后代中。为填补这一知识空白，研究人员在开花后25、35和45天（D25、D35和D45）比较高蛋白质（HP）和低蛋白质（LP）BC1F3后代，采用种子成分测定、RNA测序（RNA-seq）、差异表达分析、时序表达聚类、功能富集以及加权基因共表达网络分析（WGCNA）。HP和LP表型在种子发育晚期差异最为显著：在D45时，HP的蛋白质含量高出9.0%，而油脂含量低7.3%。HP与LP之间的差异表达基因（DEGs）数量从D25时的3386个增加到D45时的7269个，有18,394个基因表现出显著的Group × Time互作。六个时序聚类揭示了阶段依赖性的表达模式，这些模式在细胞分裂、光合作用、储备相关功能、翻译、蛋白质加工以及氧化还原相关过程中富集。WGCNA鉴定出与蛋白质含量正相关的pink和salmon模块，以及一个与蛋白质含量负相关的brown模块。整合阶段特异性DEGs、时序互作和网络指标，优先筛选出候选基因，包括Glyma.11G063900.1、Glyma.06G030400.1、Glyma.07G038500.2和Glyma.03G065500.2。与油脂含量相关的模块主要被解释为发育模块，因为油脂含量与发育阶段几乎共线。这些结果提供了一个阶段分辨的框架，旨在为功能验证和旨在提高大豆种子蛋白质含量的育种确定候选基因优先级。

论文解读

大豆（Glycine max）种子成分在食品、饲料和加工领域具有重要经济价值，蛋白质与油脂含量共同决定其商业价值。这两个性状遗传复杂且普遍呈负相关，构成了长期存在的育种挑战。最终种子成分并非单一途径或发育事件的结果。蛋白质积累依赖氮获取与同化、氨基酸代谢、储藏蛋白合成、折叠、运输与周转，而油脂积累依赖碳供应、糖酵解与质体代谢、脂肪酸合成、甘油三酯组装及油体形成。碳氮通量在种子填充过程中紧密耦合，影响运输和源-库分配的基因可同时改变两种性状。发育转录组学研究表明，这些过程的时序与平均表达同等重要。野生大豆作为栽培大豆的近缘祖先，保留了驯化过程中减少或重塑的遗传多样性，为研究种子蛋白质含量的分子基础提供了有用的遗传资源。在此背景下，研究人员提出科学问题：高蛋白表型是否与碳氮代谢、储藏过程和蛋白质稳态相关基因的阶段特异性转录差异有关，而非简单地高表达储藏蛋白基因？

为回答这一问题，研究人员以栽培大豆品种Daepung（Glycine max）与高蛋白野生大豆GWS-1887（Glycine soja）杂交构建的BC₁F₃回交后代为样本队列来源，从中选取蛋白质含量差异显著的高蛋白（HP）和低蛋白（LP）株系。2023年在温室种植条件下，于开花后25、35和45天（D25、D35、D45）采集发育种子，每个株系每个阶段设置3个独立混合生物学重复，共18个样本。研究采用的主要技术方法包括种子成分测定、RNA测序（RNA-seq）、差异表达分析、时序表达聚类、功能富集分析以及加权基因共表达网络分析（WGCNA）。种子蛋白质含量采用凯氏定氮法测定，油脂含量采用索氏提取法测定。

研究结果显示，HP和LP株系的蛋白质和油脂含量在发育过程中逐渐分化。通过SoySNP6K Illumina BeadChip芯片分析确认，两株系遗传相似度约93.5%，遗传差异主要局限于染色体20的POWR1基因座周围，这一高遗传相似性减少了背景遗传变异的混杂效应。D25时HP蛋白质含量为47.05%，LP为46.32%；D35时HP为46.07%，LP为44.37%；D45时HP为47.96%，LP为44.00%，差异3.96个百分点，HP约高9.0%，且经Holm调整后差异仍显著。D45时HP油脂含量比LP低约7.3%。双向方差分析（ANOVA）表明，株系、发育阶段及其互作对两种性状均有显著影响。

RNA-seq产生18个双端文库，每个样本测序深度约12.84至22.05百万reads，平均比对率为88.36%。基于500个变异最大的方差稳定化变换（VST）基因的主成分分析（PCA）显示，样本首先按发育阶段分离，其次按HP/LP分组，PC₁和PC₂分别解释83%和10%的方差。HP与LP之间的DEGs数量随发育急剧增加：D25时3386个，D35时4083个，D45时7269个；三个阶段共有11,703个基因至少在比较中一个差异表达，仅305个在所有三个阶段共享，表明转录差异具有强烈的阶段依赖性。HP内部发育比较（D25至D45）鉴定出14,966个DEGs，LP为8477个，表明HP的转录组重编程更为广泛。Group × Time似然比检验鉴定出18,394个具有显著株系依赖性发育轨迹的基因，占分析基因的54.4%。

研究人员将18,394个显著互作基因按平均VST表达的z-score进行k-means聚类（k=6），形成六个时序聚类。C1（5220个基因）在HP中后期显著升高，富集于自噬体组装、胚胎发育、泛素介导蛋白水解等。C2（2205个基因）在HP中向D45升高而在LP中下降，富集于核糖体、蛋白酶体、蛋白质折叠、氧化磷酸化等。C3（3474个基因）在HP中D25最高，富集于微管运动、DNA复制等。C4（2172个基因）在HP中D25至D35强烈增加，富集于营养库活性、脂质储存体成分、内质网蛋白加工等。C5（2338个基因）在LP中D25后增加而HP中下降，富集于激素响应、ABC转运蛋白等。C6（2985个基因）呈现交叉模式，富集于光合作用、糖酵解等。这些聚类表明，HP与LP在增殖性、光合、储藏、运输和蛋白质稳态程序的时序和持续性上存在差异。

对33,817个保留基因进行WGCNA分析，鉴定出14个彩色模块。pink模块（517个基因，r=0.829）和salmon模块（61个基因，r=0.836）与蛋白质含量正相关，brown模块（4124个基因，r=?0.852）与蛋白质含量负相关。五个模块与油脂含量显著相关，但油脂含量与发育时间几乎完全相关（r=0.996），因此这些模块主要反映发育进程而非HP-LP特异性差异。八个显著模块-性状对中4066个基因达到hub候选标准。

阶段特异性功能富集分析显示，D25时HP高表达基因富集于微管运动、DNA复制等细胞分裂相关功能；D35时HP高表达基因富集于核糖体、光合作用、光捕获、电子传递和ATP合成等，LP高表达基因富集于植物激素信号、ABC转运蛋白、苯丙素和黄酮类生物合成等；D45时HP高表达基因富集于胚胎发育、翻译起始、蛋白质二硫键还原酶活性、谷胱甘肽代谢等，LP高表达基因富集于转运、ABC转运蛋白、淀粉和蔗糖代谢、碳固定等。

整合证据优先筛选候选基因。HP与LP之间DEGs并集包含11,703个基因，显著WGCNA模块包含27,275个基因，两者交集10,317个基因，完整hub候选集包含4066个基因。pink模块中Glyma.11G063900.1（推测的香叶基香叶基二磷酸合酶1，GGPS1）kME=0.991，在D35和D45于HP中表达更高；Glyma.06G030400.1（质体ADP-葡萄糖焦磷酸化酶大亚基候选）同样在HP中表达更高。salmon模块中Glyma.07G038500.2（萌发蛋白样蛋白1同源物）在D45于HP中表达更高。brown模块中Glyma.03G065500.2（SWEET样/Nodulin MtN3转运蛋白）在D35和D45于HP中表达较低。

讨论部分指出，HP与LP在发育过程中逐渐分化，转录差异也逐步增加，这一模式与近期大豆研究一致，即与蛋白质相关的最强转录差异出现在中晚期种子填充阶段。六个时序聚类表明，HP表型与细胞分裂、光合、储藏、运输和蛋白质稳态程序的时序差异相关。D35时HP中高表达的光合作用相关基因应解释为HP中这些基因表达维持或更强，而非净光合通量的直接证据。高蛋白低油脂的晚期表型符合种子储备的负相关关系，但转录组结果指向分布式的资源分配而非单一拮抗开关。高SNP标记一致性表明两株系遗传背景高度相似。POWR1本身在HP与LP之间没有显著转录丰度差异，因此染色体20的多态性不应被解释为POWR1差异表达直接解释HP表型。研究限制在于未进行RT-qPCR、独立蛋白或酶测定、代谢组学分析或代谢通量测量。

结论部分翻译：在本研究中，研究人员对两个相关的回交衍生大豆株系进行了比较，种子蛋白质含量和转录组特征随发育阶段分化，在D45时观察到最大的蛋白质含量差异。差异表达分析、Group × Time分析、时序聚类、功能注释和WGCNA鉴定出与种子蛋白质含量相关的基因和模块。富集的功能类别包括光合作用相关、翻译相关、蛋白质加工相关、转运相关和氧化还原相关注释。到D45时，HP在翻译、蛋白质加工、氧化还原和成熟相关过程中表现出更强的富集，而LP在转运和碳代谢特征方面表现出更强的富集。这些观察结果支持了一种发育资源分配模型，而不是单一途径解释。RNA-seq指标的整合优先筛选了候选基因。与油脂含量相关的模块被保留为发育背景，因为油脂与发育阶段几乎共线，所以不被解释为HP-LP特异性油脂调节的证据。这些结果为独立表达确认、功能验证和其他基因型的评估提供了假设和候选靶点；它们并不确立蛋白质积累的种子因果作用。