摘要：● 本研究基于绵羊9个重要组织的多维度（RNA-Seq、ATAC-Seq、CUT&Tag和Hi-C）数据，绘制出国际首张绵羊多组织表观遗传调控图谱，对绵羊基因组进行了系统注释；并结合全球不同尾型的绵羊品种和通过10年积累构建出的“两千级”规模的羊群列队，首次在

绵羊多组织表观基因组综合图谱：复杂性状、驯化和育种的资源

iMeta主页：http://www.imeta.science

研究论文

● 原文: iMeta (IF 23.8)

●

原文链接DOI: https://doi.org/10.1002/imt2.254

● 2024年12月15日，兰州大学王维民教授团队在iMeta在线发表了题为“

Comprehensive multi-tissue epigenome atlas in sheep: A resource for complex traits, domestication and breeding

● 本研究基于绵羊9个重要组织的多维度（RNA-Seq、ATAC-Seq、CUT&Tag和Hi-C）数据，绘制出国际首张绵羊多组织表观遗传调控图谱，对绵羊基因组进行了系统注释；并结合全球不同尾型的绵羊品种和通过10年积累构建出的“两千级”规模的羊群列队，首次在绵羊上采用多维整合组学策略鉴定出影响绵羊尾脂沉积的关键基因（BMP2）及因果变异（Chr13:51760995 A > C。填补了绵羊这一物种在多组织表观遗传调控图谱领域的空白，为绵羊育种提供了宝贵的数据资源和基因素材。

● 第一作者：张德印、程江博、李晓龙

● 通讯作者：王维民（wangweimin@lzu.edu.cn）

● 合作作者：黄凯，袁律峰，赵源，徐丹，张煜坤，赵利明，杨晓斌，马宗武，徐全忠，李冲，汪晓娟，郑琛，唐德富，年芳，乐祥鹏，李万宏，田慧彬，翁秀秀，胡鹏，冯园庆，Peter Kalds，姜志华，赵云霞，张小雪，李发弟

● 主要单位：兰州大学草地农业科技学院、甘肃农业大学动物科学技术学院、华中农业大学农业动物遗传育种与繁殖教育部重点实验室、华盛顿州立大学动物科学系和生殖生物学中心、崖州湾国家实验室、宾夕法尼亚大学遗传学系、上海海洋大学水产种质资源与利用教育部重点实验室、甘肃农业大学理学院、中国农业科学院兰州兽医研究所

亮 点

● 基于9个重要组织的92个转录组和表观基因组数据集，构建了首张绵羊多组织表观遗传调控图谱;

● 鉴定出753,723个非冗余功能元件，其中60%以上为首次发现，主要包括与感知能力、免疫反应和尾部脂肪沉积相关的组织特异性启动子和增强子；

● 利用多维组学数据集鉴定出一个新的因果变异Chr13:51760995A>C，位于激活增强子区域，影响绵羊尾部脂肪沉积。

摘 要

全面的功能基因组注释对于阐明畜禽农艺性状的分子机制至关重要，但目前仍缺乏对绵羊基因组功能元件的系统注释。在本研究中，生成了9个重要组织的92个转录组和表观基因组数据集，和来自全球29个品种的2,357个个体的全基因组数据，以及与尾脂重性状相关的4,006条表型数据。从功能元件、染色质状态及其功能等方面构建了首张多组织表观基因组图谱，并探讨了功能元件在解释绵羊驯化和改良过程中表型变异方面的作用。特别是，总共鉴定了753,723个非冗余功能元件，其中60%以上为首次发现。并发现与感知能力和免疫反应相关的组织特异性启动子和增强子在受驯化影响的基因组区域中高度富集，而背最长肌组织特异性活性激活增强子和尾脂组织特异性激活启动子在受选育影响的基因组区域中高度富集。值得注意的是，基于多组学数据集，一个位于激活增强子区域的变异Chr13:51760995A>C被鉴定为影响尾脂沉积的因果变异位点。总之，本研究为基于整合多组学数据集解析绵羊复杂性状提供了基础资源和成功案例。

视频解读

Bilibili：https://www.bilibili.com/video/BV1L3BKY5EwC/

Youtube：https://youtu.be/sTnX8C75Pkc

中文翻译、PPT、中/英文视频解读等扩展资料下载

请访问期刊官网：http://www.imeta.science/

全文解读

引 言

绵羊（Ovis aries）作为最重要的农业动物之一，为人类提供了肉、奶和羊毛等重要资源。尽管绵羊具有重要的经济价值，但对其适应性进化和重要农艺性状的分子遗传机制仍知之甚少。全基因组关联研究（GWAS）已被广泛应用于解析人类疾病和畜禽复杂性状的遗传结构。但大多数与复杂性状相关的候选位点位于基因组的非编码区域，且存高度连锁，这成为解释GWAS结果的主要挑战，严重阻碍了对绵羊表型变异、适应性进化和重要农艺性状分子机制的理解。

继DNA元件百科全书（ENCODE）和表观基因组路线图（Epigenome Roadmap）项目之后，农场动物基因型-组织表达（FarmGTEx）项目和动物基因组功能注释（FAANG）计划显著推动了家畜各组织功能元件的注释，为理解复杂农艺性状的遗传机制提供了丰富的资源。然而，目前对绵羊基因组中调控元件的注释仅限于少数组织，有关绵羊基因组功能元件的全面注释落后于猪、牛、鸡等畜禽和模式生物。这有限的注释阻碍了对绵羊复杂农艺性状的分子机制的理解。因此，迫切需要构建一个全面的绵羊基因组功能元件图谱，以识别与重要经济性状相关的候选因果变异。尾脂沉积是脂尾型绵羊品种的主要经济性状，但对其研究不足。减少尾脂沉积对提高养羊业的经济效益至关重要。然而，目前仅发现少数与尾脂沉积相关的关键基因及功能突变。目前的研究主要集中在易于测量的性状上，如尾部形状和尺寸，但忽略了关键指标——尾脂重。此外，基于全基因选择信号和GWAS识别的数量性状基因座区域（QTL）候选变异数量较多而难以解释，需采用多组学策略来精细确定影响畜禽复杂性状的因果变异，但这种策略尚未用于绵羊尾脂沉积的研究。

