【无标题】

https://academic.oup.com/genetics/article/221/4/iyac091/6603118?login=true#367435093

摘要

生育酚和生育三烯酚（统称为维生素E）是脂溶性抗氧化剂，对植物适应性和人类健康都非常重要。种子油是维生素E的主要膳食来源，但这些油通常富含维生素E活性较低的生育酚异构体。生育酚和生育三烯酚的生物合成途径在植物物种中是保守的，但对于大多数谷物种子中生育酚和生育三烯酚水平自然变异的基因和机制的综合了解仍然有限。为了解决这一问题，我们利用了包含约1500个自交系的玉米Ames种质资源的高映射分辨率，这些自交系被鉴定出1220万个单核苷酸多态性，从而生成了用于遗传定位的代谢组学（成熟谷粒中的生育酚和生育三烯酚）和转录组学（发育中的谷粒）数据集。通过结合全基因组和全转录组关联研究的结果，我们共鉴定出13个候选因果基因位点，其中5个之前未与玉米谷粒中的生育酚和生育三烯酚相关联：4个生物合成基因（arodeH2同源基因、dxs1、vte5和vte7）以及一个质体S-腺苷甲硫氨酸转运蛋白（samt1）。对这13个基因位点进行表达数量性状位点（eQTL）定位表明，它们主要受顺式eQTL调控。通过联合统计分析，我们推测顺式作用变异是共定位的eQTL和GWAS关联信号的原因。我们的多组学方法提高了统计功效和定位分辨率，从而能够详细描述玉米谷粒中生育酚和生育三烯酚积累的遗传和调控结构，并为正在进行的生物强化工作提供了见解，以培育和/或设计玉米和其他谷物中的维生素E和抗氧化水平。

引言

生育酚和生育三烯酚（统称为生育酚和生育三烯酚）是一类植物合成的脂溶性抗氧化剂，它们具有源自对羟基苯丙酮酸（HGA）的色醇环和源自异戊二烯的疏水侧链。生育酚的饱和侧链源自植醇二磷酸（PDP），而生育三烯酚的侧链含有3个双键，源自牻牛儿牻牛儿二磷酸（GGDP）。生育酚和生育三烯酚以4种生物合成相关的异构体形式存在（α、β、δ和γ），它们的色醇环上甲基基团的数量和位置有所不同。在这些生育酚和生育三烯酚中，α-生育酚具有最高的维生素E活性（DellaPenna和Mène-Saffrané，2011），而生育三烯酚则具有更强的抗氧化活性（Sen等人，2006）。尽管在人类群体中，导致共济失调和肌病的严重维生素E缺乏症较为罕见（Traber，2012），但某些人群的饮食中维生素E摄入量不足（Ford等人，2006；McBurney等人，2015），并且与心血管疾病风险增加有关（Knekt等人，1994；Kushi等人，1996）。在植物种子中，生育酚和生育三烯酚含量较高，它们在种子储存和萌发过程中提供对脂质过氧化的保护（Sattler等人，2004）。然而，α-生育酚并非大多数谷物种子油中的主要生育酚和生育三烯酚，这限制了人类和动物的饮食维生素E摄入量（DellaPenna和Mène-Saffrané，2011）。

生育酚和生育三烯酚仅由光合生物合成，其生物合成途径在植物界高度保守（DellaPenna和Mène-Saffrané，2011）。在生育酚合成的决定性步骤中（见图1），对羟基苯丙酮酸植醇转移酶（VTE2）将来自莽草酸途径的PDP和HGA缩合，生成2-甲基-6-植基-1,4-苯醌（MPBQ）（Sattler等人，2004）。在单子叶植物谱系中，HGA也可以与GGDP通过对羟基苯丙酮酸牻牛儿牻牛儿转移酶（HGGT1）缩合，生成2-甲基-6-牻牛儿牻牛儿基-1,4-苯醌（MGGBQ），这是生育三烯酚合成的决定性步骤。MPBQ和MGGBQ是MPBQ/MGGBQ甲基转移酶（VTE3）和γ-生育酚甲基转移酶（VTE4）进行一系列甲基化反应以及生育酚环化酶（VTE1）进行环化反应的底物，这些反应的顺序和数量生成了生育酚和生育三烯酚的α、β、δ和γ异构体（Shintani和DellaPenna，1998；Porfirova等人，2002；Cheng等人，2003；Van Eenennaam等人，2003；Sattler等人，2004）。尽管用于生育三烯酚合成的GGDP直接来自异戊二烯途径，但用于生育酚合成的PDP的生成过程更为复杂，尚未完全明确。尽管单子叶植物和双子叶植物在叶片生育酚合成方面存在差异，但种子中生育酚的合成需要叶绿素生物合成（Diepenbrock等人，2017；Zhan等人，2019）以及VTE7的活性，VTE7是一种α/β水解酶，它与叶绿素合成相互作用以释放植醇（Albert等人，2022）。然后，植醇通过植醇激酶（VTE5）和植醇磷酸激酶（VTE6）的作用依次磷酸化为PDP（Valentin等人，2006；Vom Dorp等人，2015）。

