深度解读：发表Nat Commun后，如何再进一步？

绵羊从新月沃地驯化后，随着人类的迁徙扩散到全世界，经过对当地多样环境的适应，形成了大约1400个表型丰富的当地和培育品种，已经被培育成为一种用于产肉、产奶和产毛等用途的重要经济物种。利用高深度全基因组重测序数据，鉴定表型性状相关功能变异对指导基因组辅助育种有着重要意义。近期已有研究从基因组层面分析了绵羊的驯化和育种，然而依然存在样本量较少、测序深度低和表型数据不够充足等局限。

2020年6月，中国农大李孟华教授与新疆农垦科学院王新华研究员在Nature Communications在线发表题为“Whole-genome resequencing of wild and domestic sheep identifies genes associated with morphological and agronomic traits”的研究论文。该研究通过对16个野羊（亚洲摩弗伦）和来自世界各地的232个家羊个体进行高深度全基因组测序，全面深入地揭示了绵羊在驯化、培育改良以及表型性状选择上的遗传机制。

通过GWAS及转录组分析发现PDGFD基因是影响绵羊尾部脂肪沉积的重要因果基因。通过全基因组选择分析和基因组关联分析（GWAS），利用单核苷酸多态性（SNP）和拷贝数变异（CNV）分子标记，挖掘了一系列与驯化、重要表型形状和农业性状相关的候选基因，并检测到一些已经发生了明显的等位基因频率差异分化的重要候选基因区域；另外，通过转录组数据分析，反转录PCR（RT-PCR），实时荧光定量PCR（qPCR）和免疫印迹（Western blot）等实验验证，发现PDGFD基因是影响绵羊尾部脂肪沉积的重要因果基因。

该研究收集了覆盖亚洲、欧洲、非洲和中东的16个野羊亚洲摩弗伦、172个当地品种个体和60个培育品种个体（图1，收集样本地域覆盖情况），进行平均深度为25.7×的高深度重测序，基因组平均覆盖率达到98.27%，共计获得20.55Tb测序数据。获得单核苷酸多态性（SNP）和拷贝数变异（CNV）分子标记等变异数据。下表1为本研究收集到的样本信息。

图1. 覆盖亚洲、欧洲、非洲和中东的16个野羊亚洲摩弗伦、172个当地品种个体和60个培育品种个体分布图

表1.本研究样本信息

通过过滤低质量SNP后，保留了28.36 million个SNP，同时鉴定到13551个CNVs和28973个SVs。进一步与公共数据库中的SNP数据比较，发现增加了将近230000个新的SNP。如此高质量的基因组变异数据，为后续挖掘绵羊遗传结构和驯化研究提供了最重要的遗传数据基础。

基于28.38 million SNPs，分析三个群体中SNPs的数据情况。结果显示，在16个野羊亚洲摩弗伦（Asiatic mouflon）中发现占有23.27million SNPs，unique SNPs数量有12.06 million个；172个当地品种个体（landraces）SNP数量达到14.38 million个，unique SNP数量有1.06 million个；60个培育品种个体（improved breeds）SNP数量有14.01 million个，unique SNP有1.08 million个；分别将以上SNP数据在不同群体中开展比较分析：当地品种与培育品种共有（shared）SNP数量12.09 million，野羊与当地品种共有（shared）10.38 million个SNP，野羊与培育品种共有（shared）9.98 million个SNP，这与Fst和π的计算结果一致，证实了当地品种和培育品种间的基因组差异性均小于野羊与当地品种、野羊与培育品种的基因组差异，反应了绵羊驯化改良进程中基因组的变异规律。除了分析SNP，本研究还分析了Indel/CNVs/SVs等基因型数据，结果与SNPs数据基本一致。

图2. 不同类型群体特有和共有基因型数据韦恩图

为了研究绵羊古老的群体结构和统计学历史，研究者基于SNPs数据开展了以下分析。通过进化树分析、PCA分析和邻进树（neighbor-joining tree）分析发现，野羊与其他两类羊明显可以区分开，而且当地品种和培育品种进一步可以分为4个亚群，分别为欧洲地区、东亚地区、中东地区和非洲地区群体。其中，东亚地区和中东地区群体具有更近的遗传亲和关系。野羊LD衰减距离为2.8kb，地方品种和培育品种的的LD分别为12.1kb和17.1kb,这些结果表明人工选择和遗传分离导致了地方品种和培育品种的出现。见图3.

图3. 野羊和驯化改良品种群体结构分析