本研究生成了大规模绵羊群体列队的多组学数据集，包括9个重要组织的92个转录组和表观基因组（转座酶可及性染色质测序：ATAC-seq、染色质靶向切割并标签化测序：CUT&Tag和高通量染色体构象捕获：Hi-C）数据集、来自全球29个品种的2,357个个体的全基因组数据，以及与尾脂重性状相关的4,006条表型数据。研究目标是：（1）建立多组织表观基因组图谱；（2）探讨功能元件在解释绵羊适应性进化和复杂性状方面的作用；（3）基于整合多组学数据识别影响尾脂沉积的候选因果变异及其靶基因。总之，本研究丰富了绵羊基因组中功能元件的注释，确定了影响尾脂沉积的因果变异，为今后研究绵羊复杂经济性状的生物学机制提供了基础资源和成功案例。

结 果

绵羊多组织表观基因组图谱概述

为推进绵羊基因组的功能注释和探讨增强基因组注释对解释绵羊复杂性状的作用，研究小组对9个重要组织（软骨、盲肠、垂体、下丘脑、肝脏、背最长肌、瘤胃、脾脏和尾脂）进行了RNA-Seq、ATAC-seq、两种组蛋白修饰（H3K27ac和H3K4me3）的CUT&Tag和Hi-C测序，生成了92个表观基因组数据集（图1A），获得了149亿条clean reads，平均映射率为93.16%（表S1）。峰（peak）中reads百分比（FRiP）、相对链互相关系数（RSC）和归一化链互相关系数（NSC）均表明测序数据质量足以用于下游分析（表S2）。

为探讨增强基因组注释在解释绵羊复杂性状中的作用，整合了全球29个品种的364个全基因组测序（WGS）数据，用于群体结构和选择信号分析。同时收集了2,003只绵羊的胴体重和尾脂重表型数据以及WGS数据。经过变异检查和过滤后，共21,741,684个单核苷酸多态性（SNP）被用于后续的GWAS分析（表S3）。

在9个组织中，从ATAC-seq、H3K27ac和H3K4me3数据中分别鉴定出258,403、20,945和52,552个peaks；平均长度分别为692.92、2479.90和2029.79 bp（图S1A-D），这些peaks分别占整个绵羊基因组的3.38%、1.82%和3.72%（表S2和图S1E），ATAC peaks被注释到非编码区域，而H3K27ac和H3K4me3 peaks主要被注释到基因的启动子区域（图S1F）。此外，共鉴定出753,723个非冗余调控元件，包括40,606个潜在的启动子、158,677个潜在的增强子和846,058个开放染色质区域（图1B），其中超过84%和64%的增强子和启动子是本研究中新鉴定的（图1C），并对检测到的启动子和增强子与先前的数据比较，发现50.51%的启动子和31.35%的增强子与已发表肝脏数据中鉴定的重叠（图S1G）。

基于基因表达和表观遗传标记数据集评估了不同组织和检测方法之间的关系，层次聚类和主成分分析（PCA）结果清楚地反映了检测方法，其次是组织类型和生物学重复（图1D，S2），并且ATAC-seq、H3K4me3和H3K27ac数据集彼此之间表现出很强的正相关性（平均R = 0.68），三个活性调控标记peak分布在基因转录起始位点（TSS）周围的区域，而它们与RNA-seq数据集表现出弱正相关性（图1E和表S4）。此外，为绘制绵羊基因组的三维（3D）基因组结构，以两只绵羊的尾脂组织为代表性组织进行了原位Hi-C测序，并在25kb的分辨率下构建了全基因组互作图谱。结果发现，染色体之间的空间距离越近，相互作用程度越高，尤其是在同一个染色体中（图1F）。总之，基于上述数据集生成了高分辨率、多组织的表观基因组图谱（图1G），丰富了绵羊基因组中功能元件的注释。

图1. 绵羊多组织表观基因组图谱

（A）检测组织和数据集的示意图。左图为湖羊九个主要组织类型的示意图，中图为组织类型、试验群体和代表性品种，右图为本研究的表观基因组、表型和基因组样本数。（B）不同组织中启动子和增强子的雷达图。红线和蓝线分别代表启动子和增强子的数量。（C）新注释的启动子和增强子百分比（蓝色），与之前发表的绵羊肝组织数据重叠的百分比（橙色）。（D）不同组学、组织和两个生物学重复的Pearson相关性热图。（E）基因TSS周围平均表观遗传标记信号的分布。TSS，转录起始位点；TES，转录终止位点。（F）绵羊尾部脂肪组织的全基因组Hi-C互作热图。右上角的矩形表示11号染色体的放大图，青色点和蓝色矩形分别表示已识别的loop和TAD。（G）不同表观遗传标记在染色体上分布的Circos图。从最外层到最内层的轨道分别表示每条染色体的基因密度、SNP密度、RNA-seq、ATAC-seq和CUT&Tag（H3K4me3和H3K27ac）信号。Hi-C，高通量染色体构象捕获；TAD，拓扑关联域；SNP，单核苷酸多态性；RNA-seq，RNA测序；ATAC-seq，转座酶可及性染色质测序；CUT&Tag，染色质靶向切割并标签化测序。

跨多组织和物种染色质状态的预测和表征

