结合 GWAS 和 TWAS 鉴定玉米籽粒中生育色醇水平的候选致病基因

摘要

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

关键词

GWAS；TWAS；eQTL；代谢物；玉米

引言

生育酚（tocochromanols）包括生物合成相关的生育酚（tocopherols）和生育三烯酚（tocotrienols），是一类由植物合成的脂溶性抗氧化剂。它们含有一个源自对羟基苯甲酸（homogentisic acid, HGA）的色醇环，以及由异戊二烯衍生的疏水侧链。生育酚的饱和侧链源自植醇二磷酸（phytyl-diphosphate, PDP），而生育三烯酚的侧链含有3个双键，源自牻牛儿基牻牛儿基二磷酸（geranylgeranyl diphosphate, 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）。

翻译

生育酚（Tocochromanols）仅由光合生物合成，其生物合成途径在植物界中高度保守（DellaPenna 和 Méne-Saffrané, 2011）。在生育酚合成的决定性步骤中（见图1），对羟基苯甲酸植醇转移酶（VTE2）将植醇二磷酸（PDP）和来自莽草酸途径的对羟基苯甲酸（HGA）缩合，生成2-甲基-6-植基-1,4-苯醌（MPBQ）（Sattler 等, 2004）。在单子叶植物谱系中，对羟基苯甲酸也可以与牻牛儿基牻牛儿基二磷酸（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）指出，核心生育酚生物合成途径之外的基因在调控玉米籽粒生育酚水平中发挥作用，包括参与脂肪酸生物合成、叶绿素代谢和叶绿体功能的基因。

虽然这些综合研究为理解玉米籽粒中生育酚水平的遗传基础提供了一些见解，但这些研究受到定位面板大小和标记密度的限制。通过在包含 5000 行的美国玉米嵌套作关联定位（NAM）面板中进行联合连锁（JL）分析和 GWAS，具有高统计能力的 Diepenbrock 等（2017）识别了 50 个独特的与生育酚籽粒性状相关的数量性状位点（QTL）。其中，13 个 QTL 被解析为 7 个先前途径基因（dxs2、sds、arodeH2、hppd1、hggt1、vte3 和 vte4）以及 6 个非途径基因（por1、por2、snare、ltp、phd 和 fbn），这些基因编码的功能以前并未与生育酚性状相关联。

尽管玉米籽粒是一种非绿色、非光合组织，但发现两个原叶绿素还原酶（POR1 和 POR2）是控制总生育酚的主要位点，并假设它们是为玉米胚胎中生育酚合成提供叶绿素衍生的植醇的循环的一部分，玉米胚胎中含有极低但可检测到的叶绿素水平。por2 在玉米籽粒中生育酚积累中的作用通过转基因确认（Zhan 等, 2019）。尽管美国 NAM 面板为这 13 个 QTL 提供了无与伦比的定位分辨率，但在 NAM 面板中对生育酚具有中等效应或有限对比基因型的 37 个 QTL 无法解析到基因水平，因此，与其他更多样化的定位面板结合的其他定位方法可能在识别额外的潜在因果基因方面具有重要用途。

材料与方法

遗传图谱的实验设计

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

表型数据分析

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

基因型数据

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

全基因组关联研究

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

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

转录组分析的实验设计

2018年，我们在爱荷华州立大学种植了1023份来自阿姆斯种群（Ames panel）的玉米自交系，这些自交系曾被用于评估籽粒中的生育烷醇含量，同时还有5份美国玉米巢式关联作图（NAM）种群（Yu等，2008；McMullen等，2009）的额外亲本。这些种质最初通过以下标准筛选出256份自交系：（1）在2015年的田间试验中，籽粒代谢物表型表现出极端性（高或低）；（2）美国NAM种群的亲本；或（3）具有可用的基因组组装。此外，随机选择了771份自交系以增加遗传多样性和样本量。1028份非对照自交系根据2015年和2017年田间试验记录的授粉日期被划分并随机分配到24个增广不完全区组中，并分为2个层级。每个区组内都种植了B73对照小区，以控制田间空间变异。此外，每个区组内还随机种植了2份本地对照，以考虑跨越一个多月的新鲜收获日期的时间变异。在每个区组内，选择在前几年开花最晚的自交系作为2份本地对照之一；如果存在并列情况，则选择基于过滤后的部分填补的GBS数据集（补充数据集2）样本调用率最高的自交系。每份选定的本地对照还被种植在其相邻的晚花区组中，因此在2号至24号区组中各有2份本地对照。此外，还确定了一条早花自交系（S117）作为本地对照，并种植在1号区组中，以确保1号区组中也有2份本地对照。此外，25份本地对照（S117、C38、A508、A641、Goodman-Buckler、C31、807、LH202、764、PHG71、PHB47、SD101、A680、B93、NC292、NC280、LH208、H100、NC252、NC314、LH51、CI187-2、NC324、Mo11、NC334和M37W）被种植在单独的第三层级中，以考虑这些自交系的田间效应。