为了鉴定绵羊进化过程中，受到驯化选择的基因组区域，研究者选取了16个野羊个体和代表不同区域和遗传起源的5个地方小群体，通过XP-CLR分析，共计鉴定到302个可能存在的选择性清除位点，结合核苷酸多样性π值计算结果，筛选出144个可能存在的基因组选择区域，共覆盖261个候选基因。更进一步，整合iHS和HKA分析方法，分别鉴定到899和1503个选择清除位点。综上计算方法结果，依次将iHS和HKA与261个候选基因共分析后，分别有76和71个基因受到选择。通过比较这些驯化相关的选择位点与已知QTL结果发现，绵羊驯化进程中这些受到选择的基因组区域更多与产奶和产肉相关QTL重叠，这也反应了绵羊驯化过程中人们对羊奶和羊肉的需求。

接下来分析发现，驯化进程中受到选择的261个基因中，14个基因的功能在前期的研究中均有报道，包括SLC11A1, HOXA11, CAMK4, LEF1, TET2, KDR, CTBP1, GAK, CPLX1, PCGF3, FLT1, BCO2, CHGA, and HTRA1等基因，这些基因与繁殖育性、抗感染、骨形成、脂肪沉积等生产性状相关；22个基因报道在其他物种驯化进程中受到选择，主要包括牛、山羊、马、猪、狗、鸡、兔和小鼠等这些物种，这些基因的功能主要与动物的免疫力、色素沉着、毛色、感光体、行为/生长和繁殖特性等相关。除此之外，研究者发现5个候选基因内的非同义突变SNP等位基因频率在野羊和5个地方品种中存在显著差异，主要包PDE6B,BCO2, ADAMTSL3, NKX2-1, and an olfactory receptor 51A4-like gene LOC101108252，见图4。这些结果证实了绵羊在驯化进程中许多与农业生产性状的基因受到选择。

图4. 绵羊驯化改良进程选择位点分析

野羊驯化之后，形成各种地方品种，为了使绵羊适应不同的生长环境和特殊的生产系统，需要开展绵羊的育种和改良工作。本研究选取了36个地方品种和6个改良品种开展基因组比较分析，检测育种进程中受到正向选择的基因组区域和位点。通过Fst计算，鉴定到205个可能的受到选择基因组区域，累计23.80Mb，包含391个基因。基因功能注释证实，这些基因主要与重要表型和农业生产性状相关，包括是否存在羊角、色素沉着、繁殖和体型等。除以上基因外，也发现存在一些环境适应相关基因，如能量代谢、免疫响应等，这些基因可能是绵羊长期适应环境受到自然选择进化的结果。GO富集分析也进一步证实了以上环境适应相关的基因。其中，最高Fst值的基因组选择区域内，PAPPA2位于该选择区域，且该基因前期研究证实与人体的脂肪沉积和牛的产奶、繁殖和体型大小相关，见图5。

图5.育种和改良进程中受到选择的基因组区域和候选基因

绵羊的尾型是一种重要的表型性状，尾脂羊为满足人类早期的脂类需求提供了大量脂肪，而且也为绵羊抵御极端寒冷的环境储存了能量，但是随着人们生活水平不断提高以及对高脂类绵羊肉品需求的减少，过度尾脂沉积反而会增加饲养成本，因此，研究绵羊尾部脂肪沉积的分子调控机制对培育瘦尾以及低脂肉用绵羊具有重要意义。研究人员针对四种不同尾型（长脂尾、肥臀、短脂尾和短瘦尾）性状的家羊进行两两比较的全基因组选择分析，挖掘出PDGFD这个被共同选择的基因，进一步解析基因结构，发现该基因启动子区域的基因型分布，以及位于该基因的一个非同义突变位点（半胱氨酸/酪氨酸）的等位基因频率在肥尾/肥臀羊和瘦尾羊中差异显著。如图6。

通过后续的转录组分析以及RT-PCR、qPCR和western blot实验验证发现该基因在瘦尾羊中的表达显著高于肥尾羊。通过探索脂肪生成通路，发现PDGFD基因表达的蛋白通过形成二聚体来激活PDGFRβ通路，而该通路在脂肪扩增的过程中起着抑制作用，也进一步证实了PDGFD基因在尾脂形成中起着负调控作用。

图6. 尾型控制基因PDGFD的基因定位和遗传机制