利用ChromHMM软件的默认参数整合九种组织中的三个表观遗传标记数据集（ATAC-seq、H3K4me3和H3K27ac），共定义六种染色质状态，主要分为启动子（TssA和TssW）、增强子（EnhA和EnhAW）、ATAC岛和静止/抑制区域（图2A、B）。在九个组织中共鉴定1,936,312个非冗余染色质状态，平均长度介于906.88和25,332.1bp之间，分别覆盖整个基因组的5.9%和94.6%（图S3和表S5）。启动子状态（TssA和TssW）在TSS处表现出最高的富集度（图2C），并且与其他染色质状态相比，其基因表达平均水平最高（图S4）。此外，利用公开的DNA甲基化数据研究了组织中DNA甲基化和染色质状态之间的关系。结果发现激活启动子的甲基化水平最低（图2D），这证实了众所周知的基因表达和启动子甲基化水平之间的负相关性。同时，基于九种组织的六种非冗余染色质状态绘制了全基因组染色质状态图谱。正如预期的那样，抑制区域分布在整个基因组的大部分区域，除此之外其他染色质状态的分布相似（图2E）。还检查了不同组织中染色质状态的分布，结果表明，启动子状态（TssA和TssW）的活性变化较小，而增强子状态（EnhA和EnhAW）的活性在组织之间变化较大（图2F）。

为进一步探讨绵羊、牛和猪之间染色质状态的保守性，基于动物基因组功能注释（FAANG）和猪表观基因组项目的数据，使用相同表观遗传标记预测这些物种匹配组织类型的六种染色质状态，分析表明，这三个物种中相应的染色质状态中单个表观遗传标记表现出相似的发射概率（图S5）。此外，使用LiftOver软件（minMatch = 0.95）比较绵羊、牛和猪基因组之间激活启动子（TssA）和强激活增强子（EnhA）序列的保守性。结果显示，绵羊基因组中80.45%的TssA和81.99%的EnhA区域在牛基因组中保守，而猪基因组中只有25.87%的TssA和30.66%的EnhA保守（图2G）。同时，在绵羊、牛和猪的同源基因CD86附近均观察到保守和非保守的TssA和EnhA区域，其中H3K4me3和H3K27ac组蛋白修饰模式在这三个物种中是一致的（图2H）。这表明调控元件在进化关系更近的物种之间往往更保守，而高度保守的调控元件可能在不同物种之间发挥相似的功能作用。

图2. 九种组织的全基因组染色质状态景观

（A）染色质状态的名称、缩写和发射概率。（B）染色质状态的基因组覆盖率。（C）基因组注释中不同染色质状态的平均富集度。TSS/TES_1k，TSS和TES的上游和下游1kb区域。（D）染色质状态的平均甲基化水平。（E）不同染色质状态的全基因组图谱。TssA，激活启动子/转录本，TssW，弱启动子/转录本，EnhA，强激活增强子，EnhAW，弱激活增强子，ATAC_Is，ATAC岛，Quies，静止/抑制。（F）基于累积基因组覆盖率的染色质状态变异。虚线=0.75。（G）与绵羊相比，牛和猪基因组中TssA和EnhA的序列保守性。（H）绵羊、牛和猪CD86同源基因TssA和EnhA的保守性。H3K4me3和H3K27ac的CUT&Tag轨道中括号内的数字表示信号强度。

基因和染色质状态的组织特异性分析

调控元件是组织特异性基因表达的重要调控因子。为研究功能元件与组织特异性基因表达之间的关系，在九种组织中总共检测到8,334个组织特异性基因（图3A、B和表S6）。这些组织特异性基因的功能预测准确反映了各种组织的特定生物学功能（图3C和表S7）。基于实时荧光定量PCR（qRT-PCR）对组织特异性基因中代表基因的表达进行检测，验证了分析结果的准确性（图3D，S6）。此外与其他组织相比，组织特异性基因在相应组织中表现出更强的活跃染色质状态富集和较少的抑制区域（图3E、F）。同时，根据基因重叠EnhA的数量将基因分为三组。结果表明，与EnhA区域重叠较少的基因相比，EnhA区域重叠较多的基因表达水平较高，组织特异性较低（图S7A）。GO富集分析表明，只有一个EnhA重叠的基因与组织特异性功能相关，例如脾脏中的溶酶体活性和吞噬作用。相反，具有多个EnhA重叠区域的基因富集了更普遍的生物学功能（图S7B和表S8）。

为进一步探讨活性染色质状态的生物学功能，在九种组织中鉴定了184,019个组织特异性EnhA（图3G），发现组织特异性EnhA区域信号与组织特异性基因表达之间存在显著正相关（p

图3. 组织特异性基因和染色质状态的功能表征

（A-C）九种组织中组织特异性基因（Tau > 0.8）的数量、表达模式和基因本体论（GO）富集条目。TPM，每百万转录本。（D）组织特异性基因表达背景下的染色质状态图谱示例。显示了每个基因在九种组织中的RNA-seq信号，垂直刻度代表RNA-seq标准化信号范围从0到500。（E）9个组织的染色质状态在背最长肌组织特异性基因中的富集。红点代表背最长肌组织中的染色质状态。（F）九种组织中活性染色质状态（TssA和EnhA）在组织特异性基因中的富集。点代表来自对应组织的富集。蓝色虚线表示富集倍数=1。（G）组织特异性（TS）EnhA的数量，以及与组织特异性EnhA重叠的基因的GO富集。（H）与组织特异性EnhA重叠的基因表达模式热图。（I）不同组织的组织特异性EnhA中的motif富集及其序列的标识。（J）与不同组织特异性EnhA中富集motif相对应的转录因子表达模式热图。

染色质状态预测增强绵羊驯化和改良中表型变异的理解