在过去的十年中，通过在定位群体中进行全基因组关联研究（GWAS），已经鉴定出几个与玉米籽粒中生育酚和生育三烯酚含量和组成自然变异相关的位点。一些研究报道了VTE4基因与玉米籽粒中α-生育酚浓度之间存在强关联（Li等人，2012；Lipka等人，2013；Wang等人，2018；Baseggio等人，2019），而对于VTE1、HGGT1以及一个芳氨酸/前磷酸脱水酶与玉米籽粒中生育三烯酚水平之间的关联则相对较弱（Lipka等人，2013；Baseggio等人，2019）。Wang等人（2018）指出，核心生育酚和生育三烯酚生物合成途径之外的基因在玉米籽粒中生育酚水平方面也发挥着作用，包括参与脂肪酸生物合成、叶绿素代谢和叶绿体功能的基因。尽管这些综合研究为玉米籽粒中生育酚和生育三烯酚水平的遗传基础提供了一些见解，但这些研究受到定位群体规模和标记密度的限制。

通过在美国玉米嵌套关联作图（NAM）5000系群体中进行联合连锁（JL）分析和GWAS，Diepenbrock等人（2017）鉴定出50个独特的QTL与生育酚和生育三烯酚籽粒性状相关。其中，13个QTL被解析到7个已知途径基因（dxs2、sds、arodeH2、hppd1、hggt1、vte3和vte4）和6个非途径基因（por1、por2、snare、ltp、phd和fbn），这些非途径基因编码的预测功能之前未与生育酚和生育三烯酚性状相关联。尽管玉米籽粒是一种非绿色、非光合组织，但发现2个原叶绿素还原酶（POR1和POR2）是控制总生育酚的主要位点，并被假设是参与一个循环的一部分，该循环为玉米胚（其中含有极低但可检测到的叶绿素水平）中的生育酚合成提供叶绿素衍生的植醇。通过转基因过表达，Zhan等人（2019）证实了por2在玉米籽粒中生育酚积累中的参与。尽管美国NAM群体为这13个QTL提供了无与伦比的定位分辨率，但由于NAM群体中存在中等效应或有限对比基因型的37个QTL未能被解析到基因水平，因此结合更多样化的定位群体的其他定位方法在鉴定更多潜在候选因果基因方面可能具有显著的实用性。

材料与方法

遗传定位的实验设计

2015年和2017年，我们在爱荷华州立大学的Ames校区种植了来自Ames种质资源库（Romay等人，2013）的1815个玉米自交系，每个自交系种植一个重复。Ames种质资源库采用增广完全区组设计。对于每年的种植，设置了两个方向的区组：每个范围区组由3行相邻的小区组成，每个通道区组由8列相邻的小区组成。每个通道和范围区组内至少种植了一个B73对照小区。根据Romay等人（2013）记录的抽丝天数（开花时间），2015年的自交系被分为2个层级；根据2015年的授粉日期，2017年的自交系被分为3个层级。实验单元为1行3.05米的小区，大约有18株植株，行距为0.76米，走道宽0.76米。每小区大约有6株植株进行自交授粉，并在生理成熟时手工收获。所有干燥并脱粒的果穗的籽粒被混合，形成一个有代表性的样本。对于每个小区，25粒籽粒在IKA管式研磨机（IKA-Werke，德国斯陶芬）中研磨，研磨后的组织样本被保存在-80°C的冻存管中。

表型数据分析

从15到20毫克的研磨籽粒中提取生育酚和生育三烯酚，并通过高效液相色谱（HPLC）和荧光法进行定量，方法如之前所述（Lipka等人，2013）。简而言之，对来自1762个自交系和重复的B73对照小区的3539份籽粒样本中的生育酚和生育三烯酚进行了评估。评估的9种生育酚和生育三烯酚表型（以干种子每克微克数表示）如下：α-生育酚（αT）、δ-生育酚（δT）、γ-生育酚（γT）、α-生育三烯酚（αT3）、δ-生育三烯酚（δT3）、γ-生育三烯酚（γT3）、总生育酚（ΣT，计算为αT+δT+γT）、总生育三烯酚（ΣT3，计算为αT3+δT3+γT3）和总生育酚和生育三烯酚（ΣTT3，计算为ΣT+ΣT3）。从原始HPLC数据中识别并过滤出统计学异常值，随后进行混合线性模型分析，该模型模拟了遗传和非遗传（田间和实验室）效应，以产生1762个自交系的最佳线性无偏估计（BLUE）值（补充表1）和如补充方法中所述的遗传力估计值。考虑到极端形态的籽粒类型可能基于干样品重量具有夸大的生育酚和生育三烯酚浓度，我们谨慎地排除了265个被分类为甜玉米、爆米花玉米或其他具有胚乳突变的自交系。