实验单元为1行的小区，其尺寸和植株数量与遗传图谱实验中使用的相同。每个小区的6个授粉果穗中，有一个自交果穗在授粉后23天（DAP）手工收获，随后立即将去皮的果穗浸入液氮中，并用干冰覆盖（部分样品在-80℃下冷冻）直至脱粒。选择23 DAP这一时间点是为了捕捉生育烷醇的最大增加量以及已知途径基因的强烈表达（Diepenbrock等，2017）。为了控制时间效应，每天在新鲜收获时，从第三层级的本地对照中手工收获一个自交果穗，并在脱粒前对所有收获的果穗进行相同的液氮和干冰处理。每个冷冻果穗的中部单独在干冰上脱粒，其籽粒保存在-80℃。总共收集了1012份非对照和107份对照籽粒样本。

RNA提取与3’ mRNA测序

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

表达丰度的确定

使用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份被分类为甜玉米、爆裂玉米或其他胚乳突变体的自交系，以及未在全基因组关联分析（GWAS）中分析的另外15份自交系。最终数据集包含545份自交系中22,136个基因的最佳线性无偏估计（BLUE）表达值。

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

转录组全关联研究（TWAS）

我们对545份自交系和22,136个基因进行了转录组全关联研究（TWAS），采用混合线性模型方法（Yu等，2006；Zhang等，2010）。简而言之，我们使用R软件包rrBLUP版本4.6（Endelman，2011）中的“GWAS”函数，并将“P3D”选项设置为FALSE，拟合了每种生育烷醇表型（转换后的BLUE值，作为响应变量）与表达基因（经过异常值筛选的PEER值，作为解释变量）的组合的混合线性模型。

为了构建545份自交系的SNP标记集，我们从完整的14,613,169个SNPs中筛选出12,018,644个双等位基因SNPs（DR²≥0.80；MAF≥1%），并使用PLINK版本1.9（补充数据集4）以100kb的滑动窗口和25个SNP的步长进行连锁不平衡修剪，最终得到328,892个成对r² < 0.10的SNPs。基于这328,892个SNPs，我们按照上述方法生成了主成分（PCs）和亲缘关系矩阵。根据贝叶斯信息准则（BIC）（Schwarz，1978）确定的最佳模型包括亲缘关系矩阵，但不包括主成分，用于所有生育烷醇表型的分析。鉴于vte7位点的测序数据经过了独特处理以考虑串联重复基因（补充方法），我们单独对该位点进行了TWAS分析。

通过TWAS检测到的基因被用于基因本体（GO）术语富集分析，分析的领域为生物过程（GO本体数据库DOI: 10.5281/zenodo.6399963，发布日期为2022年3月22日）。我们使用PANTHER过表达检验（发布日期为2022年2月2日）（Mi等，2019），采用Fisher精确检验，并在假发现率（FDR）P值 < 0.05时宣布显著性。只有显著且富集倍数大于1的GO生物过程才被考虑。

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

我们根据Kremling等（2019）的方法，从混合线性模型的全基因组关联分析（GWAS）中选择最小的10%的P值对应的SNPs（共1,218,480个SNPs）来进行费舍尔联合检验（FCT）。选择前10%的SNPs是为了减少计算负担，因为GWAS中排名靠后的90%的SNPs不太可能在FCT中产生新的关联。

简而言之，我们将每个排名前10%的SNP的GWAS P值分配给其最近的基因，依据是B73 RefGen_v4基因组组装和B73 v4.59注释信息。然后，将这些GWAS P值与相应基因的转录组全关联研究（TWAS）P值进行配对。对于未在TWAS中测试的基因，其TWAS P值被设置为1，然后再与GWAS P值合并。对于每个基因，我们使用R软件包metap版本1.1（Dewey，2019）中实现的“sumlog”函数进行费舍尔联合检验。

候选基因鉴定