为探讨功能元件是否在受驯化和改良影响的基因组区域中显著富集。首先对来自全球的364只绵羊进行基于最大似然估计的群体结构分析（图4A）。结果表明，所有个体被分为三个独立的簇：亚洲摩弗伦羊、中国地方绵羊品种和改良品种（图4B、4C、S10）。为剖析绵羊的驯化，计算了野羊和家养绵羊之间的固定指数（FST）、种群扩展单倍型纯合性（XP-EHH）和核苷酸多样性比率（π比率）。根据三种统计数据的前5‰窗口，总共确定了46个候选基因组区域（图4D和表S11,12），结果发现TssA、TssW和EnhA区域在受驯化影响的基因组区域中高度富集（图4E）。在组织特异性染色质状态富集分析中，发现与免疫功能相关的脾脏特异性TssA和EnhA以及与感知能力相关下丘脑特异性TssA的富集程度较高（图4F）。这一结果与绵羊在早期驯化过程中感知能力和生存能力被优先选择的事实相一致。

此外，对亚洲摩弗伦羊与中国地方绵羊品种（WN）、亚洲摩弗伦羊与改良品种（WI）和中国地方绵羊品种与改良品种（NI）进行FST、XP-EHH和π比率的比较分析，讨论组织特异性TssA和EnhA是否在受选择的基因组区域中富集。结果发现在WN、WI和NI比较中分别鉴别到49、53和13个候选基因组区域（图S11和表S13-18）。富集分析结果显示，在NI组中盲肠、尾脂特异性TssA和背最长肌特异性EnhA高度富集，而脾脏、瘤胃和尾脂特异性EnhA在WN和WI组中高度富集（图4G、H）。这与改良绵羊品种的育种目标相一致，改良绵羊品种的育种目标是更快的生长速度和减少尾脂沉积，而地方绵羊品种则是具有较强的适应性和抵抗力，总之，这些结果表明功能元件在解释绵羊的适应性进化和复杂农艺性状方面发挥重要作用。

图4. 绵羊驯化和改良过程中受选择区域的染色质状态富集

鉴定与绵羊尾部脂肪沉积相关的基因组区域

图5. 不同尾型绵羊的全基因组选择信号分析

利用多维数据集对与尾部脂肪沉积的候选因果变异进行优先排序

为确定与尾脂沉积相关的候选因果变异并探讨其对尾脂沉积的作用，本研究整合了不同尾型的选择信号分析、尾脂重性状的GWAS结果和绵羊尾部脂肪组织的表观基因组数据（即ATAC-seq、H3K27ac、H3K4me3和Hi-C数据）（图6C）。结果表明，71.3% GWAS的SNPs位于Chr13上的受选择区域，但仅3.59%和3.14%的SNPs位于至少一个个体的开放染色质区域（OCR）和H3K27ac峰中。值得注意的是，结果仅发现两个SNPs（Chr13：51760995A>C和Chr13：51825895G>A）位于OCR和H3K27ac峰。这两个SNPs也位于与骨形态发生蛋白2（BMP2）、LOC101117953和LOC101118027基因相同的拓扑关联域（TAD）中（图6C）。这结果表明这两个SNPs可能与绵羊尾脂沉积密切相关。

图6. 整合全基因组关联研究（GWAS）和表观基因组数据，定位与绵羊尾脂沉积显著相关的候选因果变异

为进一步探讨候选变异如何影响绵羊尾部脂肪沉积，本研究统计了Chr13:51760995A>C和Chr13:51825895G>A SNPs位点在338只不同尾型绵羊群体中的等位基因频率。与脂尾型绵羊品种相比，Chr13:51760995A和Chr13:51825895G的参考等位基因在瘦尾型绵羊品种中的频率更高（图7A、B）。携带参考等位基因（A/G）个体的尾脂重和尾脂相对重（尾脂重/胴体重）显著低于携带突变等位基因（C/A，pC位点，仅BMP2的表达在携带突变等位基因C的个体和携带参考等位基因A的个体之间存在显著差异（pA的突变不影响BMP2及其附近基因的表达（图7G，H）。这些结果表明Chr13:51760995A>C的SNP位点的改变影响BMP2基因的表达。

为验证BMP2基因对绵羊尾脂肪沉积的影响，从尾部脂肪组织中分离原代前脂肪细胞，诱导其分化为成熟的白色脂肪细胞，构建过表达载体pcDNA3.1-BMP2，检测过表达效率、增殖及分化潜能。结果表明，BMP2过表达并不影响细胞增殖，但增强其成脂能力，表现为脂滴聚集和甘油三酯含量增加，以及脂肪细胞标志基因表达增加（图7I-K，S13A-C）。相反，转染BMP2基因siRNA的细胞中脂滴聚集和甘油三酯含量与阴性对照组相比明显降低（pC是一个影响绵羊尾脂重性状的可靠候选因果变异。

图7. 候选功能位点鉴定及靶基因验证

（A，B）候选位点在不同尾型绵羊品种中的等位基因频率分布。蓝色表示参考等位基因类型，红色表示突变等位基因。AL，阿勒泰羊；Ht，和田羊；LH，大尾寒羊；LL，兰州大尾羊；Cele，策勒黑羊；STH，小尾寒羊；ST，短尾羊；TY，滩羊；UJ，乌珠穆沁羊；BAR，Barag羊；BAY，巴音布鲁克羊；MG，蒙古羊；OL，藏羊（欧拉型）；Va，藏羊（山谷型）；PT，藏羊（草原型）；AW，澳洲白羊；TE，特克塞尔羊；BL，边区莱斯特羊；DP，杜泊羊；EF，东弗里生羊；SA，南非肉用美利奴羊；ML，美利奴羊；PD，无角陶塞特羊；WS，白萨福克羊；Bs，黑萨福克羊。（C-F）雄性湖羊群体（n= 2,003）中Chr13:51760995A>C和Chr13:51825895G>A位点不同基因型的尾脂重和相对尾脂重。（G,H）Chr13:51760995A>C和Chr13:51825895G>A位点不同基因型中候选基因的表达水平。红色、蓝色和橙色点分别表示为纯合变异、杂合变异和纯合参考。（I）pcDNA3.1-BMP2质粒的过表达效率。（J）过表达载体转染的前脂肪细胞中的甘油三酯水平。（K）过表达载体转染的前脂肪细胞中的脂肪生成标志基因表达水平。（L）BMP2靶向siRNA转染尾部前脂肪细胞后的干扰效率。（M）BMP2靶向siRNA转染的前脂肪细胞中的甘油三酯水平。（N）BMP2靶向siRNA转染的前脂肪细胞中的脂肪生成标志基因表达水平。数据表示为平均值±平均值的标准误差，差异通过双尾t检验进行分析。*pp

