iMeta | 中国农科院唐中林组-解析猪表型分化和复杂性状形成的遗传机制

多组学分析揭示与猪复杂性状相关的选择信号和基因座特征

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

● 结合9个高质量基因组和1081个全基因组重测序数据，构建了较为完整的猪泛基因组及遗传变异数据集，为解析猪的复杂性状奠定了基础;

●  发现了丰富的与猪脂肪沉积、肌肉生长、免疫性状相关的选择性扫描信号和基因，为未来的育种计划提供了关键的遗传标记；

●  整合27个骨骼肌发育阶段的转录组和甲基化组数据，揭示了东西方猪种骨骼肌生长差异的表观遗传和异时基因调控机制；

●  人工选择对DNA甲基化和异时基因调控有显著影响，与东西方猪种骨骼肌生长和肉质的差异密切相关。

引  言

猪（Sus scrofa）约1万年前在安纳托利亚和中国南方分别被独立驯化为重要的全球农场动物。自然选择和人工选择的综合作用促成了约600个不同品种的形成，这些品种在表型多样性和环境适应性上表现出显著差异。尤其是在过去100年中，西方猪种因体重、生长速度和肌内脂肪含量等经济性状被强烈选择，与东方猪种相比呈现显著差异。除提供丰富的农产品外，猪还作为人类疾病的生物医学模型，并在异种移植中展现出潜力，例如作为器官供体。因此，研究猪驯化和育种背后的遗传变化和选择机制具有重要的生物学、医学和经济价值。

在过去十年中，多个研究揭示了猪驯化的关键遗传结果，并鉴定出多种与表型多样化相关的数量性状位点（QTL）和基因，如EPAS1、ESR1、KIT、LCORL、MC1R、MUC13、NR6A1、PLAG1、RYR1和VRTN。这些关键基因的发现加速了猪育种进程，但仍有许多遗传决定因素未被揭示。为了填补这一空白，需要更大规模的样本和涵盖更多猪种的全面遗传调查。此外，越来越多的证据表明，由于调控元件的突变缺乏有害的多效性效应，其在驯化过程中的作用尤为重要。这些突变会影响关键发育基因的时空表达模式，即使是细微的转录变化也可能影响特定器官的发育时间，长期来看可能导致表型差异。

为了深入理解支持经济性状的（表观）遗传基础及其对猪形态发生的贡献，我们新生成了七个高质量的从头组装的基因组，并联合分析了来自1,081个个体的全基因组重测序数据以及覆盖东西方猪种的27个发育阶段的骨骼肌转录组和甲基化组数据。在这项工作中，我们构建了全面的泛基因组和遗传变异库，揭示了许多与经济重要性状相关的选择性扫描信号，强调了与肌生成和肉类生产性状相关的异时性基因及人工选择下的DNA甲基化变化。这些研究结果为东西方猪的（表观）遗传及表型分化提供了新见解，并提出了一些未来猪育种和人类生物医学研究的重要潜在候选基因。

结  果

猪泛基因组和遗传变异数据集构建

为了捕获比参考基因组Sscrofa11.1更全面的猪基因组序列，我们首先对一头亚洲野猪（AW）、三个中国品种（BMA：巴马香、TT：藏猪和LT：蓝塘猪）和三个欧洲品种（LW：大白猪、BKS：巴克夏猪和PTR：皮特兰猪）进行了PacBio高保真（HiFi）测序（图1A，S1）。获得了平均覆盖率为52.8×高精度长读长，每个样本的碱基准确度 > 99.9%（表S1）。随后，我们通过基于高通量染色体构象捕获（Hi-C）的组装生成了七个染色体水平的组装体，平均长度为2.68Gb（表S2,3），略大于Sscrofa11.1（2.68 Gb vs. 2.50 Gb）。转座子（TE）占每个基因组的38.8%，范围从35.8%到40.5%（表S4）。值得注意的是，82.33%的端粒序列和所有着丝粒重复区域都被预测出来。着丝粒重复区域的平均GC含量高于杜洛克猪基因组（49.5% vs. 41.4%），平均重复含量超过整个基因组（65.7% vs. 38.8%）（表S5,6）。与杜洛克猪的参考基因组相比，根据得到的组装统计数据和成对共线模式，这七个基因组具有更高的完整性和一致性。这些基因组将成为未来猪遗传学和基因组学研究的丰富资源（图1B、1C、S2、S3和表S7-9）。我们进一步整合了两个染色体水平的基因组，即Sscrofa11.1参考基因组和我们研究组之前报道的陆川基因组，以构建高分辨率的全基因组（表S8）。共获得134.24 Mb非冗余非参考序列，N50为287,722 bp（表S10）。134.24 Mb非冗余序列包括100.52 Mb来自中国地方品种（BMA、TT、LT和陆川），16.43 Mb来自AW，17.29 Mb来自欧洲品种猪（LW、BKS：Berkshire、PTR）。我们观察到新序列中的GC含量高于参考基因组（49.5% vs. 41.4%）（图1D）。约83.72%的非冗余非参考序列包含重复元素，高于参考基因组。基于新构建的全基因组，共检测到1,099个新的蛋白质编码基因，编码序列（CDS）N50大小为774。其中，94.54%（1,039个基因）有转录本证据支持，77.53%（852个基因）成功注释了 ≥ 1个功能项（表S11）。Sscrofa11.1组装体与推定的非冗余非参考序列一起被用作后续分析的最终参考基因组。

为了更全面、更公平地了解基因组多样性，我们使用anchorwave程序将这些组装的基因组序列锚定到参考基因组上。共发现187,927个结构变异（SV），峰值大小为300 bp（图1E），这与之前的结果一致。我们观察到，我们的SV发现结果与基于Illumina数据的SV调用集之间只有16.6%的碱基对重叠。此外，本研究中鉴定的SV中有96.8%是新的，这凸显了我们的SV检测与短读测序方法相比的灵敏度。值得注意的是，我们发现长散在核元件（LINE）和短散在核元件（SINE）是驱动结构变异的两个主要因素，主要在SV内和SV边界周围富集（图1F、G）。此外，还鉴定出7,667,617个小插入和缺失（INDEL，本研究中指≤ 50 bp）（图S4）。在序列片段分布图中，偶数高于相邻的奇数（图S4）。通过利用多组织RNA测序数据，我们发现SV通常对侧翼基因的表达水平产生负面影响（图1H）。

为了了解脂肪型猪和瘦肉型猪的遗传变异库和表型差异，我们收集了1,081个个体的全基因组测序数据，其中313个是新生成的（表S12,13），主要分布在欧洲和亚洲（图1A）。我们将所有这些个体分为五类：东方家猪（EAD），包括来自53个中国品种的665头猪（约占中国猪品种的60%）和1个韩国本土品种；西方家猪（WED），包含来自25个典型西方品种的315头猪；东方野猪（EAW），包括来自亚洲的40头野猪；西方野猪（WEW），包括46头欧洲野猪；以及来自五种不同野猪物种（猪科）的15头个体，例如Sus barbatus、Sus cebifrons、Sus celebensis、Sus verrucosus和作为外群的Phacochoerus africanus。45.3 Tb的全部测序数据（表S14）用于短序列变异发现，最终鉴定出30,143,962个单核苷酸多态性（SNP）和5,496,594个INDEL（去除低质量变异后），与美国国家生物技术信息中心（NCBI）的单核苷酸多态性数据库（dbSNP）和欧洲生物信息学研究所（EBI）的猪变异数据库进行比较，结果显示鉴定出5,277,943个新的SNP。结果显示染色体末端变异明显富集，X染色体上变异程度较低（图S5A、B）。这些变异在外显子内表现出明显的损耗，但在内含子和基因间区域内表现出高度的富集（图S5C）。大多数SNP和INDEL位于基因间（59.52%和59.08%）和内含子（37.55%和38.47%）区域，而只有0.56%的SNP和0.10%的INDEL存在于编码序列中（表S15）。

图1. 猪全基因组综合图

