基于FarmCPU与机器学习方法在田间条件下鉴定大豆锈病部分抗性相关基因组区域

《Plants》：Identification of Genomic Regions for Partial Resistance to Soybean Rust Under Field Conditions Using FarmCPU and Machine Learning Approaches António Daniel Pedro Maquil, Tonny Obua, David L. Nsibo, Mildred Ochwo-Ssemakula, Harun Murithi, Paul Gibson, Ana Luísa Garcia-Oliveira, Richard Edema, Isaac Dramadri and Phinehas Tukamuhabwa + 1 author

　　

由真菌Phakopsora pachyrhizi引起的大豆锈病（Soybean Rust, SBR）威胁着全球大豆生产，可导致高达80%的产量损失。小种特异性Rpp基因由于病原体变异只能提供短期抗性，而部分抗性（Partial Resistance, PR）则提供了持久、广谱的保护，但其遗传基础尚不清楚。本研究旨在利用固定和随机模型循环概率统一（Fixed and Random Model Circulating Probability Unification, FarmCPU）全基因组关联分析（Genome-Wide Association Study, GWAS）和机器学习（Machine Learning, ML）方法——随机森林（Random Forest, RF）与支持向量回归（Support Vector Regression, SVR），来鉴定与PR相关的基因组区域和候选基因。研究人员在乌干达六个环境的自然感染条件下，对一个包含312份大豆种质的群体进行了评估。研究人员利用由锈病严重程度和孢子形成水平推导出的锈病指数（Rust Index, RI）来估算遗传力（H2），并通过最佳线性无偏预测（Best Linear Unbiased Predictions, BLUPs）对基因型进行排序，同时使用最佳线性无偏估计（Best Linear Unbiased Estimators, BLUEs）作为GWAS的输入。经过质量控制，研究人员在±60千碱基对（kb）的连锁不平衡（Linkage Disequilibrium, LD）窗口内分析了8272个单核苷酸多态性（Single-Nucleotide Polymorphisms, SNPs）。RI的多环境分析（Multi-environmental Analysis, MEA）显示了显著的遗传效应（p< 0.01）；H2= 0.57–0.68。研究人员共检测到61个位点：其中6个由FarmCPU检测到，15个由RF检测到，41个由SVR检测到。关键基因包括Glyma.01G128100（WRKY转录因子）、Glyma.13G228000（类受体激酶）和Glyma.20G173100（WD40结构域调节因子）。将机器学习与GWAS结合使用改善了位点检测，证实了PR的多基因性质，并支持了将基因组选择和位点聚合（pyramiding）应用于实现持久锈病抗性。

论文《Identification of Genomic Regions for Partial Resistance to Soybean Rust Under Field Conditions Using FarmCPU and Machine Learning Approaches》解读

一、 研究背景、问题与目的

大豆锈病（Soybean Rust, SBR）是由真菌Phakopsora pachyrhizi引起的一种毁灭性叶部病害，在热带和亚热带地区尤为严重，不防治可导致产量损失超过80%。目前，杀菌剂是主要的防治手段，但这增加了生产成本、带来环境风险并可能促进抗药性病原群体的出现。因此，寄主植物抗性是长期控制SBR最可持续和环境友好的策略。

目前已鉴定出十个小种特异性抗性位点（Rpp1-Rpp7等），但它们仅对特定病原小种有效，且常因病原群体毒性多样性和快速进化而失效。虽然聚合多个Rpp基因可延长抗性持久性，但这依赖于对当地病原小种组成的详细了解，而在非洲，由于病害监测能力有限，这一信息常常难以获取。因此，育种工作日益关注部分抗性（Partial Resistance, PR），这是一种由多个具有微小个体效应的位点控制的、数量遗传的性状。PR通过限制病斑发展、延迟症状进程和抑制孢子形成来降低病害严重程度，具有持久、广谱保护的潜力。然而，PR的遗传基础仍不明确。早期利用简单序列重复（Simple Sequence Repeat, SSR）标记的双亲数量性状位点（Quantitative Trait Locus, QTL）作图研究分辨率不足。基因分型测序（Genotyping-by-Sequencing, GBS）能够经济高效地进行全基因组单核苷酸多态性（SNPs）发现，从而提高全基因组关联分析（GWAS）对复杂性状（如PR）的分辨率。然而，针对SBR的PR进行的GWAS研究仍然很少。传统的GWAS方法通常应用严格的显著性阈值，单独测试SNP，这降低了对小效应变异的敏感性，并忽略了上位性相互作用，导致“缺失遗传力”问题。机器学习（ML）方法，如随机森林（RF）和支持向量回归（SVR），为建模复杂的非线性高维生物数据提供了强大工具，能够捕获上位性相互作用并识别常被传统GWAS方法遗漏的小效应变异。