除了尾型，绵羊还有其他重要表型（尾脂、角型、耳朵大小、羊角数量）以及农业性状（繁殖、产毛、产奶、产肉、羊毛细度、色素沉着、生长速率等）。为了定位到驯化改良进程中受到选择的基因组区域内与表型相关的基因，研究者对三个性状进行了GWAS（垫料大小、羊角和乳头的数量）分析，结果显示垫料大小、羊角数量和羊乳头数量分别定位到600、989和1969个显著信号。其中，进一步筛选到与垫料大小相关的候选基因 12个，包括NOX4, IRF2, PDE11A, ZFAT, ZFP91, TENM1, BICC1, LRRTM3, CTNN3, SMYD3, KCNN3和CD96；羊角数量相关候选基因为HOXD1, HOXD3和 RXFP2；羊乳头数量定位候选基因除了包括乳腺癌相关基因LRP1B,GRM3, MACROD2, SETBP1和GPC3外，还包括其他基因，如PHGDH, KDM3A, GLIS3, FSHR, CSN2, CSN1S1和ROBO2等，这些基因在小鼠研究中报道与乳头发育相关。

图7. 主要3个性状统计分析

该研究揭示了绵羊的驯化、培育、各种重要表型（尾脂、角型、耳朵大小）以及农业性状（繁殖、产毛、产奶、产肉等）的遗传机制，为绵羊遗传学研究提供了宝贵的基因组资源，对未来分子辅助育种和家养绵羊的遗传改良也具有重要指导意义。

粗读完李老师和王老师这篇NC大作后，想起了2018年黄三文老师和迈维代谢首席科学家罗杰老师发表在Cell的文章，两篇文章主要有两个区别：一是研究对象不一样，NC文章研究的是绵羊，Cell研究的是番茄；二是基因组学关联载体不一样，绵羊NC文章是以基因组学+农业生产性状为主线，番茄Cell则是以基因组学+代谢组学（代谢小分子物质）为主线；细细品来，两篇文章其实思路又如出一辙，均是以基因组学数据为基础，深入研究物种驯化改良进程中受到选择的某些重要的功能基因。

接下来简单比较一下两篇文章的具体内容。绵羊NC文章是一篇基因组测序文章的经典案例，文章思路缜密，数据庞大，是一篇不可多得的组学经典文章。文章与2018年发表在Cell杂志Rewiring of the fruit Metabolome in Tomato Breeding文章的群体构成和文章分析思路较为一致：首先是群体材料较为一致，NC文章群体主要包括16个野羊亚洲摩弗伦、172个当地品种个体和60个培育品种个体，共计248个个体，Cell文章群体材料共计442个个体，包括野生番茄31份、驯化后的小果番茄124粉，以及改良后的大果番茄287份。其次是分析思路上，NC文章围绕野羊由驯化成地方品种、地方品种选育改良为培育品种进程中，分析了绵羊的遗传结构及其驯化遗传机制，还重点解析了其丰富表型的遗传机制，包括定位到了绵羊的尾型控制基因PDGFD。Cell文章也是按照野生番茄驯化为小果番茄、再改良到大果番茄，揭示了番茄驯化改良进程中遗传结构和驯化机制，定位到了一批调控番茄风味物质相关的候选基因。

图8. Cell文章主要研究思路

但是，前期科研工作者在番茄驯化改良的研究中，研究者已经把最显著的农艺性状果重和果实大小基因成功定位（fw2.2, fw3.2和fw11.3等），且番茄果实的风味是重要的品质性状，风味又是由数以万计的小分子化合物决定的，因此，研究者转换了研究对象，由原来传统农艺性状转换到看不见、摸不着的品质和代谢性状（不同含量的代谢小分子物质）。于是提出了以下科学问题：番茄在驯化和改良过程中，营养价值和口感是否发生了变化？这种变化的营养和风味口感是由哪些小分子物质决定的？这些小分子物质的遗传调控机理是什么？这些遗传位点是否受到物种驯化改良的选择？在本研究中，重点分析了野生番茄为什么有苦涩的味道，这种涩味在驯化改良后的品种中为什么消失了？

最后，在Cell文章中，研究者发现番茄中造成这种涩味是因为一类物质——生物碱在果实中积累造成的，在野生材料中，生物碱这类物质含量比较高，驯化改良后这类物质含量显著降低，最后，发现这类物质定位的位点与受到驯化选择的基因组位点较为重叠，证实了生物碱在番茄驯化改良进程中受到自然和人工的选择。最后，研究者成功定位到控制这类物质合成调控的1个基因簇，包括10个候选基因。

同样，这些问题在绵羊等其他物种中均存在，值得大家去研究，只是前期没有关注到代谢组学上来。