基因型数据

我们使用了Wu等人（2021）在补充方法中描述的基因型数据处理和BEAGLE v5.0（Browning等人，2018）插补方法。为了进行GWAS，我们生成了一个包含1462个既有基因型又有表型数据的自交系的插补标记数据集，这些数据集由玉米HapMap 3.2.1（Bukowski等人，2018）中的14,613,169个SNP组成，这些SNP位于B73 RefGen_v4坐标系中。这些SNP位点进一步经过质量过滤，以产生一个高质量的12,184,805个SNP集合，这些SNP的等位基因频率（MAF）≥1%，预测剂量r²（DR²）≥0.80（补充数据集1），用于使用混合线性模型进行GWAS。在PLINK版本1.9（Purcell等人，2007）中，使用100kb的滑动窗口和25个SNP的步长，对12,184,805个SNP的完整集合进行连锁不平衡（LD）修剪，构建了2个减少标记集（补充数据集1）：（1）7,319,895个SNP，其成对r²<0.99，用于使用多位点混合模型（MLMM）进行GWAS；（2）344,469个SNP，其成对r²<0.10，用于估计群体结构和亲缘关系。

全基因组关联研究（GWAS）

我们对1462个自交系的9种生育酚和生育三烯酚表型进行了GWAS，采用的是之前描述过的方法（Wu等人，2021）。简而言之，为纠正异方差性和误差项的非正态性，我们使用Box-Cox幂变换程序（Box和Cox，1964）来选择一个合适的lambda值，用于变换每种表型的非负BLUE值（补充表2）。鉴于在模型拟合过程中产生了几个负的BLUE值，我们在应用变换之前加上了一个常数，使所有值均为正数且不小于1E-09（补充表2）。我们使用混合线性模型（Yu等人，2006），在R包GAPIT版本2018.08.18（Lipka等人，2012）中采用先前确定的群体参数近似（Zhang等人，2010），对12,184,805个SNP与1462个自交系的变换后的BLUE值之间的关联进行了检验（补充表3）。在GAPIT中，使用344,469个SNP的简化集合，通过VanRaden方法I（VanRaden，2008）计算亲缘关系矩阵和主成分（PC）。通过贝叶斯信息准则（BIC）（Schwarz，1978）确定每个表型拟合混合线性模型时包含的主成分的最佳数量。采用基于似然比的R²统计量（R²LR）（Sun等人，2010）来近似估算SNP解释的表型变异量。在基础R版本4.0.2（R Core Team，2018）中，使用“p.adjust”函数对测试SNP的P值应用假发现率（FDR）多重检验校正程序（Benjamini和Hochberg，1995）。

为了控制大效应位点的影响，我们采用了Segura等人（2012）提出的多位点混合模型（MLMM）方法，由Wu等人（2021）实现，使用7,319,895个SNP的简化集合对每种变换后的生育酚和生育三烯酚表型进行GWAS，通过去除完全相关的SNP来减轻模型限制。在每个模型中，我们使用了与混合线性模型GWAS相同的亲缘关系矩阵和通过BIC确定的最佳主成分数。

转录组分析的实验设计