由于GWAS（全基因组关联分析）、TWAS（转录组全关联研究）和FCT（费舍尔联合检验）在统计功效和独立变量类型上存在差异，我们没有直接比较不同方法之间的P值阈值，而是参考Kremling等（2019）的方法，根据P值的排名进行选择。对于每种表型，根据GWAS结果的P值选择排名前0.02%的SNPs，选择百分比阈值的依据是美国玉米巢式关联作图（NAM）种群中每种生育烷醇籽粒表型检测到的JL-QTL（联合连锁分析数量性状位点）数量（每表型12-21个JL-QTL）（Diepenbrock等，2017）。考虑到阿姆斯种群中连锁不平衡（LD）的快速衰减（Romay等，2013），我们按照Wu等（2021）的方法，在每个位点的峰值SNP周围100kb范围内鉴定候选基因。对于每种表型，根据TWAS和FCT结果的P值选择排名前0.5%的基因，使得每种方法鉴定出的独特基因总数与GWAS相当。

鉴定候选基因的过程中，我们参考了一份包含126个已知与籽粒生育烷醇积累相关的候选基因列表（补充表5），该列表是按照Lipka等（2013）的生物信息学方法整理的。美国NAM种群中与9种生育烷醇籽粒表型相关的50个独特JL-QTL共同支持区间（CSIs）和GWAS标记的物理位置，通过Vmatch版本2.3.0（Kurtz，2010）转换为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方法中筛选出的12,018,644个SNPs分别与每个候选因果基因的PEER值进行关联分析。

在eQTL定位中，我们使用了TWAS中计算的主成分（PCs）和亲缘关系矩阵，并根据贝叶斯信息准则（BIC）（Schwarz，1978）确定最佳主成分数。为了在复杂的连锁不平衡（LD）模式和强关联信号存在的情况下严格控制I型错误率，我们采用Bonferroni校正后的5%显著性阈值（P值=4.16×10⁻⁹）来校正多重检验，峰值SNP的确定方法如Wu等（2021）所述。

变异注释

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

eQTL和GWAS因果变异鉴定

为了量化一个变异位点同时负责GWAS信号和顺式（cis）-eQTL信号的概率，我们使用了Hormozdiari等（2016）提出的eQTL和GWAS因果变异鉴定方法（eCAVIAR）。该方法在整合GWAS和顺式（cis）-eQTL结果时，考虑了连锁不平衡（LD）模式和等位基因异质性。我们将eCAVIAR方法应用于通过GWAS检测到的具有显著顺式（cis）-eQTL信号的候选因果基因位点。

对于每对基因-表型组合，我们将候选因果基因周围100kb或250kb范围内的所有显著GWAS和eQTL SNPs的t值，以及这些SNPs的成对LD矩阵（使用PLINK 1.9计算，Purcell等，2007）作为eCAVIAR软件的输入数据集。对于两个基因（aroDeH2Zm00001d014734和vte1），其峰值eQTL信号距离基因本身超过100kb，因此使用了6250kb的窗口。考虑到等位基因异质性的可能性，我们将每个位点的最大因果SNPs数量设置为3。我们使用了严格的共定位后验概率（CLPP）阈值0.01来识别在GWAS和eQTL研究中都可能是因果的SNPs。

结果

表型变异

我们评估了阿姆斯种群玉米在两个生长季收获的生理成熟籽粒样本中生育烷醇浓度的定量变异程度。通过高效液相色谱（HPLC）测定的6种生育烷醇化合物显示，γ-生育酚（γT，占55%）和γ-生育三烯酚（γT3，占23%）合计占总生育烷醇（RTT3）的近80%，而α、δ同分异构体的生育酚和生育三烯酚分别仅占RTT3的1%（δT3）到10%（αT3）（表1）。具有最高维生素E活性的α-生育酚（αT）的平均浓度最低，仅为第三低（5.83 µg/g干种子），仅占RTT3的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），阿姆斯种群中每种6种生育烷醇化合物和3种总表型的变异主要归因于遗传变异（表1和补充图2）。

表1的内容是对阿姆斯种群中评估的9种生育烷醇籽粒表型的平均值、范围、标准差（未转换的BLUE值，单位为µg/g）以及基于自交系平均值估计的遗传力及其标准误进行了总结。以下是表格内容的详细解读：

表1：阿姆斯种群中9种生育烷醇籽粒表型的平均值、范围、标准差以及遗传力估计值