图9. 生物碱类物质在3种群体中含量及GWAS定位位点

那么，在绵羊驯化改良进程中，是否有类似现象存在呢？决定绵羊羊毛颜色的物质是哪些？是否在驯化改良进程中受到选择？哪些基因控制这些物质？决定绵羊羊奶品质和口味的物质是哪些？这些物质是否在驯化改良中受到选择，受到哪些基因的调控合成？决定绵羊肉质质量和风味的物质是哪些？这些物质是否在绵羊驯化改良中受到选择，哪些基因决定这些物质的积累？绵羊驯化改良进程中能够适应环境的关键物质是哪些？这些物质是否在绵羊驯化改良中受到选择，哪些基因决定这些物质的积累？因此，对比Cell文章和以上有待进一步解决的科学问题，基于NC的结论和完整样本材料、测序数据，还有许多工作可以开展。那么，怎么找到以上那些问题的答案呢？解决问题的思路和技术就是结合当前热门的多组学技术-mGWAS。接下来，一起简单了解一下迈维代谢简单整理的mGWAS技术资料。

mGWAS（metabolome Genome-Wide Association Study）即为代谢物的全基因组关联分析，这里是把小分子物质作为代谢性状或者某些农业生产性状的代表表型。因为代谢物的种类和数量在不同品种与组织中存在很大差异。因此，可以利用二代测序技术获得群体材料的基因型数据（Genotype），结合代谢组数据（Metabolome），开展代谢物的全基因组关联分析，批量精确定位控制这些代谢表型和决定农业生产性状的小分子代谢物质（品质、营养、环境适应、功能营养分子等）相关候选基因，挖掘调控营养和品质等相关代谢通路机理，将有助于我们深入了解农业生产性状的代谢遗传机制，加快培育优质新品种。

什么是代谢组？代谢组是联系基因组和表型的桥梁。如下图10，中心法则揭示了基因组到蛋白质组的作用关系，而催化生物体内发生代谢反应的媒介大多数是酶类（多数是蛋白质），进而形成代谢小分子物质，进一步，小分子物质的积累又决定了个体性状的直接体现。如下表2，所以，基于基因组+代谢组的mGWAS研究思路是研究农业生产性状的可行方案。

图10. 代谢组是基因组与表型的桥梁

生物体的表型性状是由代谢物的积累直接体现的，当花青素的含量在甘薯块茎中过量积累时，块茎颜色将由红色变成紫色（颜色是园艺植物常见性状）；长春花长期以来被视为最有效的抗癌植物材料之一，而有效的抗癌活性物质（抗癌是一种生理活性性状）是长春花碱vinblastine（生物碱类物质）。因此，植物功能基因组的研究最终都会回落到代谢组与表型组研究中来。同样，动物研究中代谢小分子物质是动物多种性状的直接体现或间接作用，如下表2，迈维代谢小编汇总了动物研究中部分代谢物在营养品质、疾病和环境适应方向相关的作用：如尿素氮、甘氨酸、丝氨酸、胆碱和琥珀酸等物质决定了牛产奶的性状；三甲胺与鱼腥味、奶腥味相关；磷脂浓度、长链鞘磷脂含量与牛胰岛素功能紊乱疾病相关；甘油酯、甘油磷脂、鞘脂与小鼠温度适应相关。由此可看出，代谢小分子物质是农业生产性状的直接体现，这也是“基因组-代谢组-生产农业性状”多组学研究的理论基础和优势。因此，代谢物与基因组的关联分析（mGWAS）是研究农业生产性状的有效方法，越来越受到科研工作者的青睐。那如何开展mGWAS研究呢？

表2. 动物研究中部分小分子物质与农业生产性状的对应关系

简单来讲，mGWAS研究与普通GWAS的流程基本一致，主要区别在于：普通GWAS需要统计好各种农业生产性状表型数据，而mGWAS重点在于准确定量成千上万个代谢小分子物质在不同样本中的含量。当然，高水品的mGWAS文章不仅仅分析代谢物的基因定位结果，同时也会定位普通农艺性状的基因，最后把两者定位结果做整合分析，挖掘“共定位”的位点和基因。mGWAS主要思路见下图11，主要包括以下步骤：

图11. mGWAS研究技术路线

1. 选取群体材料：