2018年，我们在爱荷华州立大学种植了1815个玉米自交系中的1023个，这些自交系来自Ames种质资源库，用于评估籽粒中的生育酚和生育三烯酚，以及美国玉米嵌套关联作图（NAM）小组的5个额外的创始人（Yu等人，2008；McMullen等人，2009），每个自交系种植一个重复。最初构建的种质资源集合包括256个自交系，这些自交系至少满足以下标准之一：（1）在2015年的田间试验中，籽粒代谢物表型极端（高或低）；（2）美国NAM小组的创始人；或（3）有可用的基因组组装。另外随机选择了771个自交系，以增加遗传多样性和样本量。1028个非对照自交系根据2015年和2017年田间试验中记录的授粉日期被划分并随机分配到24个增广不完全区组中，并分为2个层级。每个区组内都种植了B73对照小区，以控制整个田间的空间变异。此外，在每个区组内随机位置种植了2个本地对照，以考虑跨越一个多月的新鲜收获日期的时序变异。在每个区组内，选择在前几年开花最晚的自交系作为2个本地对照之一；如果存在并列情况，则选择基于过滤后的部分插补GBS数据集（补充数据集2）样本调用率最高的自交系。每个选定的本地对照也种植在其相邻的晚开花区组中，因此在第2至24区组中各有2个本地对照。另外，识别出一个早开花的自交系（S 117）作为本地对照，并种植在第1区组中，确保第1区组中也种植了2个本地对照。此外，25个本地对照（S 117、C38、A508、A641 Goodman-Buckler、C31、807、LH202、764、PHG71、PHB47、SD101、A680、B93、NC292、NC280、LH208、H100、NC252、NC314、LH51、CI 187-2、NC324、Mo11、NC334和M37W）被种植在单独的第三层级中，以考虑这些自交系的田间效应。实验单元为1行小区，其尺寸和植株数量与遗传定位实验中使用的相同。在每个小区的约6个授粉果穗中，手工收获一个自交果穗，大约在授粉后23天（DAP），然后立即将剥去 husk 的果穗浸入液氮中，并用干冰覆盖（部分样本在-80°C下冷冻）直至脱粒。选择23 DAP时间点是为了捕捉生育酚和生育三烯酚的最大增加以及已知途径基因的强表达（Diepenbrock等人，2017）。为了控制时序效应，从第三层级的本地对照中，每天在约23 DAP时手工收获一个自交果穗，所有收获的果穗在脱粒前都经过相同的液氮和干冰处理。每个冷冻果穗的中部在干冰上单独脱粒，其籽粒在-80°C下保存。总共收集了1012个非对照和107个对照籽粒样本。

RNA提取和3’ mRNA测序

每份样本取8到10粒冷冻籽粒，使用液氮冷却的研磨杯在IKA管式研磨机（IKA-Werke，德国斯陶芬）中研磨，取约100毫克研磨后的组织用于RNA提取，采用改良的热硼酸法（Wan和Wilkins，1994）。按照Hershberger等人（2022）的方法，对RNA样本进行DNase处理并检查质量。将RNA样本随机分配到96孔板中，并在干冰上过夜运送到康奈尔生物技术研究所的基因组学实验室。每个板的提交样本中都包含阳性对照，即每个板的4个孔中分别分装相同的B73对照RNA池，以及每个板的4个阴性对照，由水组成。使用Lexogen QuantSeq 3’ mRNA-Seq文库试剂盒FWD（Lexogen，新罕布什尔州格林兰）构建文库，并在Illumina NextSeq 500上进行测序，产生85nt的单端读取（Illumina，圣地亚哥，加利福尼亚），每个板分成两半，每半在单独的泳道上进行测序，以达到理想的覆盖度。

表达丰度测定

使用Cutadapt版本2.3（Martin，2011）对3’ QuantSeq读取进行两轮清洗，以去除Illumina接头、前12个碱基和polyA尾。然后使用HISAT2版本2.1.0（Kim等人，2019）将读取比对到B73 RefGen_v4参考基因组（Jiao等人，2017），参数如下：–min-intronlen 20，–max-intronlen 60,000，–dta-cufflinks和–rna-strandness F。使用SAMTools版本1.9（Li等人，2009）对结果比对进行排序。然后使用HTSeq版本0.11.2（Anders等人，2015）中的htseq-count函数，使用B73版本4.59注释生成计数，参数如下：–format=bam，–order=pos，–stranded=yes，–minaqual=10，–idattr=ID，–type=gene和–mode=union。使用DESeq2的rlog函数（Love等人，2014）对1171个RNA样本（1119个籽粒样本加上52个阳性对照）的计数数据进行标准化。所有在所有样本中标准化计数小于或等于零的基因都被从最终计数矩阵中移除。基于采样问题（例如发霉的籽粒等）、比对率、样本间相关值、基因型确认评估和杂合度水平等几个严格的质控措施被实施，以过滤掉低质量样本，具体方法如补充方法所述，并总结在补充表4中，最终得到741个高质量样本用于后续分析。

表达数据分析

包含664个非对照系的665个样本和25个对照系的76个样本的表达数据集，按照Hershberger等人（2022）的方法，在基因水平上进一步严格过滤统计学异常值，以确保用于统计分析的高质量数据。过滤步骤和指标总结在补充表4中。在664个自交系中，我们排除了104个被分类为甜玉米、爆米花玉米或其他具有胚乳突变的自交系，以及另外15个未在GWAS中分析的自交系。最终数据集包含545个自交系中22136个基因的BLUE表达值。