讨 论

在本研究中，生成并描述了九种组织中功能元件的表观基因组图谱。这是迄今为止报道最全面的绵羊表观基因组参考数据集，包括753,723个非冗余调控元件和8,334个组织特异性基因。值得注意的是，与之前的数据相比，超过84%的增强子和64%的启动子是新发现的。同时，定义并表征了九种组织中的染色质状态，并对绵羊、牛和猪进行了染色质状态比较。此外，通过将功能元件与选择信号和大规模GWAS结果相结合，探讨了功能元件在解释绵羊适应性进化和尾脂沉积中的作用。该表观基因组图谱将为功能基因组注释、跨物种比较分析、GWAS结果的注释和验证、基因组选择和编辑提供重要资源。

阐明家畜适应性进化和复杂性状的遗传机制是一个活跃的研究领域。功能调控元件在适应性进化和复杂性状的遗传调控中至关重要。因此，首先探讨了在绵羊早期驯化过程，功能元件在受选择基因组区域中的潜在作用。结果表明，功能调控元件中在受驯化影响的基因组区域中显著富集，这与之前报道的结果一致。此外，在驯化过程中，脾脏特异的TssA和EnhA富集程度最高，其次是盲肠和下丘脑组织特异的TssA。该结果与先前的研究结果一致，即免疫力和感知能力可能是重要的生理表型，在早期驯化中被优先选择。此外，在改良和育种过程中，盲肠特异的TssA、瘤胃和背最长肌特异的EnhA在野羊或中国地方品种与改良品种的受选区域中富集程度最高。这可能是改良绵羊品种在生长速度和饲料效率性状方面高度选育的结果，而地方羊品种可能具有更强的适应性和抵抗力。瘤胃是反刍动物特有的器官，主要负责消化和营养吸收，因此与反刍动物的饲料效率密切相关。此外，在中国地方羊与改良品种的比较中，尾脂特异的TssA在候选区域中高度富集。这一结果可能反映了中国地方品种和改良品种之间存在明显的尾部形态差异。中国地方羊品种约60%的为脂尾型，而改良品种主要为瘦尾型。尾脂沉积是一个重要的经济性状，在改良和选育种过程中经历了高度的选择。过多的尾脂沉积过多不仅影响羊的交配和正常运动，还会影响羊的肉品质、饲料效率和经济效益。因此，为阐明尾脂沉积的分子机制，研究人员付出了巨大的努力来寻找与尾脂沉积相关的候选基因组区域或基因。然而，这些研究结果并不十分确定，可能是因为他们主要以尾长、尾型、尾椎数等间接指标作为目标性状，并采用RNA测序、全基因组选择信号扫描和GWAS方法进行研究。虽然也有少数研究将尾脂重作为目标性状，但研究样本量较小，且认为大多数候选变异的效应较小，变异之间的LD较高。

在本研究中，首先对全球168只脂尾型羊和170只瘦尾型羊进行了全基因组选择信号分析。然后，首次基于大规模试验群体（n= 2,003）对尾脂重性状进行GWAS分析，将尾脂沉积相关的基因组区域缩小至Chr13染色体上。尽管使用比较群体基因组学和GWAS缩小了基因组区域，但仍无法识别出致病变异及其靶基因。最后，通过对选择信号、GWAS结果、尾脂组织表观基因组数据和RNA-seq数据的综合分析，发现Chr13:51760995A>C的A等位基因是有利的SNP，因为它与尾脂重量减少和尾脂组织中BMP2表达降低有关。BMP2是一种转化生长因子β超家族成员，在脂肪形成分化中发挥重要作用。进一步的功能验证实验表明，BMP2过表达增加脂滴形成的潜力，但其抑制表达阻止脂滴形成，这一结果与先前的报道相一致。这些结果进一步表明，Chr13:51760995A>C位点是一个影响尾脂重性状的可靠候选致因变异。

总之，这一资源有助于确定绵羊农艺性状的候选因果变异，加速生物驱动的选择性育种以满足不断增长的全球人口的需求。值得注意的是，通过整合多组学资源确定的SNP（Chr13:51760995A>C）可能具有重要意义，特别是在脂尾羊品种中。这一结果提供了有价值的分子标记，可通过基因组选择或分子标记辅助选择选育出尾部脂肪沉积少的绵羊或瘦尾型绵羊，从而减少动物因去尾而遭受的痛苦，以提高健康和动物福利。尽管这项研究取得了重大成果，但仍存在一些局限性。首先，虽然从九种主要组织中生成了多组织表观基因组图谱，但它可能无法完全覆盖绵羊的所有生物学相关组织。其次，尽管通过综合多组学分析确定Chr13:51760995A>C位点是影响尾脂沉的候选因果变异，但需要进一步的功能验证（如基因敲除或基因抑制试验）来阐明Chr13:51760995A>C和BMP2在尾脂调控中的确切作用。还需要在其他绵羊品种或群体中进行额外验证，以评估这些发现的普遍性。

结 论

本研究建立了首张绵羊多组织表观基因组图谱，填补了现有绵羊图谱数据集的空白。这为解析绵羊的适应性进化和复杂经济性状提供了基础数据资源。同时，通过对多组学数据集的整合分析，鉴别了一个新的因果突变（Chr13:51760995A>C），该位点可作为绵羊育种计划中减少尾脂沉积的潜在遗传标记。总之，这些发现为绵羊育种提供了宝贵的资源和分子标记，有助于加速遗传改良。