根据研究目的和方向，选择一定的个体组成合理的群体结构，同时记录好农业生产性状（有助于后续mGWAS与农业性状GWAS结果共定位分析）数据；一般来讲，建议群体材料选取200个个体以上的群体，且群体包括各种代表性品种，例如野生型、地方品种、栽培品种、育成品种等，这些材料的组成决定了后续文章分析的思路。本篇绵羊NC文章中群体组成合理，包括了野羊、地方品种和培育品种，数量达到了248个个体，从群体结构来看，该群体材料非常完整，这些为后续提炼文章的“代谢生物学故事——Domestication和improved”提供了基础。

2. 构建绵羊专属代谢物质库：

mGWAS是把代谢物作为代谢表型，开展代谢物+基因组的关联分析，因此，代谢物的重要性和准确定量决定了后续数据分析的质量“好与坏”。这里的“好”，是指既要鉴定注释出到绵羊中重要的、具有代表性的、能侧面反馈农业生产性状的代谢物，也要开发出高通量、准确定量检测代谢物含量的检测技术。为了建立物质库立和准确定量代谢物含量，迈维代谢整合了市场上的非靶向代谢组学技术（untargeted metabolomics）和靶向代谢组学技术（targeted metabolomics）优缺点，于2013年由迈维代谢科学家开发了一种创新的代谢组检测技术——广泛靶向代谢组（widely -targered metabolomics）（chen.,Molecular Plant,2013）。截止目前，基于广泛靶向代谢组技术，迈维代谢建立了不同物种的代谢物质数据库（MWDB物质数据库），这为后续的功能基因组学研究提供坚实的“物质基础”。下图展示了迈维代谢动物方向部分物质库：动物广泛靶向数据库、动物广靶脂质数据库和动物氧化脂质数据库。

图12. 迈维代谢动物相关代谢物质数据库

3. 采集群体基因型数据和代谢组数据：

基于二代测序获得基因型数据，绵羊NC文章中已经产出20.55Tb高质量基因组测序数据，平均每个样本测序深度达到25.7X，这些高质量的基因型数据（包括SNP/CNV/SV等）已满足mGWAS关联分析所需的基因型数据。基于上述建立的绵羊专属代谢物质库，利用QQQ仪器定量金标准方法——MRM检测模式检测数据库所有物质在248个个体中每个代谢物的含量，并通过不同的质控（TIC重叠图/QC mix计算/标准液批次矫正）方法确证代谢组数据采集的准确性，为后续的代谢物GWAS计算提供准确的代谢表型数据，确保准确定位目标基因，降低假阳性率。

4. 功能基因定位：

基于重测序数据和代谢物含量开展GWAS分析，筛选每个物质的显著性SNP和候选基因，并结合其他数据（转录组数据/表型数据）获得候选基因，最后验证基因功能准确性。这块的重点和难点在于筛选候选基因，GWAS关联方法可以用普通的MLM模型，因为大量的代谢物检测并关联后，获得到大量的显著性SNP，每个显著性SNP在一定基因组区域内（如上下游500kb）会有一些基因，因此，这里会得到数量较大的候选基因。迈维代谢基于丰富的项目经验总结了一套挑选候选基因的方法，最终一般来讲，mGWAS文章中会列出20-40个定位到的基因，验证3-5个候选基因，可以做转基因验证和体外蛋白表达验证。

5. 代谢生物学与科学问题：

得到靶基因和相对应的代谢物质后，需要对这些基因和物质开展功能研究和数据挖掘。这里简单罗列几点：

• 筛选特定基因决定物质在不同品种含量差异积累的功能性SNP；

• 绵羊在驯化改良进程中显著改变的代谢物和对应的基因，这些基因与自然选择的关系；

• 不同物质具有不同类型的生物学功能，因此，除了研究代谢物相关基因，还可以结合绵羊各种重要表型（尾脂、角型、耳朵大小）以及农业性状（繁殖、产毛、产奶、产肉等）定位的位点，从“基因组-代谢物-农业性状”解析重要性状的遗传机制。

除此外，还可以从其他方向开展代谢组数据和故事的提炼。

1.A novel integrated method for large-scale detection, identification, and quantification of widely targeted metabolites: application in the study of rice metabolomics. Molecular plant. 2013

2.Rewiring of the Fruit Metabolome in Tomato Breeding. Cell. 2018

3.Whole-genome resequencing of wild and domestic sheep identifies genes associated with morphological and agronomic traits. Nature Communications. 2020

4.A cytochrome P450 regulates a domestication trait in cultivated tomato. PNAS. 2013

5.Tomato fruit weight 11.3 maps close to fasciated on the bottom of chromosome 11. Theor Appl Genet. 2011