（A）五大类别中所有品种和物种的地理位置。直方图显示每个类别的样本量。EAD，东方家猪；WED，西方家猪；EAW，东方野猪；WEW，西方野猪；OGP，外群（其他猪科）。（B）八种猪品种基因组组装的BUSCO评估结果。杜洛克基因组（DRC）是从Ensembl数据库下载的参考组装。AW：亚洲野猪；BMA：巴马香，BKS：巴克夏；LT：蓝塘；LW：大白猪；PTR：皮特兰；TT：藏猪。（C）描述皮特兰品种组装统计数据的蜗牛图。其他六个基因组显示在图S2中。（D）参考基因组和新的非冗余非参考序列中重复序列和GC含量的百分比。（E）推定结构变体的大小分布。SVs，结构变异。（F）重复元件中结构变异的富集。LINE，长散布核元件。LTR，长末端重复。SINE，短散布核元件。（G）重复元件在结构变异边界周围的分布。（H）多种组织中有无结构变异的基因之间的表达差异。FPKM，每百万映射读数每千碱基转录本的片段数。使用的统计方法是Wilcoxon检验。***p < 0.001。

群体结构和潜在选择特征

根据邻接法进行的群体结构分析显示出清晰的五个分支细分，尽管某些EAD和WED品种之间存在祖先混合（图2A）。这种推断关系与主成分分析（PCA）模式大致一致，表现出东方猪和西方猪之间的祖先混合（图2B）。通过去除混合个体，我们最终保留了953个个体进行后续分析。通过根据计算的复合似然比（CLR）选择前1%的推定选择性扫描，我们在EAD和WED组中确定了约227 Mb具有显著选择足迹。置换检验表明，与启动子和外显子相比，选择特征在内含子、基因间区域和TE中更为丰富（图2C）。此外，与长末端重复（LTR）（Jaccard指数 = 0.014，p < 0.001）相比，选择特征在LINE（Jaccard指数 = 0.020，p < 0.001）和SINE（Jaccard指数 = 0.017，p < 0.001）中更为丰富，表明LINE和SINE对选择性驯化的贡献大于LTR。

WED中的选择性扫描显示出比EAD组更高的信号和更集中的位置（图2D）。我们在EAD和WED特异性选择性扫描区域中分别发现了925个和700个候选基因（表S16,17）。我们还发现EAD和EAW之间有245个区域重叠，大小为2,570,028 bp。此外，WED和WEW之间有36个区域重叠，大小为415,361 bp。此外，我们基于另外两种方法FST和Pi重复了选择性扫描分析，以确保结果的稳健性。不同组的成对比较列于（表S18-20）。值得注意的是，东方猪和西方猪之间共享的选择足迹和基因数量有限（图2E），这得到了Pi方法的重叠结果的支持（表S21,22）。在所有1,625个假定的受选择基因中，我们发现了先前研究报告的几个候选基因（图2D），包括CASP10、ESR1、LCORL、NCAPG、NR6A1和PLAG1。

此外，我们还提供了一系列值得进一步探索的扫描基因的新观察结果（图2D，S6A），例如用于繁殖性状的MPP7和用于肉质的PRKAG3。发现许多转录因子（TF）（如PAX3、RUNX1、TP63和TBX19）都处于选择状态，这些转录因子参与了许多重要的生物过程。几个扫描区间远离具有特定功能的基因（如E2F2和MYOG），可能充当基因间或远端功能元件。基因本体（GO）富集分析显示，来自EAD的扫描基因集在发育生长和对生长因子的反应方面显著富集。相比之下，WED组中过度表达的条目主要指Wnt信号通路、心血管系统发育和细胞器定位（图S6B）。我们通过与猪的动物数量性状基因座（QTL）数据库整合，进一步探讨了扫描基因的表型结果，并为EAD和WED品种之间明显不同的选择方向提供了更多支持。超几何检验表明，EAD扫描基因驱动的性状与脂肪和繁殖类别有关，而WED组中选择的基因主要涉及生长、解剖和构象类别（图2F和表S23,24）。

图2. 东西方猪的种群结构和选择性扫描

（A）所有1,081个个体的种群结构。邻接系统发育树包括665头东方家猪（EAD）、315头西方家猪（WED）、40头东方野猪（EAW）、40头西方野猪（WEW）和外群的15个个体。圆形条形图表示每个样本的个体祖先系数，其中我们将K值设置为2。（B）基于所有假定的常染色体SNP的主成分分析。（C）通过Jaccard指数评分对假定的选择性扫描与启动子、外显子、内含子、基因间和转座因子之间的重叠分析。TE，转座因子。使用的统计方法是置换检验。**p < 0.01；***p < 0.001；ns，不显著。（D）东方猪和西方猪群中所有显著选择性扫描的曼哈顿图。图中标记了有希望的候选基因。这一结果表明西方家猪受到的选择压力更大。请注意，WED的Y轴范围与EAD不同，并且WED更大。（E）东方猪和西方猪之间共享和独特的基因组间隔以及选择基因。（F）针对猪QTL数据库对东方和西方扫描基因进行功能富集分析。

与选择性扫描相关的编码变异的功能意义

鉴于编码序列内的突变可能会影响蛋白质的结构和功能，我们首先探索了这些基因组变异的功能含义。总的来说，我们在全基因组范围内检测到了156,654个常染色体编码SNP。通过删除扫描基因之外的SNP或在EAD和WED组之间具有delta等位基因频率（∆AF）≤ 0.7的SNP，我们最终保留了990个几乎固定的编码变异（图3A）。Ensembl变异效应预测器（VEP）返回626个SNP（63.23%）为同义突变，339个（34.24%）为错义突变，19个为剪接突变，6个为终止增益和起始丢失突变（表S25）。先前的研究表明，NR6A1基因5号外显子中的错义替换（chr1_265,347,265_A/G、c.T575C或p.Pro192Leu）是影响椎骨数量的致病位点。我们为这一发现提供了令人信服的证据，并说明了氨基酸替换的保守性和功能（图S7）。

除了NR6A1基因外，我们还重点介绍了几个与重要经济性状相关的候选基因。我们发现两个位置仅在WED组中具有惊人的高选择性扫描信号（图3B），其中包含八个功能基因，包括BRCA1、NBR1和RPL27。值得注意的是，对该区域的进一步检查发现了位于BRCA1基因10号外显子的c.G965A（chr12_19,812,845_G/A）的编码突变。在此位点，突变型A等位基因在EAD组（AF = 91.67%）中占主导地位，但在WED组（AF = 3.27%）中几乎不存在（图S8A-C）。同时，根据Ensembl SIFT（0.06）和BLOSUM62（-1）得分以及多物种比对，该变异趋于轻度保守且可能有害（图S8D）。为了了解该变异的影响，我们比较了5头大白猪（欧洲，WED）公猪和16头民猪（中国，EAD）母猪杂交产生的589个F2个体中三种基因型的表型差异，其中在民猪和大白猪中，突变类型A的比例分布分别为73.33%（11/15，15头民猪的重测序数据）和0%（0/4，4头大白猪的重测序数据）。我们观察到该SNP对240天屠宰重量（SW）的等位基因替换效应为3.55公斤，对瘦肉率（LMP）的等位基因替换效应为0.62%。在这个基因座上，AA基因型表现出比GG和GA基因型明显更低的SW和LMP（图3C）。通过体外脂肪细胞增殖试验，G322D BRCA1的表达显示出增强的脂肪细胞生长和发育能力（图3D、E），与具有野生型G等位基因的WED品种相比，脂肪型EAD猪的脂肪沉积更为显著。

此外，与ABCA3基因重叠的扫描区域特别引起我们的关注，因为它在EAD和WED组之间的基因型分布非常不同，并且该基因肺特异性表达（图3F、S9A、S9B）。我们确定了三个非同义突变，它们彼此之间几乎完全LD。根据SIFT和BLOSUM62这两个值，我们优先考虑第19外显子中的变异（chr3_39,723,445_C/T），因为该变异对蛋白质结构和功能的影响更大（SIFT = 0.04, BLOSUM62 = -1）（图3G、S9C、S9D）。EAD组中几乎所有个体都携带此位点的替代等位基因（T），而WED猪的C等位基因几乎固定。由于ABCA3突变可导致肺表面活性剂缺乏、肺泡细胞损伤和炎症，我们假设ABCA3基因是导致品种特异性抗性/易感性病原体肺炎支原体（Mhp）的有趣遗传决定因素。体外Mhp感染实验数据显示，野生型C等位基因与IL1B基因表达水平较低相关（图3H），这表明该等位基因与免疫反应和肺损伤减轻有关。

图3. 选择性扫描驱动的编码变异的功能含义