为了考虑影响表达变异的推断混杂因素，PEER方法（Stegle等人，2012）被单独应用于545个自交系×22136个基因的BLUE表达值矩阵，如之前所述（Hershberger等人，2022）。简而言之，通过在因子相关性诊断图的曲线中找到“肘部”，确定最佳的PEER隐因子数量为11个。减去11个PEER隐因子的贡献，生成BLUE表达值的残差数据集（以下称为PEER值）。使用学生化删除残差（Neter等人，1996）来识别并从PEER值集合中移除显著异常值（Bonferroni校正α=0.05）（补充数据集3）。

全转录组关联研究（TWAS）

我们对545个自交系和22136个基因进行了TWAS，采用混合线性模型方法（Yu等人，2006；Zhang等人，2010）。简而言之，使用R包rrBLUP版本4.6（Endelman，2011）中的“GWAS”函数和将“P3D”选项设置为FALSE，为每个生育酚和生育三烯酚表型（变换后的BLUE值，响应变量）和表达基因（经过异常值筛选的PEER值，解释变量）的组合拟合混合线性模型。为了构建545个自交系的SNP标记集，从完整的14613169个SNP中筛选出12018644个双等位基因SNP（DR²≥0.80；MAF≥1%），并在PLINK版本1.9中使用100kb的滑动窗口和25个SNP的步长，将其修剪至328892个SNP，成对r²<0.10（补充数据集4）。如上所述，从328892个SNP中生成主成分（PC）和亲缘关系矩阵。所有生育酚和生育三烯酚表型的最佳模型包括亲缘关系且不包括PC，这是通过贝叶斯信息准则（BIC）（Schwarz，1978）确定的。鉴于vte7位点的读取数据经过独特处理以考虑串联重复基因（补充方法），因此单独进行了vte7位点的TWAS。通过TWAS检测到的基因被用于使用PANTHER过度表达测试（2022年2月2日发布）（Mi等人，2019）进行GO术语富集分析（生物过程）（GO本体数据库DOI：10.5281/zenodo.6399963，2022年3月22日发布），采用费舍尔精确检验，并在FDR P值<0.05时宣布显著性。只有显著且富集度>1倍的GO生物过程才被考虑。

费舍尔联合检验（Fisher’s Combined Test）

我们按照Kremling等人（2019）的方法，从混合线性模型的GWAS中选取了最小的10%的P值SNP（共1218480个SNP）来进行费舍尔联合检验（FCT）。鉴于GWAS中底部90%的SNP在FCT中不太可能产生新的关联，为了减轻计算负担，我们重点关注顶部10%的SNP。简而言之，基于B73 RefGen_v4组装和B73 v4.59注释，将每个顶部10%的SNP的GWAS P值分配给其最近的基因，然后与该基因的TWAS P值配对。对于在TWAS中未进行测试的基因，其TWAS P值在与GWAS P值合并之前被设置为1。对于每个基因，我们使用R包metap版本1.1（Dewey，2019）中实现的“sumlog”函数进行FCT。

候选基因鉴定

鉴于GWAS、TWAS和FCT在统计功效和独立变量类型方面存在差异，我们没有直接比较不同方法之间的P值阈值，而是按照Kremling等人（2019）的方法，根据P值进行排名。对于每种表型，根据GWAS结果的P值选择顶部0.02%的SNP，选择百分比阈值是根据在美国玉米嵌套关联作图（NAM）小组中检测到的每种生育酚和生育三烯酚籽粒表型的联合连锁（JL）-QTL数量（每种表型12-21个JL-QTL）来确定的（Diepenbrock等人，2017）。考虑到Ames种质资源库中连锁不平衡（LD）的快速衰减（Romay等人，2013），我们按照Wu等人（2021）的方法，在每个位点的峰值SNP周围100kb范围内鉴定候选基因。对于每种表型，根据TWAS和FCT结果的P值选择顶部0.5%的基因，从而使每种方法鉴定出的跨表型独特基因总数与GWAS相当。

我们根据Lipka等人（2013）的生物信息学方法，编制了一个包含126个与籽粒中生育酚和生育三烯酚水平积累相关的先验候选基因的列表（补充表5），以协助鉴定候选基因。通过Vmatch版本2.3.0（Kurtz，2010），将美国NAM小组中与9种生育酚和生育三烯酚籽粒表型相关的50个独特联合连锁QTL共同支持区间（CSI）和GWAS标记的物理位置提升到B73 RefGen_v4坐标系中（补充表6和7），如Wu等人（2021）所述。我们进行了默认参数的BLASTP分析，以鉴定未描述的候选因果基因在拟南芥和水稻中的最佳匹配（补充表8），如之前所述（Wu等人，2021）。

eQTL定位