表型	自交系数量	BLUE值范围 (µg/g)	平均值 (µg/g)	标准差 (µg/g)	遗传力估计值	遗传力标准误
αT	1452	1.79 - 41.36	5.83	4.59	0.87	0.006
δT	1456	0.33 - 14.32	1.74	1.62	0.85	0.007
γT	1458	1.32 - 158.91	42.19	21.34	0.86	0.006
RT	1460	1.79 - 174.84	49.95	22.65	0.85	0.007
αT3	1456	0.89 - 23.39	7.87	3.21	0.77	0.010
δT3	1454	0.01 - 17.05	0.93	0.97	0.94	0.003
γT3	1458	1.79 - 90.39	17.60	11.71	0.93	0.003
RT3	1458	2.62 - 111.01	26.55	13.32	0.91	0.004
RTT3	1460	18.13 - 205.36	77.04	28.08	0.87	0.006

籽粒生育烷醇水平的遗传分析

我们通过费舍尔联合检验（FCT）整合了GWAS（全基因组关联分析）和TWAS（转录组全关联研究）的结果。FCT是一种集成方法，已被证明比单独的GWAS或TWAS具有更强的统计功效，用于检测与玉米籽粒生育烷醇表型自然变异相关的因果基因（Kremling等，2019）。我们将FCT（排名前0.5%）、GWAS（排名前0.02%）和TWAS（排名前0.5%）的结果与美国玉米巢式关联作图（NAM）种群中相同籽粒表型的遗传定位结果（表2）进行了整合，旨在进一步将NAM种群中发现的位点解析到因果基因水平（图2和补充图3）。

在每种分析中，我们共鉴定出720个独特基因，涉及GWAS中的121个位点、TWAS中的676个基因以及FCT中的918个基因，涵盖了9种生育烷醇籽粒表型（补充表9-12）。其中，330个（GWAS）、299个（TWAS）和646个（FCT）基因位于NAM种群中9种表型的联合连锁分析数量性状位点（JL-QTL）共同支持区间（CSIs）内（Diepenbrock等，2017）。在对TWAS鉴定出的676个基因进行基因本体（GO）术语富集分析时，我们没有发现任何GO生物过程在5%的假发现率（FDR）和>1倍富集水平上显著。

在Diepenbrock等（2017）通过美国NAM种群鉴定出的与籽粒生育烷醇相关的13个基因位点中，有5个（por1、por2、vte4、hggt1和hppd1）在NAM种群中倾向于表现为大效应位点，这些基因通过FCT在阿姆斯种群中被检测到与一种或多种表型相关（表2）。在这5个基因中，por1、por2、vte4和hggt1同时被GWAS和TWAS检测到，而hppd1仅被GWAS检测到。相比之下，dxs2是Diepenbrock等（2017）鉴定出的另外13个基因之一，仅被TWAS检测到。aroDeH2的两个拷贝（Zm00001d014734和Zm00001d014737）位于cT3和RT3的GWAS峰值SNP 100kb范围内，其中Zm00001d014734曾被Diepenbrock等（2017）鉴定为与aT3和RT3遗传控制相关的小效应基因。然而，在阿姆斯种群中，Zm00001d014737被FCT和GWAS检测到，而Zm00001d014734仅被GWAS检测到。总体而言，我们重新鉴定出Diepenbrock等（2017）研究中的13个基因中的7个，并暗示aroDeH2的第二个拷贝可能参与调控生育三烯酚的变异。

这8个基因的检测，其中包括一个新的关联，展示了我们整合遗传定位方法在基因水平上的解析能力；因此，该方法被应用于更好地解析NAM JL-QTL CSI，并在阿姆斯种群中检测到新的位点。总体而言，4个NAM JL-QTL CSI被更精细地剖析，导致与3个位点（samt1、vte7和dxs1）的新关联，以及在尚未完全解析的美国NAM种群中更精确地定位了第四个位点（vte1）。编码S-腺苷甲硫氨酸转运蛋白的基因（samt1，Zm00001d017937）被FCT、GWAS和TWAS检测到。Zm00001d017937编码的蛋白与拟南芥S-腺苷甲硫氨酸转运蛋白1（SAMT1，AT4G39460）具有77%的同源性（补充表8），它将S-腺苷甲硫氨酸（SAM，生育烷醇的甲基供体）运输穿过质体包膜，并在沉默（N. benthamiana）或敲除（拟南芥）时对叶片中的生育酚水平产生负面影响（Bouvier等，2006；Palmieri等，2006）。