（A）EAD组和WED组之间delta等位基因频率（∆AF）大于0.7的编码变异的注释。SNP，单核苷酸多态性。（B）EAD、EAW、WED和WEW组之间BRCA1区域的复合似然比（CLR）值和基因型模式。虚线紫色和红色分别代表前5%和1%选择性扫描的阈值。（C）BRCA1基因的编码SNP导致240天屠宰重量和瘦肉比例存在显著差异。使用的统计方法是单因素方差分析（ANOVA）。**p < 0.01；***p < 0.001；ns，不显著。（D）野生型和突变型等位基因之间猪脂肪细胞EdU增殖测定的比较。突变型A等位基因的Edu阳性细胞明显高于野生型G等位基因，表明BRCA1中突变型A等位基因比野生型G等位基因具有更高的脂肪细胞增殖能力。数据以平均值±标准差（SD）表示。所用的统计方法是Student’s t检验。**p < 0.01。（E）通过定量逆转录聚合酶链反应（qRT-PCR）检测PCNA标志物的表达增加，表明突变等位基因的增殖能力增强。PCNA是一种重要的增殖标志基因。突变型A等位基因细胞中PCNA基因表达明显高于野生型G等位基因细胞，表明BRCA1中G到A的突变增加了脂肪细胞的增殖能力。值以平均值±SD表示。所用的统计方法是Student’s t检验。**p < 0.01。（F）EAD、EAW、WED和WEW组ABCA3区域的基因型景观。在此区域中发现了三个几乎固定的编码变体。（G）十种代表性哺乳动物物种的氨基酸多序列比对。（H）根据qRT-PCR检测的IL1B标记物表达降低，野生型等位基因的免疫损伤减轻。IL1B是评估免疫损伤的常用标记基因。免疫损伤与IL1B表达水平呈正相关。值以平均值±SD表示。使用的统计方法是Student’s t检验。**p < 0.01。

扫描基因调节EAD猪和WED猪骨骼肌的发育异时性

扫描基因与猪QTL数据库的富集分析表明，EAD组和WED组之间的生长和胴体特征存在很大差异（图2F）。为了研究EAD和WED猪之间生长和胴体特征差异逐渐形成的分子机制，我们以通城猪（TC，东方脂肪型品种）和长白猪（LDR，西方瘦肉型品种）为研究模型，比较了两个品种27个阶段的骨骼肌发育轨迹。表型记录显示，这两个品种的生长速度和肌肉质量差异很大，160日龄时LDR猪的体重是TC猪的两倍多（图4A）。因此，我们推断长白猪和通城猪的骨骼肌生长发育应该是不同步的。随后，我们对比了TC和LDR骨骼肌27个时间点的RNA-seq数据，包括15个产前阶段和12个产后阶段（表S26）。综合分析表明，TC和LDR之间骨骼肌发育的转录分析具有全局一致性（图S10,11）。随后，我们应用加权基因共表达网络分析（WGCNA）来探索与骨骼肌发育有关的关键功能模块。共发现18个模块（图S12A），其中两个与假定的扫描基因有显著重叠，主要与细胞分裂、形态发生和发育过程有关（图S12B、C）。

此外，我们基于TimeMeter中的动态时间规整（DTW）算法进行了比较转录组分析，以剖析TC和LDR之间骨骼肌发育过程中的时间偏移模式。我们提供了骨骼肌生成过程中轻度异时性的线索，其特点是LDR在产前阶段成肌细胞增殖更先进、更持久，产后肌肉生长更快，而TC在产前阶段成肌细胞成熟速度更快（图4B，S13A）。总的来说，我们在产前和产后阶段分别检测到了4,947个和1,184个差异进展基因（DPG），包括MYOG、TP53INP2和MEF2B（图4C，S13B和表S27,28）。根据计算出的进展优先分数（PAS），正值表示TC中发育进展的表达模式比LDR中更先进，定义为TC加速基因，而负值表示LDR中加速变化（LDR加速基因）。功能富集分析强调，产前阶段的LDR加速基因在细胞周期和DNA复制中显著过度表达。相反，TC加速基因主要促进细胞分化和成熟（图4D）。在出生后时期，LDR进展更先进的基因持续富集于细胞发育和生物合成过程，而TC先进基因与免疫反应有关。

为了进一步阐明品种间发育时间差异的遗传基础，我们评估了扫描基因在调节骨骼肌生成中的贡献。与猪基因组中的所有基因相比，选择基因在每个品种的所有阶段都表现出更高的表达水平（图S14），这意味着它们在特定的重要生物过程中发挥着重要作用。后续的置换检验报告了DPG和扫描基因之间存在高度显着的重叠，尤其是在产前阶段（图4E）。我们对这些扫描基因进行了相同的DTW分析，并发现了与所有表达基因的结果高度一致的动态模式（图4F）。

图4. 通过扫描基因研究猪骨骼肌发育的品种间异时性

（A）长白猪（LDR）和通城猪（TC）在多个发育阶段的体重差异。（B）LDR猪和TC猪骨骼肌发育阶段的对应关系。（C）产前和产后阶段差异进展基因的示例。根据TimeMeter的定义，正的进展高级分数（PAS）表示TC中的更快表达模式（TC加速基因），而负值表示LDR中的加速变化（LDR加速基因）。（D）在产前和产后阶段使用LDR和TC加速基因预测的GO富集分析。（E）差异进展基因（DPG）和扫描基因之间存在显著重叠。DPG集是根据LDR和TC猪的时间过程RNA测序数据预测的。Y轴是通过1,000次置换检验确定的观察到的（紫色菱形）和预期的（箱线图）重叠数。所用的统计方法是置换检验。*p < 0.05；**p < 0.01；***p < 0.001。（F）仅基于LDR和TC猪的扫描基因的发育异时性分析。使用LDR作为参考，TimeMeter在所有RNA-seq数据中实施动态时间规整算法，以根据阶段转录组相关性选择时间序列之间的最佳比对。这一发现表明，TC和LDR之间骨骼肌发育的差异可以追溯到肌肉生成的非常早期的阶段，其中LDR显示出更提前和更长的成肌细胞增殖时间。

选择性扫描通过重塑DNA甲基化模式调节基因表达

DNA甲基化是调节真核生物基因表达的主要表观遗传因素之一。为了阐明DNA甲基化变异在猪驯化和肌肉生长发育中的贡献，我们对27个发育阶段的骨骼肌进行了全基因组亚硫酸盐测序（WGBS）（表S29），与上述具有RNA-seq数据的样本相对应。我们在TC和LDR中生成了全局DNA甲基化组景观，发现TC中的CpG甲基化水平更高（图S15、S16A-C）。此外，基于WGBS的主成分分析（PCA）揭示了与RNA-seq结果相似的分类模式，但显示出更强的品种效应（图S16D）。甲基化组和转录组之间的PC1和PC2均表现出强烈的曲线关系（r_PC1 = 0.858，r_PC2 = 0.873），表明发育轨迹和品种的分类清晰（图5A）。这两类数据集之间的功能相关性再次证实了启动子和基因体区域的结果（图S17-20）。通过成对差异甲基化分析，我们获得了166,265个差异甲基化区域（DMR）（表S30），与LDR相比，TC中约65%的差异甲基化区域为低甲基化。DMR与五个基因组特征的交集显示，高甲基化和低甲基化区域主要存在于内含子和TE区域中（图5B）。UCSC基因组浏览器中的可视化展示了具有不同分布模式的五个重要DMR基因座（图5C）。

为了了解选择性扫描对DNA甲基化的影响，我们分别基于这些扫描基因、所有表达基因和随机采样基因对WGBS和RNA-seq数据进行了相关性分析。当将仅基于这些扫描基因的结果与基于所有表达基因和随机采样基因的结果进行比较时，我们发现启动子的调控关系降低，相关系数降低（LDR为0.08和0.09，TC为0.14和0.14），但基因体的模式增强，相关系数增加（LDR为0.08和0.09，TC为0.09和0.10）（图S21）。这一发现强调了由于DNA甲基化变化，自然和人工选择下的序列环境可以重塑整个发育过程中的整体表达模式。与所有表达的基因相比，推定的扫描基因的启动子区域在TC和LDR品种中通常具有较低的DNA甲基化水平，而基因体总体上显示出较高的DNA甲基化水平（图5D、S22、S23）。根据DMR域与最近的SNP之间的绝对距离，我们发现大多数DMR（98.39%）位于已识别SNP的50kb范围内，约55%的DMR位于5 kb范围内（图5E）。此外，结果表明，这些几乎固定的SNP（∆AF大于0.7）更容易位于DMR附近，而∆AF低于0.7的SNP则不然。为了进一步阐明DNA甲基化在猪驯化中的重要性及其与遗传选择的关系，我们首先根据发育时间点将DMR分为三组，分别在产前和产后阶段存在显著差异。我们发现，在几乎所有阶段（产前11 ~ 15个阶段和产后9 ~ 12个阶段）都存在持续DNA甲基化变异的类别比其他两组具有更高比例的Hypo-DMR域（图5F，上），支持了我们上面提到的结果。通过置换检验确定的富集分析表明，具有持续变化的DMR容易在产前和产后时期被选择性扫描所覆盖（图5F，下图）。