我们对鉴定出的候选因果基因进行了表达数量性状位点（eQTL）定位。为了进行eQTL定位，我们在R版本4.0.2（R Core Team，2018）中使用GAPIT版本2018.08.18（Lipka等人，2012）实现的混合线性模型，将TWAS方法中筛选出的12018644个SNP分别与每个候选因果基因的PEER值进行关联测试。在TWAS中使用的计算出的主成分（PC）和亲缘关系矩阵也用于eQTL定位，其中最佳主成分数由贝叶斯信息准则（BIC）（Schwarz，1978）确定。为了在复杂的连锁不平衡（LD）模式和强关联信号存在的情况下严格控制I型错误率，我们采用5%的Bonferroni校正显著性阈值（P值≤4.16E-09）来考虑多重检验，峰值SNP的位点鉴定如Wu等人（2021）所述。

变异注释

按照Diepenbrock等人（2021）的方法，进行SNP效应分析（SnpEff；Cingolani等人，2012），以预测位于候选因果基因内的与GWAS相关的SNP的效应。为了定量评估变异位点是否在进化上保守，从两项早期研究（Kistler等人，2018；Ramstein等人，2020）中提取了候选因果基因内相同SNP位点的基因组进化速率分析（GERP；Davydov等人，2010）得分。

eQTL和GWAS因果变异鉴定

为了量化一个变异同时负责GWAS和顺式eQTL信号的概率，我们使用了Hormozdiari等人（2016）提出的eQTL和GWAS因果变异鉴定方法（eCAVIAR）。该方法在整合GWAS和顺式eQTL结果时，考虑了连锁不平衡（LD）模式和等位基因异质性，通过概率模型来实现。我们将eCAVIAR方法应用于通过GWAS检测到的具有显著顺式eQTL信号的候选因果基因位点。对于每对基因-表型，我们将候选因果基因上下游100kb或250kb范围内的所有显著GWAS和eQTL SNP的t值，以及这些SNP在PLINK 1.9（Purcell等人，2007）中计算出的成对LD矩阵作为eCAVIAR软件的输入数据集。只有2个基因（arodeH2 Zm00001d014734和vte1）的峰值eQTL信号距离其相应基因超过100kb，因此仅对这两个基因使用了±250kb的窗口。鉴于等位基因异质性的可能性，将每个位点的最大因果SNP数量设置为3。使用严格的共定位后验概率（CLPP）截止阈值≥0.01来识别在GWAS和eQTL研究中都可能是因果的SNP。

表型变异

我们评估了从玉米Ames种质资源库的两个扩展试验中收获的生理成熟籽粒样本中生育酚和生育三烯酚浓度的定量变异程度。通过高效液相色谱（HPLC）测量的6种生育酚和生育三烯酚化合物显示，γT（约55%）和γT3（约23%）合计占总生育酚和生育三烯酚（ΣTT3）的近80%，而α和δ异构体的生育酚和生育三烯酚各自分别占ΣTT3的约1%（δT3）到10%（αT3）（表1）。具有最高维生素E活性的生育酚和生育三烯酚化合物αT，其平均浓度（5.83微克/克干种子）在所有化合物中排名第三低，仅占ΣTT3的约8%。在化合物类别内，生育酚（r=0.67）和生育三烯酚（r=0.62）的δ和γ异构体之间的BLUE值（α=0.05）的显著成对相关性最强，而化合物类别之间最强的相关性出现在αT与αT3之间（r=0.45）。然而，尽管存在共同的生物合成途径，其他所有化合物对之间的相关性都很弱（-0.15到0.19）（补充图1）。正如从基于系平均值的高遗传力估计（0.77到0.94，表1）所推断的那样，Ames种质资源库中6种生育酚和生育三烯酚化合物以及3种总表型的大部分变异可归因于遗传变异（表1和补充图2）。

表1. 在Ames种质资源库中评估的9种生育酚和生育三烯酚籽粒表型的未变换BLUE值（以微克/克计）的平均值、范围和标准差，以及基于系平均值估计的遗传力及其标准误，跨越2年。

表型a	系数	BLUEs	遗传力
平均值	范围	标准差	估计值
αT	1452	5.83	-1.79到41.36
δT	1456	1.74	-0.33到14.32
γT	1458	42.19	-1.32到158.91
ΣT	1460	49.95	1.79到174.84
αT3	1456	7.87	0.89到23.39
δT3	1454	0.93	0.01到17.05
γT3	1458	17.60	-1.79到90.39
ΣT3	1458	26.55	2.62到111.01
ΣTT3	1460	77.04	18.13到205.36
a
3种总表型的计算方法如下：ΣT=αT+δT+γT；ΣT3=αT3+δT3+γT3；ΣTT3=ΣT+ΣT3。

籽粒中生育酚和生育三烯酚水平的遗传分析