针对以上问题，本研究旨在利用FarmCPU GWAS模型与RF、SVR机器学习方法，在一个多样化的乌干达田间大豆种质群体中，鉴定与SBR部分抗性相关的基因组区域和候选基因，以阐明其遗传结构，并为培育持久抗性大豆品种提供遗传资源和分子标记基础。本论文发表在期刊《Plants》上。

二、 主要研究方法概要

研究人员开展了一项田间与生物信息学相结合的研究。首先，他们在乌干达六个自然感染SBR的环境下，对一个包含312份来自不同地理起源（乌干达、美国、台湾、日本、津巴布韦、尼日利亚等）大豆种质的群体进行了表型评估，计算了锈病指数（RI）。通过多环境混合模型分析评估了表型变异、遗传力、基因型与环境互作等参数，并使用最佳线性无偏估计（BLUEs）和预测（BLUPs）进行基因型排序和GWAS输入。在基因型方面，研究人员利用DArTseq平台对群体进行基因分型测序（GBS），经质量控制后获得8272个高质量SNPs用于后续分析。通过主成分分析（PCA）和亲缘关系矩阵分析评估了群体结构。核心分析整合了三种方法：1） 基于混合线性模型的FarmCPU GWAS；2） 基于特征重要性的随机森林（RF）分析；3） 基于支持向量回归（SVR）的机器学习分析。最后，研究人员在连锁不平衡（LD）定义的窗口内（±60 kb）挖掘候选基因，并进行了基因本体（GO）功能富集分析。

三、 研究结果

1. 表型分析

研究人员对312份大豆基因型在六个乌干达环境下的锈病指数（RI）变异进行了评估。基因型分为两个成熟期组：MG1（80份，全为美国起源）和MG2（232份，多地理起源）。日本起源的基因型RI值最高，而美国起源的基因型最低。基于BLUPs的基因型排序显示，在MG1和MG2中分别鉴定出高PR和感病的极端基因型。多环境分析（MEA）表明，MG2中基因型、环境及基因型与环境互作（G × E）效应均显著，而在MG1中G × E互作不显著。遗传力（H²）估计值为0.57（MG1）至0.67（MG2），表明PR表达主要受遗传因素控制。

2. 全基因组SNP图谱与群体结构

经质量控制后，8272个高质量SNPs被保留用于分析。SNP在染色体上的数量分布不均，平均次要等位基因频率（MAF）为0.246。主成分分析（PCA）显示，PC1和PC2分别解释了15.3%和7.4%的遗传变异，成功将群体按成熟期组（MG1和MG2）和地理起源区分开。亲缘关系矩阵分析进一步识别出六个亚群（SP1–SP6），与PCA聚类模式一致。群体间遗传分化（FST）分析显示，MG1与MG2之间的平均FST为0.309，其中美国-津巴布韦（0.548）和美国-日本（0.459）的遗传分化最高。

3. 连锁不平衡（LD）分析

基于配对r²估计的LD衰减分析显示，r²在约60 kb处下降至0.20，在200 kb以远接近基线水平（<0.10）。这一LD衰减距离被用于后续候选基因挖掘的窗口定义。

4. 全基因组关联分析（GWAS）

研究人员采用三种方法进行了关联分析：

• •
  FarmCPU：检测到6个与PR显著相关的SNPs，分别位于1、12、13、16和17号染色体，其中1号染色体上的位点（ss14979891）效应最大，可解释24.03%的表型变异。
• •
  随机森林（RF）：基于特征重要性得分，检测到15个SNPs，分布在1、7、9、15、19和20号染色体。
• •
  支持向量回归（SVR）：检测到41个SNPs，集中在2、3、7、9、10、11、12、14、16、17、19和20号染色体。
  RF和SVR识别出了许多FarmCPU未检测到的小效应位点，证明了机器学习方法在捕获非线性关联和更多变异方面的优势。

5. 显著SNPs的等位基因效应