在考虑GWAS中与δT最相关的排名前0.02%的SNPs时，发现vte7位点由串联重复基因（Zm00001d006778和Zm00001d006779）组成，距离一个相关SNP 64kb（图2）。这个相同的SNP是与δT相关的位点的峰值，该位点包含45个在5%假发现率（FDR）下显著的SNPs（补充表12），为在阿姆斯种群中检测到vte7提供了比NAM种群更强的证据。另外两个已知途径基因，dxs1（TWAS）和vte1（FCT和GWAS），与一种或多种生育烷醇籽粒表型相关（表2）。除了更精细地剖析NAM JL-QTL CSI外，我们还通过GWAS单独检测到vte5与RTT3的显著关联，这是首次报道该位点与玉米中任何生育烷醇籽粒性状的关联。因此，阿姆斯种群不仅为现有的NAM JL-QTL CSI提供了基因水平的解析，还使得我们能够鉴定出在NAM种群中未被检测到的位点。

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

关联SNPs的功能注释

我们采用了两种方法来识别可能改变编码氨基酸的候选因果基因中的潜在因果变异。SnpEff工具预测了位于候选因果基因内的320个GWAS相关SNPs的功能效应，结果在por1、por2、vte1、vte4、samt1和hggt1中预测到15个错义变异注释，这些变异被认为具有中等效应（补充表14）。在15个错义变异中，有9个在两项早期研究（Kistler等，2018；Ramstein等，2020）中至少有一项提供了GERP评分。这9个位点的GERP评分均大于0，表明这些编码区内的位点在进化上受到约束（Davydov等，2010；Rodgers-Melnick等，2015）。其中，por1、vte1、samt1和hggt1中的6个错义变异具有最高的正GERP评分（>2），表明它们可能更具破坏性（Yang等，2017；Lozano等，2021）。

鉴于vte1未被TWAS检测到，且在6个可能更具破坏性的变异中有3个，这些氨基酸变化可能具有高度的功能重要性，因为它们可能会影响生育酚环化酶的活性。然而，需要通过实验评估VTE1的同工酶的酶活性水平，以确定这些或其他变异是否为补偿性突变。

候选因果位点的eQTL（表达数量性状位点）定位

为了深入了解通过GWAS、TWAS和FCT在阿姆斯种群中鉴定的位点（表2）的调控模式，我们对13个鉴定出的候选因果基因位点分别进行了eQTL（表达数量性状位点）定位分析（图3和补充图4）。在13个位点中，除了1个基因（aroDeH2 Zm00001d014737）外，其余基因均鉴定到了顺式（cis）-eQTL（峰值SNP位于基因1Mb范围内），而总共鉴定到5个反式（trans）-eQTL，涉及4个基因（vte5、por2、dxs1和dxs2）（补充表15）。顺式eQTL的峰值SNP通常位于相应基因的100kb范围内，但aroDeH2 Zm00001d014734（227kb）、dxs2（808kb）和vte1（220kb）是例外（补充表15）。总体而言，顺式eQTL的P值通常小于反式eQTL，但dxs2是一个例外，其反式eQTL的P值小于顺式eQTL。这个dxs2的反式eQTL位于第6染色体上，其峰值SNP距离八氢番茄红素合成酶1（psy1，Zm00001d036345）1.5kb，该基因编码类胡萝卜素生物合成中的第一个且关键的步骤（Hirschberg，2001）。

GWAS和eQTL信号的共定位