我们通过费舍尔联合检验（FCT）整合了全基因组关联研究（GWAS）和全转录组关联研究（TWAS）的结果，这是一种集成方法，已被证明在检测与玉米籽粒中生育酚和生育三烯酚表型自然变异相关的因果基因方面，比单独的GWAS或TWAS具有更强的统计功效。我们将FCT（排名前0.5%）、GWAS（排名前0.02%）和TWAS（排名前0.5%）的结果与美国玉米嵌套关联作图（NAM）小组中相同籽粒表型的遗传定位结果进行了整合（表2），目的是进一步将NAM小组中之前发现的位点解析到因果基因的水平（图2和补充图3）。在每种分析中，从9种生育酚和生育三烯酚籽粒表型中，GWAS共鉴定出121个位点的720个独特基因，TWAS鉴定出676个基因，FCT鉴定出918个基因（补充表9-12）。其中，330个（GWAS）、299个（TWAS）和646个（FCT）基因位于9种表型的NAM联合连锁（JL）-QTL共同支持区间（CSI）内。在对TWAS中鉴定出的676个基因进行GO术语富集分析时，我们没有发现任何GO生物过程在5%的FDR和>1倍富集水平上显著。

表2. Ames种质资源库中9种生育酚和生育三烯酚籽粒表型的遗传定位结果。

基因ID	位点a	染色体	基因起始	基因结束	GWASb	TWASc	FCTd	NAM JL-QTL CSI Ide
Zm00001d032576	por1	1	231,120,510	231,123,615	ΣTf	δT, γT, ΣT, ΣTT3	δT, γT, ΣT, ΣTT3	5
Zm00001d001896	vte5	2	2,509,567	2,511,414	ΣTT3	–	–	–
Zm00001d006778	vte7g	2	216,443,683	216,448,352	δT	–	–	12
Zm00001d006779		216,461,133	216,468,803
Zm00001d013937	por2	5	25,431,430	25,434,346	δT, γT, ΣT, ΣTT3	αT, δT, γT, ΣT, ΣTT3	αT, δT, γT, ΣT, ΣTT3	24
Zm00001d014734	arodeH2	5	61,099,110	61,100,192	γT3, ΣT3	–	–	25
Zm00001d014737	arodeH2	5	61,117,986	61,119,432	γT3, ΣT3	–	δT3, γT3, ΣT3	25
Zm00001d015356	hppd1	5	86,084,655	86,086,755	γT3, ΣT3	–	δT3, γT3, ΣT3, ΣTT3	26
Zm00001d015985	vte1	5	136,805,708	136,822,194	δT3	–	δT, δT3	26
Zm00001d017746	vte4	5	205,825,586	205,829,216	αT, αT3, δT, γT, γT3, ΣTT3	αT, αT3, γT3	αT, αT3, δT, γT, γT3	28
Zm00001d017937	samt1	5	210,385,310	210,401,948	αT3, δT	δT	δT, δT3	29
Zm00001d038170	dxs1	6	150,418,144	150,422,431	–	γT3	–	32
Zm00001d019060	dxs2	7	14,494,700	14,497,925	–	δT3, γT3, ΣT3	–	35
Zm00001d046558	hggt1	9	95,895,575	95,899,061	αT3, δT3, γT3	δT3, γT3, ΣT3	αT3, δT3, γT3, ΣT3, ΣTT3	45
a
玉米基因位点的名称和符号用小写、斜体字符表示，以符合玉米遗传学命名法。
b
全基因组关联研究。
c
全转录组关联研究。
d
费舍尔联合检验。
e
玉米NAM小组中分析的9种生育酚和生育三烯酚籽粒表型的联合连锁-数量性状位点（JL-QTL）结果的共同支持区间（CSI）（补充表6），其中包含基因的开放阅读框。
f
3种总表型的计算方法如下：ΣT=αT+δT+γT；ΣT3=αT3+δT3+γT3；ΣTT3=ΣT+ΣT3。
g
由于缺乏围绕重复的α/β水解酶基因（Zm00001d006778和Zm00001d006779）的非同源SNP，vte7位点无法在FCT中进行测试。

在Diepenbrock等人（2017）通过美国玉米嵌套关联作图（NAM）小组鉴定的与籽粒中生育酚和生育三烯酚相关的13个基因位点中，有5个（por1、por2、vte4、hggt1和hppd1），这些位点在NAM小组中倾向于产生大效应，通过FCT在Ames种质资源库中检测到与一种或多种表型相关（表2）。在这5个基因中，por1、por2、vte4和hggt1也被GWAS和TWAS鉴定出来，而hppd1仅通过GWAS检测到。相比之下，dxs2是Diepenbrock等人（2017）鉴定的另外13个基因之一，仅通过TWAS检测到。aro-deH2的两个拷贝（Zm00001d014734和Zm00001d014737）位于γT3和ΣT3的GWAS峰值SNP的100kb范围内，其中Zm00001d014734之前被Diepenbrock等人（2017）鉴定为一个与αT3和ΣT3的遗传控制相关的具有小等位基因效应的基因。然而，在Ames种质资源库中，Zm00001d014737通过FCT和GWAS被检测到，而Zm00001d014734仅通过GWAS检测到。总的来说，我们重新鉴定出了Diepenbrock等人（2017）中的13个基因中的7个，并暗示aro-deH2的第二个拷贝控制生育三烯酚的变异。