图5. 猪驯化和肌肉生长发育过程中DNA甲基化组的差异

（A）PC1和PC2的甲基化组与转录组的相关性分析。PC1和PC2值分别从RNA-seq和全基因组亚硫酸盐测序（WGBS）数据中计算得出。（B）启动子、外显子、内含子、TE和基因间区域等五种不同基因组特征中差异甲基化区域（DMR）比例的堆叠图。（C）UCSC基因组浏览器展示图，显示六个代表性阶段中五个显着的差异甲基化区域（DMR）位点，选择性扫描区域被着灰色。E，产前阶段；D，产后阶段。这表明TC和LDR猪骨骼肌发育中甲基化水平的差异受到驯化和选择的影响。（D）EAD选择性基因、WED选择性基因和所有表达基因之间的DNA甲基化水平比较。（E）全基因组水平上DMR与最近的SNP之间的绝对距离。（F）出生前和出生后骨骼肌发育过程中DNA甲基化变异与选择性扫描之间的关系。在第二幅图中：我们首先从基因组中提取区域，并计算与DMR重叠且属于扫描区域的区域数。然后，我们计算DMR和扫描区域之间的重叠区域数，并使用置换检验对结果进行显著性检验。倍数变化表示从随机抽样获得的中位数与DMR和扫描区域之间重叠区域的实际数量之间的比率。使用的统计方法是置换检验。*p < 0.05；**p < 0.01；***p < 0.001。

骨骼肌发育和产肉性状候选基因的遗传基础

我们首先探索了结构变异对候选基因的调控作用。一个有趣的候选基因是BDH1基因，它主要参与骨骼肌、肝脏和心脏组织中酮代谢过程的调节。BDH1基因在出生后阶段的TC骨骼肌中表现出比在LDR中更快的表达模式（图6A），我们在东方猪品种中发现了位于该基因上游的272 bp插入。该插入被具有强信号的选择性扫描覆盖（图6B）。根据3D相互作用热图（图6B），我们推测该突变通过潜在的染色质相互作用对BDH1基因发挥了正调控作用。荧光素酶报告基因检测表明该插入的增强子活性增加，定量逆转录聚合酶链反应（qRT-PCR）实验显示中国地方猪品种的BDH1表达水平高于LDR品种（图6C、D）。进一步的过表达/抑制分析显示该基因使肌肉增殖能力增强，但成肌细胞分化率较弱（图6E、F、S24）。

图6. 与骨骼肌发育相关的BDH1基因综合分析

（A）TimeMeter分析显示BDH1在TC出生后阶段的加速表达。（B）四组中的选择性扫描区域。紫色虚线和红色虚线分别表示BDH1区域中前5%和1% CLR值的阈值。（C）HEK293T细胞中的荧光素酶报告基因检测，以比较插入之间的增强子活性。值以平均值±SD表示。使用的统计方法是Student's t检验。**p < 0.01。（D）通过qRT-PCR检测中国地方猪品种中BDH1的基因表达较高。分别将中国地方猪品种中BDH1的表达与LDR中的表达进行比较。值以平均值 ± SD表示。使用的统计方法是Student's t检验。***p < 0.001。（E）通过qRT-PCR分析在C2C12细胞中评估BDH1过表达后PCNA、Ki67、MyoG和MyHC标志物的表达水平的细胞增殖。值以平均值±SD表示。使用的统计方法是Student’s t检验。*p < 0.05；**p < 0.01。（F）BDH1过表达后小鼠C2C12细胞的增殖和分化测定。值以平均值 ± SD表示。使用的统计方法是Student's t检验。**p < 0.01。

接下来，我们重点研究了在假定的选择性扫描中从甲基化CpG到TpG/CpA的CpG-SNP突变的功能后果，并探讨了它们对异时性基因和骨骼肌发育的调控功能。此外，我们发现胚胎期GHSR基因在LDR中的进展（PAS = -0.814）比TC更先进（图7A）。我们选择该基因主要基于以下五点：1）该基因是扫描基因。2）该基因在LDR和TC之间有DMR。3）GHSR中的突变位点靠近DMR。4）EAD和WED组之间GHSR中SNP的等位基因频率差异（ΔAF）高于0.7。5）该基因与生长性状和产肉性状有关。这一结果表明GHSR应该是加速成肌细胞增殖的关键因素。多个实验支持GHSR可以促进成肌细胞增殖，但抑制其在体内和体外的分化和融合（图7B、C、S25）。同时，在GHSR的内含子中检测到明显的选择信号，其中包含两个显著的Hypo-DMR（图7D）。对该区域的进一步分析表明，在EAD群体和TC品种中，突变等位基因的CpG-SNP（chr13_111,051,076_C/T）几乎完全固定，由于从C到T等位基因的转变，TC中的甲基化状态完全消失（图7E）。我们通过基序扫描分析发现，SNP通过替代的T等位基因创建了TBX21的新转录因子结合位点。随后的荧光素酶报告基因检测显示T等位基因具有更高的增强子活性，这可能是延迟TC中GHSR表达的关键因素（图7F）。该SNP对D240屠宰重量的估计等位基因替代效应为4.22公斤，携带一个或两个显性T等位基因拷贝的个体的屠宰重量明显较低（图7G）。基于公开的猪骨骼肌Hi-C数据，我们发现CpG-SNP可能调节GHSR基因的表达，因为两者都位于同一拓扑相关域（TAD）中（图7H）。此外，基于类似策略，我们还提出异时候选基因CD36和NCAPG-LCORL簇是猪骨骼肌发育和屠宰重量的关键因素（图S25-28）。

图7. 与骨骼肌发育和产肉性状相关的GHSR基因的综合分析

（A）TimeMeter分析显示LDR产前阶段GHSR的加速表达。（B）GHSR过表达后小鼠C2C12细胞的增殖和分化测定。值以平均值 ± SD表示。使用的统计方法是Student's t检验。**p < 0.01。（C）分别在体内注射AAV介导的非靶向对照和pcDNA3.1-GHSR后10天小鼠骨骼肌的H&E染色和肌纤维数量。代表性图像以20倍放大率显示（比例尺 = 100 μm）。使用的统计方法是Student's t检验。*p < 0.05。（D）选择性扫描重塑LDR和TC猪六个代表性发育阶段的DNA甲基化水平。紫色和红色虚线分别表示GHSR区域中前5%和1% CLR值的阈值。（E）GHSR中固定的C到T替换导致DNA甲基化水平的去除和TBX21新转录因子结合位点的产生。（F）HEK293T细胞中的荧光素酶报告基因检测以比较两个等位基因之间的增强子活性。值表示为平均值±SD。使用的统计方法是Student's t检验。***p < 0.001。（G）CC基因型在第240天的屠宰重量更高。使用的统计方法是单向方差分析。***p < 0.001。（H）推定的拓扑相关域（TAD）表明GHSR基因座和候选CpG-SNP之间存在潜在相互作用。

讨  论

长期的人工选择促成了许多具有惊人表型多样性的猪品种的形成。在过去的几十年里，现代育种方法对杜洛克猪、长白猪和大白猪等猪品种施加了巨大的选择压力，以改善其理想性状。这些商业瘦肉型品种主要被塑造为胴体产量高、生长速度快，而东方猪品种（主要是中国猪）则以肉质优良、脂肪含量高和成熟期早而闻名。然而，表型变异背后的遗传结构，尤其是对形态发生的影响，尚未完全了解。我们基于九个高质量基因组和1,081个个体进行了全基因组分析，这是迄今为止最广泛的数据集，使我们能够利用编码和调控特征对几个重要性状的贡献，特别是骨骼肌发育进行研究。当前的研究不仅建立了丰富的资源来探索基于基因组序列、基因表达和DNA甲基化数据的功能变异，还可以通过重塑猪的基因组、转录组和表观基因组来增进我们对人工选择如何决定表型变化的机制基础的理解。