方 法

样本收集和数据来源

为构建多组织表观基因组图谱，收集了两只6月龄半同胞湖羊公羔的18个组织样本（软骨、盲肠、垂体、下丘脑、肝脏、背最长肌、瘤胃、脾脏和尾脂）。

对于群体结构和选择信号分析，获得了364个基因组数据集，包括本研究中生成的5个品种50个个体的WGS数据（10个个体来自用于GWAS分析的群体），以及来自我们之前研究和公开数据库的314个个体（26个亚洲摩弗伦羊、13个中国地方绵羊品种的167个个体和10个改良品种的121个个体）的WGS数据（表S24）。

对于GWAS分析，本研究分7个批次（从2018年到2022年）构建了一个由2,003只湖羊公羔组成的综合队列，用于性能测定和样本采集。从国家绵羊、山羊核心育种场和规模较大湖羊场随机挑选出生日龄相近的羔羊，在断奶后转移至兰州大学民勤试验场，在相同的管理条件和饮食下饲养。直到6月龄，通过颈静脉穿刺采集全血样本于-20℃下储存直至基因组DNA的提取。然后，屠宰以测量其胴体重和尾脂重。本研究中的所有组织样本均用冰冷的磷酸盐缓冲盐水（PBS）清洗，并在屠宰后30分钟内收集。样品立即浸入液氮中并运送到实验室。最后，将样品放入-80℃超低温冰箱中保存直至进一步使用。

文库构建和测序

对于RNA-seq，按照制造商的说明，使用FastPure Cell/Tissue Total RNA Isolation Kit V2（RC112-00,Vazyme）从九个组织（每个组织两个重复）中分离总RNA。使用NanoDrop 2000分光光度计（Thermo Fisher Scientific, USA）评估RNA纯度和浓度，并且仅使用高质量RNA样本构建链特异性RNA-seq文库。使用Agilent 2100 Bioanalyzer（Agilent Technologies, Santa Clara, CA, USA）评估文库质量，合格的文库于武汉影子基因技术公司（中国武汉）的Illumina NovaSeq 6000平台上进行测序。

对于ATAC-seq，取约5mg冷冻组织样本在液氮中粉碎，置于1mL冰冷PBS中均质成细胞悬浮液。按照之前描述的方法并稍作修改，提取细胞核，使其与Tn5转座酶反应混合物在37℃下孵育1小时；采用DNA Purification and Concentration Kit（TD413; Genstone Biotech）进行纯化。接着使用NEBNext High-Fidelity 2X PCR Master Mix（M0541L; NEB）对转座后的DNA片段进行PCR扩增，再用Kapa Pure Beads（KS8002; Kapa Biosystems）进行纯化，最后于武汉影子基因技术公司（武汉，中国）的Illumina NovaSeq 6000平台上进行测序。

对于Hi-C测序，选择尾部脂肪组织样本，并根据先前描述的方法稍作修改后进行处理。简而言之，取约1g样品粉碎用1%甲醛（252549；Sigma-Aldrich）在室温下固定20分钟。再用0.125 M甘氨酸溶液（V900144；Sigma）室温下孵育5分钟以终止反应，交联的组织样品与Hi-C裂解缓冲液混合冰浴15分钟后，加入0.3%十二烷基硫酸钠62℃孵育10分钟后，加入Triton X-100（28314; Thermo Fisher Scientific）停止反应。再用AluI内切酶（R0137L;NEB）在37°C下消化染色质7小时后，加入dATP（N0440s;NEB）在37℃下孵育1小时，然后以含有生物素接头（M0202L;NEB）的T4 DNA连接酶16℃过夜连接DNA，最后，提取、纯化DNA并在Illumina NovaSeq 6000平台上测序。

对于WGS，按照制造商的说明，使用EasyPure Blood Genomic DNA Kit（TransGen Biotech，中国）从2,053个血液样本中基因组DNA。使用琼脂糖凝胶电泳和NanoDrop 2000分光光度计（Thermo Fisher Scientific, USA）评估DNA完整性和质量。为每个DNA样本构建双端测序文库，并在Illumina NovaSeq 6000平台（诺禾致源科技有限公司）上进行测序。

原始测序数据分析

总共生成并分析了2,135个新数据集，包括来自9个组织的RNA-seq（n=18）、ATAC-seq（n=18）、CUT&Tag（H3K27ac、H3K4me3和IgG，n=54）、尾部脂肪组织的Hi-C（n=2）和全基因组重测序（n=2,043）数据集。我们还统一分析了之前研究的314个全基因组重测序数据和三个肝脏组织的全基因组亚硫酸盐测序（WGBS）数据。详细分析方法如下：

对于RNA-seq数据，我们使用TopHat（版本2.1.1）软件的默认参数进行比对，使用Samtools（版本1.12）对BAM文件进行排序，利用featureCounts（版本2.0.3）量化比对到每个基因上的reads数。接下来，将原始count矩阵标准化为每百万转录本（TPM）以供后续分析。根据先前描述的tau指数方法鉴定其中的组织特异性表达基因；tau值 > 0.8的基因被定义为组织特异性基因。

对于ATAC-seq和CUT&Tag数据，使用Bowtie2（版本2.2.4）对测序数据进行比对；使用Samtools（版本1.12）和Picard（版本3.1.1）过滤低质量的reads、未比对的reads和PCR重复。利用bedtools（版本2.26.0）软件过滤线粒体reads。然后，使用MACS2（版本2.2.6）软件进行call peak，显著peak（p

对于Hi-C数据，使用HiC-Pro（版本3.1.0）流程对来自不同文库的序列进行比对和过滤。然后合并多个文库生成接触矩阵，并执行矩阵平衡迭代校正和特征向量分解（ICE）归一化。为研究全基因组染色质相互作用，在25kb分辨率下使用Juicebox（版本1.11.08）构建和可视化接触矩阵。使用HiCExplorer（版本3.7.2）检测显著的TAD，显著性阈值设置为0.05的错误发现率，并使用fanc（版本0.9.28）的fancplot函数可视化相应的交互热图和TAD边界。