通过采用概率方法（eCAVIAR；Hormozdiari等，2016）整合GWAS和eQTL定位结果，我们测试了在具有显著顺式（cis）-eQTL信号的每个GWAS鉴定的候选因果基因位点上，一个变异是否同时负责GWAS信号和顺式eQTL信号。这些位点包括：aroDeH2 Zm00001d014734、hggt1、hppd1、por1、por2、samt1、vte1、vte4、vte5和vte7。在分析的23对基因-表型组合中，有18对存在1到6个SNPs，其共定位后验概率（CLPP）值——即一个变异同时是GWAS信号和eQTL信号的因果变异的概率——超过了严格的阈值0.01（补充图5和补充表16）。总体而言，这些分析选出了24个独特的SNPs，这些SNPs位于8个被研究基因的内部（por1和hppd1）或附近（aroDeH2 Zm00001d014734、hggt1、por2、samt1、vte4和vte5）。在候选因果基因内的3个SNPs中，只有一个被注释为错义变异（精氨酸->甘氨酸；GERP = 2），该变异是从por1中与总生育酚相关的SNPs中选出的（补充表14）。其他21个被选出的SNPs与它们各自的候选因果基因的中位距离为9.3kb。共定位的证据最为强烈的是，位于por2基因603bp处的一个SNP（5_25434949）具有最高的CLPP值（0.12-0.17），分别对应于γT、RT和RTT3。这与该SNP同时是por2的GWAS（混合线性模型：γT、RT和RTT3）和顺式eQTL信号的峰值标记是一致的（图4和补充表9、15）。我们还观察到，在hggt1的δT3中存在额外的统计学证据支持等位基因异质性（即具有CLPP值≥0.01的独立SNPs），但在vte4中没有观察到，因为对于αT、αT3、δT、γT和γT3，vte4的同一个SNP（5_205853870，距离vte4 25kb）被选出，其CLPP值为0.06-0.14。另外两个在GWAS中检测到且具有显著顺式eQTL信号的位点是vte1和vte7，但这两个位点的SNPs的CLPP值均<0.01。鉴于这一发现以及这两个位点在TWAS中未被检测到，vte1和vte7的因果变异可能与籽粒生育烷醇的遗传调控和基因表达变异不同。

讨论

通过结合跨越生物组织层次的互补数量遗传学方法，我们在包含近1500份自交系的玉米阿姆斯种群中进行了全面研究，以阐明籽粒中生育烷醇水平自然变异的遗传基础。对阿姆斯种群进行1200万个SNP标记的填补，生成22000个基因在籽粒发育的生物学信息阶段的转录本丰度，以及对3500多个成熟籽粒样本进行高效液相色谱（HPLC）分析，使得我们能够实施GWAS（全基因组关联分析）、TWAS（转录组全关联研究）、FCT（费舍尔联合检验）和eQTL（表达数量性状位点）分析方法。这些方法使我们能够鉴定出13个候选因果基因位点，这些基因编码的功能负责调节籽粒中生育烷醇的含量和组成。这13个基因包括一个S-腺苷甲硫氨酸转运蛋白（samt1）以及参与生育烷醇尾部合成（por1、por2、vte5、vte7、dxs1和dxs2）、芳香族头部基团合成（2个aroDeH2拷贝和hppd1）和核心生育烷醇途径（vte1、vte4和hggt1）的基因。在鉴定出的13个基因中，有8个被3种（GWAS、TWAS和FCT）遗传定位方法中的2种或更多方法检测到（表2）。这表明调控变异的重要性，因为13个基因中有7个通过TWAS被检测到，其中2个基因（dxs1和dxs2）仅通过TWAS被检测到。总体而言，整合的GWAS和TWAS方法增强了我们在阿姆斯种群中以基因水平解析生育烷醇籽粒表型的努力。

我们的多组学方法应用于具有高等位基因多样性和快速连锁不平衡（LD）衰减的阿姆斯种群（Romay等，2013），显示出高分辨率和统计功效。在鉴定出的13个基因中，有6个（samt1、aroDeH2 Zm00001d014737、dxs1、vte1、vte5和vte7）未在美国玉米NAM种群（Diepenbrock等，2017）中被鉴定到。值得注意的是，除了1个基因（vte1）外，其余5个基因之前均未与玉米籽粒中的生育烷醇相关联（Li等，2012；Lipka等，2013；Diepenbrock等，2017；Wang等，2018）。vte5的关联信号在阿姆斯种群中是新的，这可能是由于美国NAM种群仅从26个多样化的亲本中采样，限制了其等位基因多样性的代表性。

在美国NAM种群中，samt1、dxs1和vte7分别位于尚未解析的NAM联合连锁分析数量性状位点共同支持区间（JL-QTL CSI）内。然而，aroDeH2 Zm00001d014737（通过GWAS和FCT检测到）位于一个与aroDeH2 Zm00001d014734相关的NAM JL-QTL CSI内（Diepenbrock等，2017）。vte1位于一个包含hppd1的NAM JL-QTL CSI内，但由于其所在的着丝粒周围区域重组受到抑制，限制了定位精度，因此无法明确鉴定vte1位点（Diepenbrock等，2017）。然而，阿姆斯种群提供了更高的分辨率，能够在基因水平上解析hppd1和vte1。此外，美国NAM种群中鉴定出的其他6个较小效应的籽粒生育烷醇相关位点（sds、vte3、snare、ltp、phd和fbn）在阿姆斯种群中未被鉴定到，这表明尽管阿姆斯种群具有更高的物理分辨率，但在鉴定小效应QTL方面，其统计功效可能低于NAM设计。