越来越多的研究探索了猪在驯化和繁殖过程中的全基因组选择性扫描特征，但这些研究中的品种数量和样本量仍然有限。凭借这里介绍的9个组装基因组和1,081个个体（包括78个不同的家养品种），我们提供了全面的全基因组和遗传变异资源，并大大扩展了当前遗传信息的目录。我们的工作为东西方猪种群之间的表型和遗传分化提供了进一步的支持。西方几个品种的密集育种实践巩固了更集中的区域，而东方猪的选择性扫描显示出整个基因组中更分散的模式，这可能承受的选择压力较小。此外，EAD猪和WED猪仅共享极少数的扫掠区域和基因，支持这两类猪相对独立的驯化和育种过程。丰富的潜在候选基因库为相关研究和未来猪的基因组设计育种提供了宝贵的资源。

所有假定遗传变异的富集/耗竭分析均显示外显子内显著耗竭，表明对编码多态性（尤其是具有潜在有害影响的错义突变）的纯化选择信号强烈。然而，许多编码位点被有意识地保留下来，并通过现代育种程序进一步分散在种群中，受益于它们对若干经济利益性状的有利影响，如MC1R、NR6A1、PRKAG3、RYR1和SERPINA6基因中的致病突变。除了这些众所周知的基因外，我们还通过生物信息学和实验方法确定了两个有希望的候选基因BRCA1和ABCA3。BRCA1基因编码一种肿瘤抑制因子，该基因在人类中的突变与因细胞周期途径异常而导致的乳腺癌和卵巢癌易感性有关。我们观察到WED组中BRCA1基因选择性扫描的信号非常强烈。有趣的是，BRCA1也被证明在鸡的驯化过程中受到选择，特别是在商业肉鸡养殖中。这一一致的发现反映了人类驱动的选择过程中对农场动物特定优选育种目标的趋同模式，例如更高的肉类产量和更快的肌肉生长。猪地方性肺炎（PEP）是一种慢性高度传染性呼吸道疾病，造成相当大的经济损失，主要由Mhp引起。一般来说，中国地方品种对Mhp的敏感性高于西方瘦肉型品种。我们的结果表明，突变的T等位基因导致非同义替换，从而导致肺泡巨噬细胞减少和宿主防御进一步减弱，这得到了体外细胞培养中的Mhp感染实验的支持。此外，我们还发现了几个编码突变的基因，如ALS2、EXOC2、PRDM15和SEMA3B/E，它们在肌肉生长、脂肪沉积和神经元发育中起着至关重要的作用。总之，这些功能性编码变体为未来的养猪育种提供了宝贵的信息。

不同的人工选择导致中国本土和西方瘦肉型猪品种的产肉量存在巨大差异，但这种差异背后的发育程序尚不清楚。我们推断，表型结果归因于肌肉生成过程的时间。异时性被广泛定义为生物体与其祖先或其他生物体相比，发育事件的速率或时间在基因控制下的变化，这可能导致大小、形状和某些器官特征的巨大差异。多项研究探索了分子异时性，主要涉及多个物种的大脑和视网膜，并确定了许多驱动有机发育、形态发生和进化的枢纽基因。人们普遍认为肌肉生长可分为两个主要时期：产前肌肉生成，由肌纤维数量的增加（增生）决定；产后生长，通过增加肌纤维的大小（肥大）实现。尽管相似性很高，但我们的研究结果表明，LDR在产前阶段的肌纤维形成时间比TC相对较长，但在产后阶段肌肉生长更快，瘦肉增加更多。我们发现产前阶段有更多差异进展基因，与推定的候选基因有更大的重叠，这意味着人工选择对骨骼肌发育的影响可以追溯到肌肉生成的非常早期的阶段。综上所述，我们的数据表明，脂肪型猪和瘦肉型猪种之间猪肉产量的巨大差异归因于骨骼肌发育的异时性，人工选择是肌肉生长的关键驱动力，尤其是对于产前肌肉的发育。

DNA甲基化是一种经常出现在CpG二核苷酸上的表观遗传标记，在塑造动物驯化表型多样性方面起着重要作用。与以前的研究结果一致，我们发现长白猪比通城猪显示出更多的高甲基化区域，这是由于广泛的繁殖造成的。选择性扫描似乎更有可能与在多个阶段检测到的DMR重叠，这意味着周围的固定等位基因可能对DNA甲基化变异产生持续影响。此外，异时性基因，尤其是那些被选择性扫描覆盖的基因，容易受到DMR的调节，这表明具有差异甲基化水平的位点对于控制发育进程至关重要。这一发现为理解人工选择如何通过重塑表观基因组来塑造器官发生的时空模式提供了新的见解。

总体而言，甲基化胞嘧啶表现出较高的突变率，许多研究报告称CpG到TpG/CpA的转变与发育障碍之间存在显著关联。大多数受选择的遗传基础位于非编码区内，通常作为调控元件来控制基因表达。我们发现一个插入调节东方猪BDH1的增强子，这导致东方和西方品种的肌纤维增殖和分化存在差异。BDH1是一种重要的限速酶，负责酮代谢和ATP合成。插入诱导的BDH1表达增高可能在肌纤维增殖和分化中发挥关键作用。除了结构变异之外，本研究表明，（几乎）固定的SNP是品种间甲基化模式差异的重要来源，可以作为与感兴趣性状相关的调控序列变体。GHSR是经过筛选的强有力候选基因，通常与牲畜的生长和胴体性状相关，我们对调控SNP（chr13_111,051,076_C/T）的关联结果也支持了这一点。此外，妊娠期间GHSR和生长素释放肽浓度较高会导致胎儿体重增加和出生后体重增加。有趣的是，我们发现GHSR基因仅在LDR猪的产前阶段表现出更高级的表达模式，这可能是DNA甲基化状态改变的结果。这一发现意味着生长素释放肽主要发生在妊娠早期，并且在LDR中持续时间更长，从而导致肌肉质量和体重增加。这些数据证明了猪作为生物医学模型在探索脂肪型和糖尿病等某些人类代谢疾病的分子机制方面的重要性。

尽管我们的研究比较全面，但仍存在一些局限性。首先，依赖有限数量的基因组（尽管质量很高）可能无法完全捕捉所有猪品种中存在的遗传多样性。其次，虽然我们整合了多个组学数据，但基因调控的复杂性意味着其他因素（例如环境影响和DNA甲基化以外的表观遗传修饰）也可能在表型分化中发挥重要作用，但本研究并未详细阐述。此外，已确定的调控元件的功能验证及其对性状的影响仍有待在未来的研究中探索。

结  论

本研究阐明了猪的复杂性状和表型分化的遗传基础，揭示了非编码区和特定基因（如BDH1和GHSR）的关键作用。我们发现了显著的结构变异和SNP，这些变异和SNP导致了感兴趣的性状的品种间差异，特别是在肌纤维增殖和分化方面。基因组、转录组和甲基组数据的整合突出了这些表型变异背后的调控机制。这些发现不仅增强了我们对猪遗传学的理解，还为旨在优化牲畜性状的遗传改良计划提供了宝贵的资源。

方  法

样本采集和测序

为了捕捉世界各地猪的比较完整的遗传多样性，我们首先收集了30个不同品种的313个个体的耳组织，这些品种代表了中国气候和地理极端的种群。使用标准苯酚-氯仿法从这313个样本的耳组织中提取基因组DNA用于二代重测序。我们的合作者共享了来自11个本地品种的共98个个体的测序数据。此外，我们从SRA数据库下载了670个高深度全基因组测序数据（表S12），以建立全面的猪基因组序列数据集。

本研究新采集了通城猪（TC）27个发育阶段的81个骨骼肌（背最长肌）样本，包括胚胎第33、40、45、50、55、60、65、70、75、80、85、90、95、100和105天（缩写为E33、E40、E45、E50、E55、E60、E65、E70、E75、E80、E85、E90、E95、E100和E105）以及出生后第0、9、20、30、40、60、80、100、120、140、160和180天（缩写为D0、D9、D20、D30、D40、D60、D80、D100、D120、D140、D160和D180）。进一步提取了这81个样本的总RNA和DNA，并在Illumina HiSeq X Ten平台（Illumina，CA）上进行WGBS和RNA-seq测序，每个发育阶段的三个生物学重复。长白猪（LDR）的转录组（81个样本）和甲基化组（81个样本）数据来自我们团队之前的研究。最终，我们使用162个转录组和162个DNA甲基化数据来研究27个阶段的骨骼肌生长发育。值得注意的是，长白猪和通城猪在同一个猪场饲养，饲养环境、饲养饲料和管理模式相同。