对于WGBS数据，使用trim_galore（版本0.6.10）的默认参数过滤原始reads。过滤后的reads用Bismark（版本0.24.2）与参考基因组进行比对，并利用picard MarkDuplicates（版本2.27.5）删除PCR重复。使用bismark（版本0.24.2）的bismark_methylation_extractor函数提取CpG位点的甲基化水平。然后基于bedtools map函数（版本2.27.1）计算不同染色质状态的平均甲基化水平。

对于WGS数据，使用Burrows-Wheeler aligner（版本0.7.15）默认参数将过滤后的序列比对到绵羊基因组（Oar_rambouillet_v1.0）上。将获得的比对文件转换为BAM文件并使用Samtools（版本1.10）进行排序，使用Samtools和Picard删除PCR重复。使用Genome Analysis Toolkit（GATK）中的HaplotypeCaller函数进行SNP识别，得到每个个体的基因组变异识别格式（gVCF）文件。然后使用CombineGVCF将所有的gVCF文件合并，使用GATK中的GenotypeGVCF对SNP进行联合识别。最后使用VCFTools（版本0.1.16）对SNP数据集进行过滤，过滤参数为：-remove-indels，-minDP 5，-min-alleles 2，-max-alleles 2，-maf 0.05，-max-missing 0.8。所有SNPs过滤完成后，利用ANNOVAR基于绵羊参考基因组（Oar_rambouillet_v1.0）的基因注释进行注释，分为外显子（即同义或非同义SNP）、内含子、基因间、上游、下游区域以及剪接位点的变异。

染色质状态划分

使用ChromHMM（版本1.24）对9种组织的CUT&Tag（H3K27ac和H3K4me3）和ATAC-seq数据集进行整合分析；基于多元隐马尔可夫模型对绵羊基因组的染色质状态进行注释。由于六种状态可以代表三种表观遗传标记的最多组合，根据表观遗传修饰的富集预测了六种染色质状态模型：激活启动子（TssA）、弱启动子（TssW）、强激活增强子（EnhA）、弱激活增强子（EnhAW）、ATAC岛（ATAC_Is）和静止/抑制区（Quies）。使用OverlapEnrichment参数计算不同染色质状态和基因组区域的倍数富集程度。

跨物种染色质状态的保守性分析

基于动物基因组功能注释（FAANG）和已发表的文章中下载牛（肝脏、脾脏、脂肪、下丘脑）和猪（骨骼肌、肝脏、脾脏、脂肪）的H3K27ac和H3K4me3的ChIP-seq和ATAC-seq数据。牛和猪的数据分别比对到susScr11和bosTau9参考基因组。所有数据处理和染色质状态划分均按照上述绵羊的相同方法进行。为探索跨物种染色质状态的保守性，使用带有默认参数（minMatch=0.95）的LiftOver工具将绵羊基因组（Oar_rambouillet_v1.0）位置中的TssA和EnhA转换为牛/猪基因组（bosTau9/susScr11）。成功转换为绵羊基因组坐标的TssA和EnhA被视为序列保守。此外，如果它们对应的猪/牛同源序列中与猪/牛TssA和EnhA重叠，则它们被归类为使用保守。

组织特异性染色质状态的鉴定

为探索染色质状态的组织特异性，使用bedtools中的merge参数，在9种绵羊组织中生成非冗余染色质状态。然后，主要对TssA和EnhA进行特异性分析，将不同组织的染色质状态与非冗余染色质状态重叠。如果相交区域≥1，则该组织的非冗余染色质状态计算为1；否则，视为0。单个组织中的非冗余染色质状态1被认为是组织特异性染色质状态，所有组织中的非冗余染色质状态1被认为是常见染色质状态。对与组织特异性调控元件重叠的基因进行GO功能富集分析。使用HOMER（版本4.11）识别在组织特异性调控元件中显著富集的motif，并使用R包ComplexHeatmap展示p值最高的基序。使用RNA-seq获得的表达矩阵分析motif结合转录因子（TF）的表达。

组织特异性基因表达分析

为鉴定组织特异性表达的基因，根据先前文献中描述的方法，通过计算每个基因的Tau值来确定基因的组织特异性。Tau值的计算方法如下，其中x代表组织基因表达的TPM值，n代表组织数量。

使用R包ComplexHeatmap（版本2.14.0）对组织特异性基因簇进行可视化，在DAVID网站（https://david.ncifcrf.gov/）上对组织特异性基因进行基因本体（GO）富集分析。

群体结构分析

为探讨所有野羊和家养绵羊的群体结构，使用VCFTools（版本0.1.16）将VCF文件转换为PLINK文件，并使用PLINK（版本1.9）中的-indep-pairwise 50 5 0.2参数对常染色体的SNP进行过滤。过滤后的SNP用于系统发育树重建、群体遗传结构和PCA分析。使用PLINK计算成对遗传距离矩阵，并使用在线工具iTOL（https://itol.embl.de/）构建和可视化邻接树。使用默认设置的ADMIXTURE进行群体结构分析。假定祖先种群（K）的数量设置为2和3。此外，使用PLINK进行PCA分析，并使用R包ggplot创建第一和第二特征向量的图。

染色质状态的选择特征富集分析

GWAS

转录组测序和分析

基于Chr13:51760995A>C和Chr13:51825895G>A位点的基因型（Oar_rambouillet_v1.0），从2,003只雄性湖羔羊的综合队列中随机选取74个尾部脂肪组织样本。使用上述方法，从每个组织中提取总RNA并评估RNA质量。然后，由上海派森诺生物技术有限公司（中国上海）在Illumina NovaSeq 6000平台上制备RNA-seq文库并进行测序。使用上述过滤和分析方法获得每个个体的每百万转录本（TPM）值，并使用Prism（版本8.0.1；GraphPad）统计和可视化Chr13：51760995A>C和Chr13：51825895G>A位点不同基因型的候选基因表达。