这8个基因的检测，其中一个是新的关联，展示了我们综合遗传定位方法的基因水平分辨率；因此，它被应用于更好地解析NAM联合连锁（JL）-QTL共同支持区间（CSI），并在Ames种质资源库中检测新的位点。总共，4个NAM JL-QTL CSI被更精细地剖析，导致与3个位点（samt1、vte7和dxs1）的新关联，以及第四个（vte1）的更精确定位，这在NAM小组中并未完全解析。编码SAM转运蛋白的基因（samt1，Zm00001d017937）被FCT、GWAS和TWAS检测到。Zm00001d017937编码一种与拟南芥S-腺苷甲硫氨酸转运蛋白1（SAMT1，AT4G39460）有77%同源性的蛋白（补充表8），该蛋白将SAM（生育酚和生育三烯酚的辅底物，用于VTE3和VTE4甲基转移酶）通过质体外膜运输，并在沉默（弯叶烟草）或敲除（拟南芥）时对叶片生育酚水平产生负面影响（Bouvier等人，2006；Palmieri等人，2006）。在考虑与δT最相关的GWAS中排名前0.02%的SNP时，由串联重复基因（Zm00001d006778和Zm00001d006779）组成的vte7位点被发现在一个相关SNP的大约64kb处（图2）。这个相同的SNP是与δT相关的位点的峰值，该位点由45个在5% FDR下显著的SNP组成（补充表12），与NAM小组相比，这为在Ames小组中检测到vte7提供了更强的证据。除了更精细地剖析NAM JL-QTL CSI外，我们还通过GWAS单独检测到vte5与ΣTT3的显著关联，这是首次报道该位点与玉米中任何生育酚和生育三烯酚籽粒性状的关联。因此，Ames小组不仅为现有的NAM JL-QTL CSI提供了基因水平的分辨率，还使得在NAM小组中未检测到的位点的鉴定成为可能。

我们还采用MLMM（多位点混合线性模型）方法进行了GWAS分析，这使我们能够更好地解析由大效应位点支撑的关联信号。在通过混合线性模型的GWAS检测到的11个基因位点中，有8个基因（por2、vte1、两个aro-deH2拷贝、hppd1、vte4、samt1和hggt1）位于至少1个MLMM选择的SNP的100kb范围内，这些模型解释了3%到37%的表型变异。虽然定位精度稍低，但por1基因距离多个MLMM选择的SNP之一约162kb，这些SNP与生育酚表型相关。尽管vte5和vte7在混合线性模型的GWAS中被检测到，但它们都没有通过MLMM方法被检测到。在距离vte4 100kb范围内的MLMM选择的SNP中，分别有2到4个SNP被选为与αT、αT3和γT相关，而每个δT和γT3只选了1个SNP。相比之下，MLMM从包含hggt1的约1.2Mb基因组区域中为δT3、γT3和ΣT3选择了2到3个SNP，但只有2个MLMM选择的SNP位于hggt1的100kb范围内。正如之前在Goodman-Buckler小组中对vte4的假设，MLMM选择多个独立SNP表明，在Ames小组中，vte4和hggt1位点存在多个因果变异（即等位基因异质性）。

我们采用了两种方法来识别可能改变编码氨基酸的候选因果变异。SnpEff工具预测了位于候选因果基因内的320个GWAS相关SNP的功能效应，结果为por1、por2、vte1、vte4、samt1和hggt1预测出15个错义变异注释，这些变异被认为具有中等效应。在15个错义变异中，有9个来自两项早期研究（Kistler等人，2018；Ramstein等人，2020）中的一个或两个的GERP评分可用。这9个位点中，每个位点至少有1个GERP评分大于0，表明这些编码区内的位点在进化上受到约束。por1、vte1、samt1和hggt1中的6个错义变异具有最高的正GERP评分（>2），表明它们可能更具破坏性。鉴于vte1未被TWAS检测到，并且有3个可能是更具破坏性的变异，这些氨基酸变化可能具有较高的功能重要性，因为它们可能影响生育酚环化酶的活性。然而，需要通过实验评估VTE1的同工酶的酶活性水平，以确定这些或其他变异是否为补偿性突变。