基因组组装

使用PacBio Sequel平台和循环共识测序（CCS）模式对七个猪品种（一个AW、三个中国地方猪品种（BMA、TT和LT）和三个西方猪品种（LW、BKS和PTR））进行了测序。利用Hifiasm（v0.16.1_r375）软件使用默认参数从HiFi CCS reads生成组装。从相同样本生成插入片段大小为300-700 bp的Hi-C片段文库，并在Illumina平台上测序。Hi-C文库中使用的酶是DpnII，它在“GATC”处切割DNA。使用fastp（v0.19.5）程序修剪原始reads的接头序列，并去除低质量的双端reads以获得干净的reads。Bowtie2（v2.2.9）用于将干净的读段与组装的重叠群对齐。我们使用HiC-Pro流程（v3.1.0）和默认参数过滤了低质量读段。使用Juicer（v1.6）和3D-DNA流程（v180419）将有效读段用于锚定染色体。使用包含哺乳动物odb10基因集（9,226个BUSCO基因）的Benchmarking Universal Single-Copy Orthologs（BUSCO）程序（v5.3.2）评估组装完整性。报告了单拷贝加重复的完整BUSCO基因计数。使用与近乎完整的绵羊基因组研究中类似的方法，应用21聚体的Merqury（v1.3）来计算质量值（QV）。

泛基因组构建

使用minimap2（-a -x asm 10）（v2.24-r1122）将七个HiFi组装体和我们之前的陆川基因组与杜洛克猪参考基因组（Sscrofa11.1）比对。可靠比对定义为长度超过300 bp且序列同一性高于90%的连续比对。没有可靠比对的序列保留为未比对序列。接下来，使用MUMmer（v4.0.0rc1）使用参数-maxmatch将未比对序列映射到Sscrofa11.1基因组，然后使用delta-filter -I 90 -l 300 -1参数和一对一比对块选项来过滤比对结果。使用blastn（参数为-word_size 20 -max_hsps 1 -max_target_seqs 1 -dust no -soft_masking false -evalue 0.00001）（v2.9.0）将所得序列与Sscrofa11.1基因组比对，并删除与Sscrofa11.1序列同一性 > 90%的序列。将剩余序列根据200 bp内的相邻区域合并，并删除长度 < 500bp的序列。接下来，使用BLASTN（v2.9.0）将生成的序列与GenBanknt数据库比对，参数为“-evalue 1e-5 -best_hit_overhang 0.25 -perc_identity 0.5 -max_target_seqs 1”。具有来自其他物种的最佳匹配或已知动物线粒体基因组覆盖的序列可能是污染物，因此被移除。随后，将从所有组装中获得的剩余序列组合在一起。为了消除冗余，使用minimap2（v2.24-r1122）（-X -a -x asm 10）进行了全对全比对。

转座子注释

转座子的识别使用了两种方法：从头重复识别和在现有数据库中搜索已知重复。RepeatModeler（v2.0.1）用于预测基因组中的重复序列，然后使用RepeatMasker（v4.0.7）在从头转座因子（TE）库中搜索基因组。RepeatMasker（v4.0.7）和Repbase（v15.02）数据库用于识别组装基因组中已知的转座因子TE重复序列。RepeatMasker用于DNA水平识别，RepeatProteinMasker用于进行蛋白质水平识别。我们在这里使用的注释是Ensembl release-94。我们采用了与长白猪研究中描述的类似的方法来识别七个新创建的基因组中的端粒。根据之前的分析，获得了预期染色体位置中的着丝粒重复区域。如果无法通过同源比对确定着丝粒区域，我们采用quartTeT（v1.1.8）方法进行识别。

基因预测

对于同源预测方法，从Ensembl数据库（Ensembl版本105）下载了六种哺乳动物基因组（包括人类、小鼠、牛、狗、山羊和杜洛克猪）的基因组序列和注释文件。此外，从中国国家基因库（CNGB；https://db.cngb.org/）下载了陆川猪蛋白质序列，登录号为CNP0001159。使用GeMoMa（v1.8）将这些序列与基因组组装进行比对，以检测同源肽。对于基于RNA测序的预测方法，使用了九种组织的原始Illumina短读段以及一个RNA库（即皮下脂肪、肾脏、心脏、肺、背最长肌、肝脏、腰大肌、脾脏、卵巢）。登录号：SRR3160015、SRR3160012、SRR3160008、SRR3160011、SRR3160014、SRR3160009、SRR3160017、SRR3160013、SRR3160010、SRR3160016）从NCBI下载以进行进一步分析。所有原始读段均使用fastp（v0.19.5）和默认设置进行评估。使用HISAT2（v2.0.1）将干净的读段与基因组组装对齐，以识别假定的外显子区域和剪接点。然后使用StringTie（v1.2.2）将映射的读段组装到基因模型中，并通过拼接对齐程序（PASA）（v2.4.1）进行验证。具有正确结构但同源性预测缺失的PASA支持基因被添加到基因集中。最后，使用PASA确定非翻译区和可变剪接区。

基因功能注释

使用五个公共数据库对参考基因集进行功能注释，包括NCBI非冗余蛋白质序列数据库、SwissProt、京都基因和基因组百科全书（KEGG，84.0）、欧洲分子生物学实验室翻译和基因本体论（GO）。使用InterProScan确定这些基因的推定域和GO术语，而使用Diamond程序（v0.9.30）将猪基因组的蛋白质序列与其余四个公共数据库进行比较，E值截止值为1e-05。

SV检测

我们采用了与Munasinghe、Read等人中描述的类似方法来检测SV。简而言之，使用anchorwave （v1.0.1）通过“genoAli”命令和“-IV”参数将八个猪基因组（AW、BMA、BKS、LT、LC、LW、PTR、TT）中的每一个与杜洛克猪（Sscrofa11.1）进行比较。使用maf-convert.py（https://gitlab.com/mcfrith/last/-/blob/main/bin/maf-convert）将输出的.maf文件转换为.sam格式。最后，应用samtools mpileup（v1.10）来调用SV。

将全基因组测序数据比对到泛基因组

我们将假定的非冗余非参考序列和Sscrofa11.1参考组装合并为猪泛基因组。首先使用fastp（v0.19.5）使用默认参数修剪fastq数据。接下来，使用BWA-MEM管道（v0.7.17-r1188）将所有干净reads（包括我们新生成的样本）与泛基因组对齐。然后对成功比对的reads进行排序，并使用Picard工具（v2.23.3）和SAMtools（v1.10）删除重复项。

全基因组SNP和INDEL鉴定

使用GATK UnifiedGenotyper（v3.7-0-gcfedb67）和-glm BOTH -rf BadCigar --sample_ploidy 2选项对每个样本进行全基因组变异检测。为确保变异检测的准确性，筛选了QD < 2.0 || FS > 60.0 || MQ < 20.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0的SNP。使用ANNOVAR（v2013-06-21）包和Ensembl Variant Effect Predictor包，根据全基因组注释进行基于基因的SNP注释。根据基因组注释，SNP被分为外显子区域、5’或3’非翻译区、内含子区域、剪接位点（剪接点2bp以内）、上游和下游区域（转录起始位点上游或下游1kb区域内）或基因间区域。编码外显子中的SNP进一步分为同义SNP或非同义SNP。为了评估错义突变的序列保守性，我们首先从NCBI或Ensembl数据库下载了九种相关哺乳动物的氨基酸序列，包括人类、小鼠、牛、马、山羊、虎、猫、狗和兔。随后，使用MUSCLE程序进行多重序列比对。同时，使用SWISS-MODEL预测野生型和突变型的蛋白质结构模型。

系统发育和种群遗传分析