BMP2 siRNA干扰和过表达

使用玛雅生物技术公司合成的两对小干扰RNA（siRNA）来敲低BMP2（表S25）。根据其干扰效率，选择si-BMP2-1序列进行后续分析。为构建BMP2过表达载体，根据克隆获得的BMP2编码区序列设计含有BamHI和EcoRI位点的正向和反向引物。将扩增的片段与线性pcDNA3.1连接。本研究中使用的siRNA和引物序列列于表S25中。

细胞分离与培养

CCK-8检测前脂肪细胞增殖

油红O染色

转染BMP2过表达载体或siRNA后，用PBS洗涤前脂肪细胞三次，并用10％甲醛固定60分钟。接下来根据制造商的说明，用油红O染色试剂盒对细胞进行染色，以观察脂滴聚集。在倒置荧光显微镜下观察染色细胞的形态变化并进行成像。然后，加入异丙醇以提取脂滴内容物；在510 nm处测量反映脂滴含量的显色光密度（OD）。

细胞甘油三酯含量分析

转染过表达载体或siRNA的前脂肪细胞，分化48小时后收集，用预冷的PBS洗涤三次。使用三酰甘油测定试剂盒（南京建成生物工程研究所，中国南京）评估细胞甘油三酯含量。在微孔板读数仪（Bio-Rad）上测量显色的OD；然后将值标准化为总蛋白质含量（以µg/mg为单位）。

RNA提取和定量实时逆转录PCR（qRT-PCR）

统计分析

使用R软件（版本4.2.1）和GraphPad Prism 8.0.1（GraphPad Software）进行统计分析和数据可视化。使用R中的cor函数计算组织、生物重复和测定之间的Pearson相关性。使用双尾t检验或单因素方差分析（ANOVA）评估组间差异。统计显著性用以下符号表示：*p≤0.05；**p≤0.01；***p≤0.001；****p≤0.0001。结果以平均值±SEM表示。

代码和数据可用性

引文格式：

Deyin Zhang, Jiangbo Cheng, Xiaolong Li, Kai Huang, Lvfeng Yuan, Yuan Zhao, Dan Xu, et al. 2024. “Comprehensive multi-tissue epigenome atlas in sheep: A resource for complex traits, domestication and breeding.” iMeta3. e254. https://doi.org/10.1002/imt2.254.

作者简介

张德印（第一作者）

● 兰州大学草地农业科技学院在读博士研究生。

● 研究方向为绵羊复杂性状的遗传调控机理解析，以第一作者在iMeta、Zoological research、Journal of Integrative Agriculture、Animal等期刊发表SCI论文10篇。

程江博（第一作者）

● 研究方向为绵羊经济性状表观遗传调控机制，以第一作者（含共同）在iMeta、Journal of Animal Science等期刊发表SCI论文7篇。

李晓龙（第一作者）

● 研究方向为绵羊功能基因组，以第一作者在iMeta，Frontiers in Veterinary Science，Frontiers in Genetics期刊发表SCI论文。

王维民（通讯作者）

● 兰州大学草地农业科技学院教授，博士生导师。国家高层次人才特殊支持计划青年拔尖人才、国家重点研发计划青年科学家项目主持人、“西部之光”人才培养计划-西部青年学者。

● 研究方向为绵羊功能基因组与生物育种。设计并牵头实施绵羊双万羊基因组计划，研发出“华羊芯”基因组育种芯片及配套遗传评估技术，选育出特色绵羊新品系4个和新品种群2个。兼任国家肉羊种业科技创新联盟秘书长、中国畜牧兽医学会畜禽遗传标记学分会理事，iMeta青年编委、《遗传》编委等。主持国家级、省部级项目18项，在Nature Genetics、iMeta、Microbiome等期刊发表论文60余篇，授权国家发明专利11件，获省级一等奖1项、二等奖2项、三等奖3项。

更多推荐

（▼ 点击跳转）

高引文章 ▸▸▸▸

iMeta | 引用15000+，海普洛斯陈实富发布新版fastp，更快更好地处理FASTQ数据

高引文章 ▸▸▸▸

iMeta | 兰大张东组：使用PhyloSuite进行分子系统发育及系统发育树的统计分析

高引文章▸▸▸▸

iMeta | 唐海宝/张兴坦-用于比较基因组学分析的多功能分析套件JCVI

1卷1期

1卷2期

1卷3期

1卷4期

2卷1期

2卷2期

2卷3期

2卷4期

3卷1期

3卷2期

3卷3期

3卷4期

3卷5期

1卷1期

1卷2期

期刊简介

“iMeta” 是由威立、宏科学和本领域数千名华人科学家合作出版的开放获取期刊，主编由中科院微生物所刘双江研究员和荷兰格罗宁根大学傅静远教授担任。目的是发表所有领域高影响力的研究、方法和综述，重点关注微生物组、生物信息、大数据和多组学等前沿交叉学科。目标是发表前10%(IF > 20)的高影响力论文。期刊特色包括中英双语图文、双语视频、可重复分析、图片打磨、60万用户的社交媒体宣传等。2022年2月正式创刊！相继被Google Scholar、PubMed、SCIE、ESI、DOAJ、Scopus等数据库收录！2024年6月获得首个影响因子23.8，位列全球SCI期刊前千分之五(107/21848)，微生物学科2/161，仅低于Nature Reviews，学科研究类期刊全球第一，中国大陆11/514！

“iMetaOmics” 是“iMeta” 子刊，主编由中国科学院北京生命科学研究院赵方庆研究员和香港中文大学于君教授担任，是定位IF>10的高水平综合期刊，欢迎投稿！

来源：微生物组