本研究鉴定出首个与籽粒生育烷醇相关的代谢物转运蛋白Zm00001d017937，它编码玉米中拟南芥SAMT1的同源蛋白，负责将S-腺苷甲硫氨酸（SAM）运输到质体中（Bouvier等，2006；Palmieri等，2006）。玉米samt1与δT、δT3和αT3相关（表2），这与两种质体定位的生育酚甲基转移酶（VTE3和VTE4）依赖其从细胞质中导入它们的辅底物SAM是一致的（Bouvier等，2005；DellaPenna，2005）。在拟南芥中敲除SAMT1以及在N. benthamiana中沉默SAMT1，导致叶片中αT和总叶绿素（其合成也需要SAM）水平降低，而γT含量增加（Bouvier等，2006），这再次与SAM导入减少限制了VTE4将γT转化为αT（以及将δT转化为βT）是一致的（Bouvier等，2006）。鉴于TWAS揭示的samt1表达与δT浓度之间的强负相关关系（补充表17），我们推测samt1的弱表达通过降低VTE4和/或VTE3的SAM依赖性活性，导致δT的积累增加（图1）。

鉴定出的por1、por2、vte7和vte5位点主要影响生育酚，这与它们都参与生育酚合成的PDP（植物二烯磷酸）生成是一致的（Valentin等，2006；Diepenbrock等，2017；Albert等，2022）。这两个por位点是玉米NAM种群中籽粒总生育酚（RT）的两个最大效应的联合连锁分析数量性状位点（JL-QTL）的基础（Diepenbrock等，2017），它们编码叶绿素生物合成中一个高度调控活性的关键酶。Diepenbrock等（2017）发现por1和por2是RT的表达和效应QTL（ceeQTL）；即这些QTL对该性状的联合连锁等位基因效应估计值与相应基因在多个籽粒发育时间点的表达水平显著相关。与这些发现一致的是，por1和por2在本研究中也是与RT通过TWAS和FCT最强烈相关的两个基因。在阿姆斯种群中，通过GWAS（5% FDR）以及更弱地通过TWAS（最显著基因的0.6%）定位到δT的一个关联信号位于vte7。最后，vte5是生成PDP所需的两个激酶之一（Valentin等，2006；Vom Dorp等，2015），在GWAS中与玉米中的总生育烷醇（RTT3）相关。在拟南芥中，VTE5是与种子RT最强相关的基因之一（Albert等，2022）。本研究中por1、por2、vte7和vte5与玉米籽粒生育酚性状的关联，以及拟南芥中VTE7和VTE5与种子生育酚性状的关联，进一步加强了生育酚与叶绿素生物合成在种子中的联系，并表明其是单子叶植物和双子叶植物种子生育烷醇变异的主要调控点。

dxs1和dxs2基因仅在TWAS中与生育三烯酚性状相关，而与生育酚性状无关。这两个基因编码甲基赤藓醇-4-磷酸（MEP）途径的第一个且关键的步骤，表明它们是提供用于生育三烯酚合成的IPP（异戊二烯焦磷酸）的主要调控点。除了我们发现dxs2转录本丰度与生育三烯酚水平相关外，Diepenbrock等（2017）还发现dxs2是生育三烯酚性状（δT3、γT3和RT3）的ceeQTL（表达和效应QTL）。我们对dxs2表达变异的eQTL定位分析鉴定出一个位于第6染色体上的主要反式（trans）-eQTL（补充表14），该eQTL定位到八氢番茄红素合成酶1（psy1，Zm00001d036345）。psy1编码类胡萝卜素生物合成中的第一个且关键的步骤（Hirschberg，2001）。在拟南芥中，通过PSY活性对DXS蛋白水平的反馈调控来控制进入类胡萝卜素途径的流量（Rodríguez-Villalón等，2009）。因此，尽管我们没有发现psy1与生育三烯酚显著相关，但psy1可能间接影响玉米籽粒中的生育三烯酚水平。