为了深入了解不同猪品种之间的系统发育关系，我们根据所有常染色体高质量双等位基因SNP（call rate ≥ 90%且次要等位基因频率 ≥ 5%）进行了全面的基因组研究。使用TreeBeST（http://treesoft.sourceforge.net/treebest.shtml）构建了邻接树，并进行了200次引导重复。使用交互式生命之树（iTOL）（v6）显示该树。为了推断种群结构，我们使用了ADMIXTURE（v1.3.0），它实现了块松弛算法。我们还通过测试Hardy-Weinberg平衡（HWE）违规（p值 > 10^-4）来过滤SNP，并重建基于模型的聚类分析。为了确定最佳遗传簇K，对从2到10的每个K值测试了交叉验证误差。PCA使用程序GTAC（v1.92）进行。

选择性扫描和分化分析

我们首先根据种群结构的结果排除东方猪和西方猪之间的潜在混合样本，以确保假定选择性扫描的准确性。当K = 2时，如果西方猪的血统在东方品种中超过20%，它将被过滤掉，西方品种则反之亦然。从此分析中去除了次要等位基因频率低于5%的SNP。随后，使用计算先进的复合似然比（CLR）检验通过SweeD软件（v4.0.0）检测所有种群的全基因组选择性扫描。计算所有常染色体上每个10 kb非重叠窗口的CLR分数。使用Pi比率（π野生/π驯化，对西方和东方种群分别进行测试）检测野生种群和驯化种群之间潜在的选择性扫描，并使用VCFtools（v0.1.13）以50 kb的滑动窗口和10 kb的步长获得东方和西方猪之间的FST值。包含少于10个SNP的窗口被排除在进一步分析之外。前1%值的窗口被视为候选区域。对于前1%的CLR分数，将删除驯化和野生种群之间重叠的窗口。将上述每种方法的结果结合起来研究猪的驯化和分化。候选基因和区域的连锁不平衡（LD）块大小由LDBlockShow（v1.30）计算。

为了评估候选变异与猪的经济重要性状之间的关联，我们利用了合作伙伴先前研究中589个F2个体的基因型和表型数据，这些数据来自5头欧洲大白猪（WED）和16头中国民猪（EAD）杂交。我们首先使用自定义Python脚本根据出生顺序对谱系信息进行排序，然后建立相对亲属系数矩阵。最后，我们使用矩阵作为协变量的混合线性模型（MLM）比较了240天屠宰重量和瘦肉率的差异。使用Bonferroni方法校正p值，校正后的p值小于0.05被视为显着差异。为了详细分类假定的选择性扫描的功能类别，我们估算了这些扫描区域与5个基因组特征（启动子、外显子、内含子、基因间和TE）之间的Jaccard指数，这些特征有1,000个排列。使用Metascape工具对扫描的基因进行GO富集分析。校正p值 < 0.05的GO术语被认为是显著富集的。为了探索由假定的选择特征驱动的富集表型，我们通过超几何检验对猪QTL数据库进行了性状/QTL富集分析。我们重点关注置信区间小于1 Mb的QTL，因为QTL置信区间太大，无法在后处理中有效使用。

转录组测序

用TRIzol试剂（Invitrogen，上海，中国）提取总RNA。每个样本的RNA总量为3 μg，用作RNA样本制备的材料。用Epicentre Ribo-zeroTM rRNA去除试剂盒（Epicentre，威斯康星州）去除核糖体RNA，用乙醇沉淀法扫描无rRNA残留物。用Qubit荧光计（Invitrogen）监测RNA浓度，在文库制备前用Agilent Bioanalyzer 2100系统（Agilent Technologies，加利福尼亚州）评估RNA质量。随后，使用NEBNext UltraTM Directional RNA Library Prep Kit for Illumina（NEB，马萨诸塞州）生成测序文库。最后在Illumina NovaSeq 6000平台（Illumina，CA）生成双端150 bp测序reads。

全基因组重亚硫酸盐测序

将含有0.5 ng lambda DNA的100 ng基因组DNA用Covaris S220超声破碎成200 – 300 bp。使用EZDNA甲基化GoldTM试剂盒（Zymo Research）对这些DNA片段进行重亚硫酸盐处理，并由诺禾致源公司（北京，中国）构建文库。随后，在Agilent Bioanalyzer 2100系统上评估文库质量，并在Illumina Novaseq 6000平台（Illumina，CA）上对每个样本进行150 bp的双端测序。

RNA-seq分析

原始数据通过fastp （v0.19.5）以默认参数进行过滤。使用STAR（v2.7.1）将过滤得到的reads比对到构建的全基因组和注释文件中。同时计算基因表达水平，并将其用作DESeq2包（v1.38.3）的输入。基于DESeq2软件中rlogTransformation函数对这些标准化计数进行PCA、层次聚类和k均值聚类。随后，我们应用TimeMeter工具（v1.0.5）来评估时间基因表达相似性并检测LDR和TC品种之间的差异进展基因。PAS表示TC中的表达模式更快（TC加速基因），而负数表示LDR中的变化加速（LDR加速基因）。使用Metascape进行差异表达基因（DEG）预测的生物过程分析的GO术语分析。

全基因组重亚硫酸盐测序分析

首先使用Bismark工具将我们构建的基因组转换为重亚硫酸盐基因组。然后根据默认参数使用Bismark（v0.22.0）对质量控制后的reads比对。过滤重复reads后提取每个胞嘧啶位点的甲基化信息。使用MethylKit软件（v0.99.2）识别差异甲基化胞嘧啶位点（DMC）和DMR。对于DMC分析，去除覆盖率小于10X的胞嘧啶位点。采用滑动窗口法计算DMR，窗口和步长均为1000bp，去除覆盖度小于5的胞嘧啶位点，保留至少3个胞嘧啶位点的区域用于下游分析。DMC和DMR均以Bonferroni校正q值 <  0.05和meth.diff > 30为标准定义。我们重点关注基因体区域（从转录起始位点（TSS）到转录终止位点（TES））、启动子区域（TSS上游和下游2kb）和转座子元件来评估甲基化信息。

基因表达与DNA甲基化的相关性分析

相关性分析前，我们剔除rlog-normalizedcounts小于0和大于20的基因，然后计算每个表达基因的基因体和启动子区域的DNA甲基化水平，在R软件（v4.2.1）中计算不同特征的DNA甲基化水平与基因表达之间的Pearson相关系数（r）。

Motif扫描分析

使用Homer软件（v4.11）进行Motifs扫描分析，使用findMotifsGenome.pl函数在默认参数下进行，以SNPs及其上下游10 kb为区域制作bed文件，然后根据基因组中的序列进行motif扫描分析。

体外脂肪细胞增殖试验

在无菌条件下收集新生猪颈部和背部的皮下脂肪组织，用含有高浓度青霉素和链霉素的PBS缓冲液洗涤三次。去除可见血管及结缔组织，将脂肪组织切成1 mm³的小块，用I型胶原酶溶液（800 U/ml）在37 ℃下消化1 h，每5 min搅拌一次，脂肪组织消化物经双层尼龙筛（100 μm、25 μm）过滤，滤液1500 r/min离心10 min，获得猪脂肪细胞初步培养物，将脂肪细胞悬浮于无血清DMEM/F12（DMEM/F12 -Dulbecco's Modified Eagle Medium）培养基中，1000 r/min再次离心5 min，除去残留消化物。最后，将脂肪细胞沉淀物悬浮在完全培养基（DMEM/F12 90% + 致死牛血清（FBS，Gibco，NY）10%+青霉素-链霉素（PS，Thermo Scientific，MA）中，并接种到T25培养瓶中进一步生长。对于脂肪细胞的增殖分析，采用了两种主要方法：EDU incorporation和qRT-PCR。EDU（5-乙炔基-2'-脱氧尿苷）测定法用于测量DNA合成，其作为细胞增殖的指标。除了EDU，还采用qRT-PCR来量化增殖相关基因的表达水平。

细胞培养

小鼠C2C12成肌细胞和HEK-293T细胞来自北京协和医院，并在添加了10% FBS和1% 青霉素-链霉素（PS, ThermoScientific, MA）的杜氏改良Eagle培养基（DMEM, Gibco, NY）中培养。此外，猪脂肪细胞在添加了10% FBS和1% 青霉素-链霉素（PS, Thermo Scientific, MA）的DMEM/F12中培养。所有细胞均在37 °C的5% CO2培养箱中培养。具体而言，当细胞达到约80-90%时，用分化培养基（添加了2%马血清的DMEM）替换生长培养基来诱导C2C12成肌细胞分化。在4-6天内每48小时更换一次培养基。为了监测分化效果，我们使用了形态学和分子标记。形态学变化通过细长、多核肌管的出现来评估，这是成功分化的关键指标。在分子水平上，我们通过qRT-PCR监测肌生成素和肌球蛋白重链（MHC）等成肌标记物的表达，这些标记物被广泛认为是C2C12分化的标记物。