对FarmCPU识别的SNPs进行等位基因效应分析显示，其加性效应范围从-0.138到0.189。基因型特异性表型分布图验证了关键SNPs（如ss14979891和ss14980636）的不同基因型间RI存在显著差异。

6. 候选基因挖掘与功能富集

研究人员在显著SNPs上下游±60 kb的LD窗口内挖掘候选基因：

• •
  FarmCPU：优先选出6个候选基因，包括WRKY转录因子Glyma.01G128100、一个参与基础免疫的MLO样膜蛋白、一个细胞壁生物合成相关基因等。
• •
  RF：优先选出多个激酶相关基因、突触结合蛋白和富含亮氨酸重复序列（LRR）蛋白等。
• •
  SVR：识别出更广泛的基因库，以类受体激酶（RLKs）、核苷酸结合富含亮氨酸重复序列蛋白（NLRs，包括TIR-NBS-LRR免疫受体）以及WRKY、NAC家族转录因子为主。其中Glyma.20G173100编码一个WD40重复蛋白，可能在免疫信号传导的翻译后调控中发挥作用。
  跨方法重叠的基因有限，仅Glyma.01G128100（FarmCPU + RF）和Glyma.09G197000（RF + SVR）。基因本体（GO）富集分析显示，与蛋白磷酸化、蛋白激酶活性和防御反应相关的条目得到富集，但经过Benjamini–Hochberg校正后，仅“转移酶活性”这一分子功能条目保持显著。

四、 讨论与结论

讨论部分总结：

本研究通过整合FarmCPU GWAS与机器学习（RF和SVR）方法，阐明了大豆锈病部分抗性（PR）的遗传结构。这种整合方法比单一模型分析更能改善真实关联的检测并减少假阳性。研究揭示了PR的显著表型变异和数量遗传性质。环境异质性导致了显著的基因型与环境互作，但通过多环境混合模型和BLUEs进行GWAS，有助于鉴定在不同环境下效应稳定的位点。本研究中观察到的中等至高的遗传力估计值与其他研究一致，支持PR主要由遗传因素控制且具有多基因性质。群体结构分析揭示了与成熟期组和起源相关的遗传聚类，这在GWAS中需要严格校正以避免假阳性。本研究中观察到的连锁不平衡衰减（~60 kb，r²=0.2）支持基因水平的关联作图。RF和SVR模型检测到了许多FarmCPU未发现的小效应位点，证明机器学习能够捕获非线性关联，提升对复杂性状变异的解析。候选基因分析揭示了一系列与植物防御相关的基因，包括WRKY转录因子（如Glyma.01G128100/ GmWRKY4）、类受体激酶（如GmLMM1）和WD40结构域蛋白（如Glyma.20G173100），表明PR涉及转录调控、激酶信号传导和泛素介导的翻译后控制等多种通路。未在已知Rpp基因富集的18号染色体上检测到显著关联，这可能反映了本群体中PR（多基因、小效应）而非单基因抗性占主导。研究也承认了其局限性，如依赖自然感染、中等标记密度、使用单一参考基因组等。未来研究可整合泛基因组、高密度测序、转录组分析和受控接种实验来验证关键位点，并利用基因组选择（GS）和标记辅助聚合育种来培育具有持久广谱抗性的品种。

结论部分翻译：

本研究表明，田间条件下对大豆锈病的部分抗性（PR）由多基因结构控制，该结构可通过整合多基因座统计GWAS（FarmCPU）与基于机器学习的关联分析（RF和SVR）得到有效解析。这个组合分析框架一致地识别出了与数量抗性相关的基因组区域，支持了跨方法检测位点的稳健性。

关键的抗性相关区域，包括含有Glyma.01G128100（GmWRKY4）、GmLMM1和Glyma.20G173100的区域，表明转录调控、类受体激酶信号传导和泛素介导的翻译后控制是应对Phakopsora pachyrhizi的部分抗性（PR）反应的核心组成部分。这些结果表明，持久的田间抗性源于协调的调控和信号通路，而非单个主效抗性基因。

通过将统计学支持的基因组区域与生物学相关的防御机制联系起来，这项工作为解释大豆中的数量抗性提供了机制基础。本文鉴定的位点代表了下游功能验证的优先候选区域，并为旨在稳定跨环境抗性的标记辅助育种策略提供了一个基因组框架。