与dxs1和dxs2类似，鉴定出的3个参与芳香族头部基团生物合成的基因也仅与生育三烯酚性状相关，包括2个aroDeH2拷贝（aroDeH2 Zm00001d014734和aroDeH2 Zm00001d014737）和hppd1。aroDeH2催化芳氨酸氧化脱羧生成酪氨酸，后者再转氨生成对羟基苯丙酮酸（HPP），即HPPD的底物。这两个aroDeH2基因相隔18kb，它们的GWAS关联峰值SNP相同（5_61159296），该SNP距离Zm00001d014737有40kb，距离Zm00001d014734有60kb。只有Zm00001d014737通过了FCT阈值（<0.3% vs >1.2%用于Zm00001d014734），主要是因为它更接近具有较低P值的GWAS SNPs。Zm00001d014734被鉴定出一个顺式（cis）-eQTL（补充表15），而Zm00001d014737缺乏显著的顺式或反式eQTL。目前来自美国NAM种群和阿姆斯种群的研究证据尚不能排除这两个aroDeH2基因都参与头部基团生物合成的可能性。

hppd1编码的酶能够产生HGA（对羟基桂皮酸），这是合成生育酚和生育三烯酚所用的芳香族头部基团。与Diepenbrock等（2017）的发现一致，我们发现hppd1与生育三烯酚的关联性很强，但在美国NAM种群中，hppd1与生育酚的微弱关联在阿姆斯种群中未能重现，这可能是因为阿姆斯种群中生育酚性状的等位基因效应大小未达到显著水平。与Diepenbrock等（2017）的发现一致，他们发现hppd1并非ceeQTL（表达和效应QTL），在阿姆斯种群中，hppd1也未通过TWAS在0.5%的阈值下被检测到。

鉴定出的3个基因（vte1、vte4和hggt1）编码生育烷醇途径中的核心活性，其对生育烷醇表型的影响与其酶活性一致（Shintani和DellaPenna，1998；Cahoon等，2003；Hunter和Cahoon，2007；DellaPenna和Mène-Saffrané，2011）。VTE4能够将γ-和δ-生育酚甲基化，生成α-和β-生育酚（Shintani和DellaPenna，1998；DellaPenna和Mène-Saffrané，2011）。vte4在阿姆斯种群中与γ和α同分异构体的关联性很强，这与之前的GWAS报告一致（Li等，2012；Lipka等，2013；Diepenbrock等，2017；Wang等，2018；Baseggio等，2019），并且与δT（GWAS和FCT）和RTT3（GWAS）的关联性较弱。同样，hggt1编码生育三烯酚生物合成中的关键步骤，并且与所有生育三烯酚性状显著相关，这再次确认了hggt1是控制生育三烯酚变异的关键基因（Lipka等，2013；Diepenbrock等，2017；Baseggio等，2019）。然而，之前在美国NAM种群中检测到的hggt1与生育酚的相对较弱的关联在阿姆斯种群中未能重新识别（Diepenbrock等，2017）。vte4和hggt1均被Diepenbrock等（2017）指定为ceeQTL，并且在阿姆斯种群的TWAS中，它们是检测到的与αT、αT3和RT3相关的最高排名基因，这突显了核心生育烷醇途径中调控变异的重要性。VTE1能够将MPBQ和MGGBQ分别转化为δT和δT3，我们鉴定出vte1与δT和δT3的关联在更高分辨率上（即更接近基因的显著变异）被识别，这比其他多样性种群的研究结果更精确（Lipka等，2013；Baseggio等，2019）（表2和补充表9）。

通过采用概率方法整合GWAS和eQTL结果的统计分析，我们鉴定出一组24个独特的SNPs，位于8个候选因果基因位点上，这些SNPs同时贡献了GWAS信号和顺式（cis）-eQTL信号，为通过eQTL介导的机制影响生育烷醇水平提供了支持。其中21个被选出的SNPs位于这8个基因的外部（最远达136.5kb），其中18个位于相应基因5’端的上游。尽管这些选出的SNPs本身可能并非因果变异，且可能与包括拷贝数变异在内的因果变异处于强连锁不平衡（LD）中，但这些结果暗示，共定位信号的潜在因果变异主要位于启动子和上游顺式调控元件附近。

尽管顺式（cis）-eQTL几乎总是被检测到的候选因果基因表达变异中效应最大的QTL，但我们的研究受到统计功效和多重检验负担的限制，因此无法更完整地解析反式（trans）-eQTL（Albert等，2018）。然而，顺式效应在不同环境下的稳定性通常高于反式效应，因此在基因组预测模型中以单倍型水平纳入时，顺式作用的因果变异有可能在不同种群之间更具可转移性（Giri等，2021）。