质粒构建

通过PCR扩增小鼠CD36、GHSR和NCAPG基因的全长编码序列区。总反应体积为50 µL，包含2.5 ul模板DNA、2.5 ul含有BamHI和XhoI位点的正向和反向引物、25 ul 2 × Vazyme Lamp Master Mix（Vazyme，中国）和17.5 ul无核酸酶水。PCR程序包括初始变性（95 °C，5分钟），然后进行35个扩增循环（95 °C，30秒（变性）、56 °C，30秒（退火）和72 °C，1分钟（延伸）。反应以72 °C，7分钟的最终延伸结束，以确保PCR产物完全合成，然后使用琼脂糖凝胶电泳分析PCR产物以确认扩增片段的大小正确和存在。

PCR扩增后，使用BamHI和XhoI（美国赛默飞世尔科技）对产物进行双消化。消化在总反应体积为50 µL，其中包含10 µL PCR产物（约0.5 - 1 µg）、与两种酶兼容的1X反应缓冲液、每种酶1 µL（10单位）和无核酸酶水。反应液按照厂家说明书37 ℃反应2 h，酶切产物用MolPure® PCR纯化试剂盒（YESEN，中国）纯化后进行克隆。酶切产物与载体（pcDNA3.1（+）载体（Invitrogen，中国））纯化后，按照载体：插入片段的摩尔比1:3进行连接反应，总反应体积为20 µL，其中载体DNA约50 ~ 100ng，插入片段DNA按摩尔比计算，10X T4 DNA连接酶缓冲液（含ATP）2 µL，T4 DNA连接酶（Thermo Fisher Scientific，美国）1 µL（400单位）。反应液4 ℃反应过夜，连接后产物转化大肠杆菌进行后续筛选。

为了验证CD36、GHSR、NCAPG和BDH1基因中候选SNP和插入的双荧光素酶活性，首先扩增以SNP和插入为中心的400 bp DNA片段，然后插入PGL4.23双荧光素酶表达载体（Promega，威斯康星州麦迪逊）。引物采用Primer Premier 5设计，列于表S31中。

RNA干扰

根据NCBI提供的目标基因的CDS区域，使用siDirect网站（http://sidirect2.rnai.jp/）设计目标基因的siRNA，并在GenePharma（中国上海）合成。siRNA序列列于表S32中。

细胞转染和双荧光素酶报告基因检测

细胞培养12小时后，根据制造商的方案，使用Lipofectamine 3000和Opti-MEM将适当的质粒或寡核苷酸转染到C2C12成肌细胞或HEK293T细胞中。转染期为24小时后收集这些细胞。使用双荧光素酶报告基因检测系统（Promega）按照制造商的说明量化荧光素酶活性。将荧光素酶活性标准化为海肾荧光素酶活性。

定量RT-PCR（qRT-PCR）

对于CD36、BRCA1、GHSR和NCAPG的表达分析，从培养细胞中分离总RNA。在BDH1表达分析的情况下，从猪的不同骨骼肌（LDR、TC和BM）以及培养细胞中分离总RNA。使用TRIzol试剂（Invitrogen，CA）从培养细胞中分离总RNA，然后根据制造商的方案使用PrimeScript™ RT Master Mix（Perfect Real Time）（Takara，北京，中国）逆转录为cDNA。根据SYBR Premix Ex TaqTM说明，在StepONE Real-Time PCR系统（Applied Biosystems）上进行qRT-PCR。每个反应在20 μl反应混合物中进行，其中包含10 μl ChamQ SYBR qPCR Master Mix（Vazyme Q311-02，南京，中国）、0.5 μl基因特异性引物、2 μl模板cDNA和7 μl无菌水。qRT-PCR的实施条件如下：初始变性95 °C 10分钟，40个循环（变性95 °C 15秒，退火60 °C 30秒，延伸72 °C 30秒），然后进行熔解曲线分析。所有反应使用来自三个独立个体的cDNA重复三次，结果使用2^-ΔΔCT方法分析，以管家基因为标准。

5-乙炔基-2'-脱氧尿苷（EDU）测定

将C2C12成肌细胞和脂肪细胞接种于12孔板中。当细胞生长至50%汇合密度时，用过表达质粒、siRNA或miRNA模拟物转染。转染48小时后，将成肌细胞暴露于50 μ MEDU（RiboBio，中国）中，37 °C下2小时。随后，将细胞在4％多聚甲醛中固定30分钟，用2 mg/ml甘氨酸溶液中和，然后通过添加0.5％TritonX-100进行通透。加入含有EDU（Apollo Reaction Cocktail；RiboBio，中国）的溶液，并将细胞在室温下孵育30分钟。然后加入核染料Hoechst 33342，继续孵育30分钟。使用荧光显微镜（DMi8，Leica，德国）捕捉三个随机选择的视野以观察EDU染色细胞的数量。

免疫荧光分析

将培养细胞用多聚甲醛固定在6孔板中，用0.5％ Triton处理，用山羊血清封闭1小时，并与抗MyHC（1：500，DSHB MF20，上海，中国）抗体孵育2小时。接下来，将细胞与山羊抗小鼠二抗孵育。最后，加入DAPI（1：1000，Invitrogen D3571，上海，中国），并使用Leica DMI3000 B显微镜（Leica）观察细胞。

AAV介导的小鼠体内GHSR和CD36过表达

为了在体内验证GHSR和CD36基因的功能，我们首先从HANBI（上海，中国）获得了腺相关病毒9（AAV9）血清型pcDNA3.1-GHSR、pcDNA3.1-CD36或pcDNA3.1-Control。同时，从华福康公司（北京，中国）获得了7周龄小鼠。每只小鼠的右腿胫骨前肌（TA）肌肉注射带有pcDNA3.1-GHSR或pcDNA3.1-CD36的AAV9病毒（100 µL，滴度 ≥ 1×10¹³ vg/ml），左侧TA肌肉注射pcDNA3.1-Control。用异氟烷对AAV9载体处理的小鼠进行麻醉，并通过颈椎脱位处死以收集肌肉。对于GHSR，我们在注射AAV9 10天后收集TA肌肉。为了验证CD36调节加速肌肉成熟，我们设计了一个肌肉损伤和再生模型。我们首先在注射心脏毒素（CTX）前10天注射过表达CD36或对照的AAV9病毒，然后在注射CTX后第0、1、5和10天收集TA肌肉。在最佳切割温度包埋后，将冷冻切片在4%多聚甲醛中固定过夜，并根据苏木精-伊红/HE染色试剂盒（Solarbio，北京，中国）进行H&E染色。

肺炎支原体攻毒实验

将猪肺泡巨噬细胞（PAMs）的细胞密度调整为1×10⁴细胞/孔，置于24孔板中，在0.5 ml RPMI 1640培养基中培养，培养基中添加10% FBS和2%青霉素-链霉素溶液，37 ℃培养。培养12小时后，用PBS清洗细胞，以去除任何未附着的细胞，然后根据制造商的说明，使用QIAGEN转染试剂（货号301005）转染0.4 μg DNA。转染24小时后，再次转染细胞。随后，用肺炎支原体处理具有野生型和突变型等位基因的PAMs。孵育6小时后，收获这些感染细胞进行qRT-PCR。PAMs和猪肺炎支原体材料来自江苏省农业科学院兽医研究所。

数据统计分析

使用R（v4.2.1）和GraphPad Prism（v9）进行统计分析。进行置换检验以评估选择特征在不同基因组区域的富集情况。使用置换检验进行富集分析以确定DMR与产前和产后持续变化和选择性扫描之间的关联。使用R中的“cor”函数和“Pearson”方法对基因表达水平和DNA甲基化模式进行Pearson相关性分析。用于多种组织中有无结构变异的基因表达差异的统计方法是Wilcoxon检验。对于D240时的瘦肉比例和屠宰重量，使用单因素方差分析（ANOVA）在多个组之间进行统计比较。对于细胞测定、qRT-PCR结果、EDU和免疫荧光，使用Student’s t检验确定两组之间的差异显着性。当*p < 0.05、**p < 0.01和***p < 0.001时，具有统计学意义。

（▼ 点击跳转）

“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的高水平综合期刊，欢迎投稿！