基因組和DNA甲基化組聯合分析篩選豬肉質性狀關鍵基因

趙真堅，王凱，陳棟，申琦，余楊，崔晟頔，王俊戈，陳子旸，禹世欣，陳佳苗，王翔楓，唐國慶

基因組和DNA甲基化組聯合分析篩選豬肉質性狀關鍵基因

趙真堅，王凱，陳棟，申琦，余楊，崔晟頔，王俊戈，陳子旸，禹世欣，陳佳苗，王翔楓，唐國慶

四川農業大學動物科技學院/農業農村部畜禽生物組學重點實驗室/畜禽遺傳資源發掘與創新利用四川省重點實驗室/豬禽種業全國重點實驗室，成都 611130

【背景】豬的肉質性狀是重要的經濟性狀。研究影響各類肉質性狀的分子機制，發掘關鍵基因，從而指導豬的遺傳改良，對改善豬肉品質具有重要意義。當前肉質相關的機制研究主要是基于DNA的基因組研究，而對肉質性狀的DNA甲基化研究以及結合基因組和甲基化組的綜合分析卻鮮有報道。【目的】通過基因組和DNA甲基化組聯合分析，篩選、鑒定了影響豬肉質的潛在關鍵基因。為豬肉品質的遺傳改良研究提供借鑒。【方法】檢測了140頭大白豬背最長肌的28個肉質性狀，通過表觀基因組關聯分析（EWAS）與全基因組關聯分析（GWAS）篩選了各性狀顯著關聯的CpG和SNP位點。隨后以SNP為協變量對GWAS和EWAS重疊的顯著關聯位點進行條件關聯分析，進一步篩選具有獨立效應的CpG位點。然后以CpG位點的甲基化水平為因變量，以SNP為自變量進行關聯分析從而鑒定甲基化數量性狀位點（meQTL）。最后使用順式甲基化數量性狀位點（cis-meQTL）作為工具變量進行孟德爾隨機化分析，從而推斷cis-meQTL與表型之間的因果關系，同時對位點進行注釋，鑒定潛在的關鍵基因。【結果】（1）在屠宰45 min黃度值（b45min）、滴水損失（DL）、二十二碳六烯酸（C22:6n-3）3個肉質性狀上，EWAS和GWAS在相同基因組區域鑒定到顯著關聯位點。（2）b45min的7個CpG位點在條件關聯分析后仍保持顯著，DL有1個CpG位點在條件關聯分析后仍保持顯著，而C22:6n-3的3個位點在使用SNP作為協變量分析后不再顯著，表明EWAS鑒定的b45min的7個CpG位點和DL的1個CpG位點的顯著關聯不受附近的顯著SNP影響。（3）b45min的7個CpG位點和DL的1個CpG位點共鑒定了10個meQTL，但絕大多數是trans-meQTL，只有一個CpG位點（SSC12:44 254 675 bp）鑒定到一個cis-meQTL，表明該位點可能受到近距離SNP調控。（4）孟德爾隨機化分析顯示該CpG位點（SSC12:44 254 675 bp）與b45min表型存在一定的因果關聯。（5）對該位點注釋發現，距離CpG位點（SSC12:44 254 675 bp）和其cis-meQTL最近的基因是NOS2，且CpG位點位于NOS2基因內。【結論】綜合DNA甲基化組合基因組數據聯合分析結果，可以推測NOS2基因是肉色性狀關鍵候選基因，其DNA甲基化、SNP共同作用調控基因表達，進而影響肉色性狀相關基因表達。

表觀基因組關聯分析；全基因組關聯分析；肉質性狀；DNA甲基化；甲基化數量性狀位點

0 引言

【研究意義】豬的肉質性狀是重要的經濟性狀，是肉品的營養學、衛生學、組織學和加工學特性的總和[1]。肉質性狀既包含了客觀特征如pH值、肉色、系水力等，也包括主觀特征如風味、口感等[2]。豬肉的肉質好壞直接決定豬肉在市場的接受程度。因此，研究影響各類肉質性狀的分子機制，發掘關鍵基因，從而指導豬的遺傳改良，對改善豬肉品質具有重要意義。【前人研究進展】國內外眾多研究者對影響肉質性狀的基因展開大量研究，發現了一些對肉質有明顯影響的主效基因，如氟烷基因[3]（halothane，HAL）、酸肉基因[4]（rendement napole，RN）、RKAG3[5]（AMP-activated protein kinase γ 3）等。此外也鑒定了與豬肉質性狀相關的重要候選基因，如瘦素基因[6-7]、脂肪酸結合蛋白基因家族[8]（fatty acid binding protein，FABP）、FASN[9-10]編碼脂肪酸合成酶、MYOD[11-13]基因家族等。研究者通過全基因組關聯分析（genome- wide association study，GWAS）對肉質性狀的遺傳機制進行解析[14]。Gao等[15]對582頭商業豬中對pH值、肉色、大理石花紋、導電率等肉質性狀進行GWAS分析，研究發現一些關鍵的候選基因如與pH值相關的PPP1R3B基因、與肉色相關的SLC15A4等。Liu等[16]對杜洛克和二花臉豬F2群體進行GWAS分析，研究發現了GALNT15、GALNTL2和HTATIP2等與pH值和滴水損失顯著相關的候選基因。然而，豬的肉質性狀是由多基因控制的復雜性狀，單一組學數據不能全面揭示影響肉質性狀的分子機制。近些年來，表觀遺傳學和甲基化測序技術的發展為畜禽重要經濟性狀的表觀遺傳研究提供基礎，例如豬[17]、牛[18]、雞[19]、綿羊[20]和山羊[21]等。目前，對豬肉質性狀的DNA甲基化研究主要集中于不同品種不同組織的全基因組DNA甲基化差異的比較分析。唐韶青等[22]發現，在同一個物種的不同組織中，其DNA甲基化模式存在特異性，其中組織的甲基化水平一般要高于血液。肌肉組織和脂肪組織是研究肉質性狀最重要的組織。LI等[23]基于長白豬、藏豬和榮昌豬2個部位的肌肉組織和8個部位脂肪組織的甲基化測序數據，描繪了豬肌肉和脂肪組織的甲基化圖譜，揭示了不同品種間和不同組織間的差異甲基化位點，并證實了甲基化在肥胖的發生過程中起重要作用。張小麗[24]對3頭長白豬不同部位的脂肪組織進行全基因組DNA甲基化分析，研究發現基因外顯子區有更高的甲基化程度，而基因啟動子區域的甲基化會抑制基因的轉錄。總體而言，研究者對豬肉質性狀的研究主要集中于基于DNA的基因組研究，而對肉質性狀的DNA甲基化研究以及結合基因組和甲基化組的綜合分析卻鮮有報道。【本研究切入點】多項研究表明遺傳變異，通常為單核苷酸多態性（single nucleotide polymorphism, SNP），可以對DNA甲基化產生重大影響[25]，這些變異位點即為甲基化數量性狀位點（methylation quantitative trait loci，meQTL）。通過meQTL可以探究復雜性狀潛在的生物學機制，而順式meQTL（cis-meQTL）的識別是大多數研究的熱點。【擬解決的關鍵問題】基于純種大白豬背最長肌的DNA甲基化組和基因組數據，利用表觀基因組關聯分析（epigenome-wide association analysis, EWAS）和全基因組關聯分析鑒定了與肉質性狀相關的關鍵位點。隨后對GWAS鑒定的顯著SNPs與EWAS鑒定的顯著CpG位點進行比較分析，通過條件關聯分析篩選具有獨立顯著信號的CpG位點，再進行meQTL鑒定和孟德爾隨機化驗證，最后鑒定了影響肉質的潛在關鍵基因。本研究為應用甲基化組和基因組聯合分析，篩選肉質性狀關鍵基因提供了新的見解。

1 材料與方法

1.1 試驗動物與表型測定

選取來自浙江天蓬豬場的140頭純種大白豬（包括85頭公豬和55頭母豬）。所有試驗豬在相同條件及環境的圈舍中飼養，并且執行相同的飼養標準。試驗豬只屠宰時體重為（111.71±12.92）kg，在屠宰前禁食24 h。在屠宰時收集每頭豬的左半邊胴體第3—4肋處背最長肌樣品，并測定pH值、肉色、滴水損失、肌內脂肪、脂肪酸組成等肉質性狀。肉質性狀的測定參考李娜[26]、秦一禾[27]以及張俊杰[28]的學位論文。測定具體性狀包括：屠宰45 min pH值（pH45min）、屠宰24 h pH值（pH24h）、屠宰45 min亮度值（L45min）、屠宰45 min紅度值（a45min）、屠宰45 min黃度值（b45min）、屠宰24 h亮度值（L24h）、屠宰24 h紅度值（a24h）、屠宰24 h黃度值（b24h）、滴水損失（DL）、剪切力（SF）、肌內脂肪（IMF）、癸酸（C10:0）、月桂酸（C12:0）、肉蔻酸（C14:0）、棕櫚酸（C16:0）、硬脂酸（C18:0）、花生酸（C20:0）、棕櫚油酸（C16:1）、油酸（C18:1n-9）、順式-11-二十碳一烯酸（C20:1n-9）、亞油酸（C18:2n-6）、α-亞麻酸（C18:3n-3）、γ-亞麻酸（C18:3n-6）、順式-11,14-二十碳二烯酸（C20:2）、順式-11,14,17-二十碳三烯酸（C20:3n-3）、順式-8,11,14-二十碳三烯酸（C20:3n-6）、花生四烯酸（C20:4n-6）、二十二碳六烯酸（C22:6n-3）。

1.2 組織樣本采集及DNA提取

在對140頭試驗豬進行屠宰測定的同時采集背最長肌樣品于凍存管中，加入75%的乙醇，于-20 ℃條件下保存備用。樣品DNA的提取參照OMEGA組織DNA提取試劑盒說明書所述步驟進行。

1.3 簡化甲基化測序及質量控制

簡化甲基化測序（reduced representation bisulfite sequencing，RRBS）是一種準確且經濟的DNA甲基化測序方案。RRBS通過酶切MspI（methylation- insensitive restriction enzyme）富集啟動子及CpG島區域，并進行Bisulfite測序。本研究采集140頭大白豬的背最長肌樣品，基于RRBS技術進行甲基化測序。首先檢測DNA樣品質量，通過瓊脂糖凝膠電泳結果判斷DNA降解程度以及是否有RNA污染，然后利用Nanodrop檢測DNA濃度，最后利用Qubit對DNA濃度精確定量。合格的DNA樣品加入一定比例的陰性對照（lambda DNA），接著使用MspI對DNA樣品進行酶切處理（37 ℃，16 h）。隨后對酶切后的DNA片段采用末端修復、加Ａ尾，并連上胞嘧啶均經過甲基化修飾的全部測序接頭。通過切膠選擇插入片段長度在40—220 bp范圍的DNA片段進行Bisulfite處理，經處理后沒有發生甲基化的C變為U（經過PCR擴增后變為T），而甲基化的C保持不變。最后進行PCR擴增，獲得最終的DNA文庫。文庫構建完成后進行質量檢驗，使用Qubit2.0對構建的文庫進行初步定量，將文庫稀釋至1 ng·μL-1，隨后使用Agilent 2100檢測文庫的插入片段長度，達到預期片段長度后，使用qPCR方法對文庫的有效濃度進行準確定量，保證文庫質量，文庫質量檢驗合格后，進行上機測序。甲基化測序交由北京諾禾致源科技公司完成。

測序完成獲得原始下機數據raw data，對raw data進行初步過濾以去除質量差的數據，使用Trimmomatic軟件[29]對raw data過濾，去除測序數據的測序接頭和低質量片段，從而獲得高質量clean data。首先采用滑動窗口的方式去除低質量reads，4個堿基為一個窗口，若該窗口平均堿基質量低于15，則在此截去reads，參數設置為SLIDINGWINDOW: 4:15。去除尾部質量小于3或者包含無法確定的堿基信息的reads，參數設置為TRAILING:3。隨后去除接頭污染的reads，去除接頭有兩種模式，包括simple alignment mode（seed與接頭序列比對得分超過7）的模式和palindrome mode（當read1與read2的重疊區域堿基評分超過30時，截去seed部分序列）的模式，參數設置為ILLUMINACLIP:adapter.fa:2:30:7:1:true。最后去除不能形成paired的reads。

經過初步過濾后獲得clean data被用于后續分析，然后使用FastQC軟件對clean data進行再次質量統計，以保證數據質量。數據質量合格后，將clean data比對到豬參考基因組上（sus scrofa 11.1）上進行分析。使用Bismark v0.22.1軟件[30]的bismark_genome_ preparation函數準備參考基因組；調用bowtie2 v2.3.5.1軟件[31]對雙端測序文件進行比對，生成BAM文件；使用Bismark v0.22.1軟件的deduplicate bismark函數對BAM文件去除重復序列；使用Bismark v0.22.1軟件的bismark methylation extractor函數提取甲基化統計信息文件。提取的甲基化統計信息文件中包含覆蓋率信息以及在一個基因組位置上甲基化reads和未甲基化reads的比例信息。

1.4 表觀基因組關聯分析

在EWAS分析之前，首先對每個個體的甲基化統計信息文件進行質量控制以獲得用于EWAS分析的高質量位點。質量控制標準包括：在所有個體中只保留CpG位點；每個個體中只保留覆蓋度大于10 reads的位點；所有個體中保留至少在30個個體中檢出的CpG位點。過濾后，基于R軟件包CpGassoc v2.60對CpG位點甲基化水平和肉質性狀進行關聯分析。在EWAS分析時，首先使用Box-Cox轉換將表型數據校正至符合正態分布，然后用出生年、出生月、性別和批次作為固定因子對表型進行校正，獲得表型校正值，然后應用混合線性模型進行關聯分析，模型如下：

=+++

模型中，表示表型向量；表示平均值；表示CpG甲基化水平矩陣；表示CpG效應；表示隨機效應；表示殘差，服從（～(0,2)）分布，表示單位向量。

使用Bonferroni校正法來設置全基因組顯著閾值，即EWAS結果中當CpG位點的關聯值達到＜0.05/時，表示該位點為顯著位點，其中表示EWAS分析中使用的CpG位點數量。

1.5 SNP分型和基因型填充

基于Illumina平臺（Illumina，美國）標準流程使用Neogen GeneSeek Porcine50K芯片對140頭大白豬進行基因分型，芯片包含50 697個SNPs，隨后對芯片數據進行填充。基因型填充參考模板是本課題獲取的40頭大白豬重測序數據。對重測序的原始數據使用FastQC軟件進行嚴格的質量控制以獲得clean data，然后對clean data進行序列比對和SNP檢測。使用BWA v0.7.17軟件[32]將clean reads比對到豬參考基因組上（sus scrofa 11.1），設置比對參數mem -t 4 -k 32 -M；使用SAMtools v1.19軟件[33]將比對獲得的SAM格式文件轉換為BAM格式文件并進行排序；使用Picard v1.119軟件的MarkerDuplicates去除BAM文件中的PCR重復序列；使用GATK v4.1.9.0軟件[34]的HaplotypeCaller函數對BAM文件進行變異檢測，從而獲得原始SNP信息；使用GATK v4.1.9.0軟件對原始SNP信息文件進行質量控制以去除結果中的假陽性位點，過濾標準為：QualByDepth （QD）＜2.0，FisherStrand（FS）＜60.0，RMSMappingQuality（MQ）＜40.0，MappingQualityRankSumTest＜-12.5以及ReadPosRankSumTest＜-8.0。

填充完成后利用 Plink v1.09軟件[35]對填充后的SNP數據進行質量控制以獲得高質量SNP位點，過濾標準為：去除實際等位基因與最可能等位基因相關性平方（R2）小于0.9的SNP位點；去除MAF小于0.01的SNP位點；去除HWE檢驗小于1×10-6的SNP位點；去除性染色體上的SNP位點。經過質量控制的SNP位點被用于后續GWAS分析。

1.6 全基因組關聯分析

使用GEMMA v0.98.5軟件[36]中的線性混合模型對28個肉質性狀進行GWAS分析。在GWAS分析前，先使用Box-Cox變換將表型數據校正至符合正態分布，然后以出生年、出生月、性別、批次作為固定因子對表型進行校正，得到的校正值被用于GWAS分析。GWAS所使用的線性混合模型如下：

=++

模型中，表示表型值向量；表示SNP效應；表示剩余多基因效應，服從（～(0,2)）分布，表示分子親緣關系矩陣（genomic relatedness matrix，GRM）；表示殘差，服從（～(0,2)）分布，表示單位向量；和分別表示和的關聯矩陣。GRM根據如下公式計算：

使用Bonferroni校正法來設置全基因組顯著閾值，即GWAS結果中當SNP的關聯值達到＜0.05/時，表示該SNP為顯著SNP位點，其中表示GWAS分析中使用的SNP數量；建議顯著閾值設置為＜1/。使用R軟件中GenABEL包[37]中estlambda函數計算基因組膨脹系數λ用于判斷是否有假陽性信號。

1.7 鑒定EWAS與GWAS重疊的顯著關聯位點

本研究首先對140頭試驗大白豬進行50K芯片基因分型，經質量控制和基因型填充后進行GWAS分析。同時基于RRBS對所有個體進行了甲基化測序，經過質量控制之后用于EWAS分析。相關結果發表在意大利動物科學雜志[38]。本文對肉質性狀的GWAS鑒定的顯著SNPs與EWAS鑒定的顯著CpG位點進行比較，從而鑒定GWAS與EWAS重疊的顯著關聯位點。如果GWAS的顯著SNPs與EWAS顯著CpG位點在2 Mb范圍內，則被認為是二者重疊的顯著關聯位點[39]。

1.8 條件關聯分析

為了確定GWAS與EWAS重疊的顯著關聯位點是由SNPs和CpG共同導致或者是由同一位置的SNPs或CpG的獨立信號導致，基于GWAS與EWAS重疊的顯著關聯，通過條件關聯分析，使用重疊的顯著關聯位點的SNP作為協變量重新進行EWAS分析[39]。模型如下：

=++++

模型中，表示表型向量；表示平均值；表示CpG甲基化水平矩陣；表示CpG效應；表示SNP取值；表示SNP效應；表示隨機效應；表示殘差，服從（～(0,2)）分布，表示單位向量。

使用Bonferroni校正法來設置全基因組顯著閾值，即EWAS結果中當CpG位點的關聯值達到＜0.05/時，表示該位點為顯著位點，其中表示EWAS分析中使用的CpG位點數量。如果使用SNP作為協變量重新進行EWAS分析后，原本顯著的CpG位點的顯著性消失，表明該EWAS中CpG位點的顯著性受到GWAS結果的混淆；如果重新進行EWAS分析后，該CpG位點仍保持顯著性，表明該重疊的顯著關聯位點是由CpG位點的顯著信號產生的，EWAS鑒定的顯著CpG位點沒有受到GWAS的混淆。

1.9 鑒定甲基化數量性狀位點

以SNP為自變量，以條件關聯分析中仍保持顯著的CpG位點的甲基化水平作為因變量并進行關聯分析，從而鑒定甲基化數量性狀位點（meQTL）。MeQTL可以根據遺傳變異與CpG位點的距離而分為順式meQTL（cis-meQTL）和反式meQTL（trans-meQTL）。其中，cis-meQTL是靠近CpG位點的遺傳變異（二者距離在1 Mb范圍內），而trans-meQTL是與CpG位點距離較遠（二者距離超過1 Mb）或者位于不同染色體的遺傳變異[40]。在關聯分析中使用混合線性模型，具體模型如下：

=++

模型中，表示CpG甲基化水平；表示基因型向量；表示SNP效應；表示剩余多基因效應的關系矩陣；表示剩余多基因效應；表示殘差，服從（～(0,2)）分布，表示單位向量。

使用Bonferroni校正法來設置全基因組顯著閾值，即關聯結果中當SNP的值達到＜0.05/時，表示該SNP為顯著SNP位點，其中表示關聯分析中使用的SNP數量，即全基因組顯著閾值為＜0.05/10961097。

1.10 孟德爾隨機化分析

孟德爾隨機化是指個體的遺傳變異導致不同程度的DNA甲基化或不同水平的特定表型[41]。使用單樣本孟德爾隨機化分析來將-meQTL與EWAS整合分析，可以推斷-meQTL的位點和表型之間是否存在因果關系。使用R軟件的MendelianRandomization v0.5.0包中的逆加權方差方法（INVERSE-variance weighted，IVW）對-meQTL的位點進行孟德爾隨機化分析。其中SNP與表型性狀關聯的效應值(β)和標準誤(se(β))以及SNP與CpG位點關聯的效應值(β)和標準誤(se(β))由關聯分析結果獲得。

2 結果

2.1 GWAS和EWAS重疊的顯著關聯位點

為了從基因組水平驗證EWAS結果，首先對EWAS顯著CpG位點和GWAS顯著SNP位點進行比較分析。如果GWAS的顯著SNP位點與EWAS顯著CpG位點距離在2 Mb范圍內，則被認為是GWAS與EWAS重疊的顯著關聯位點[38]。基于肉質性狀的EWAS和GWAS結果，在28個肉質性狀中，共在3個性狀（b45min、DL和C22:6n-3）上發現了GWAS與EWAS重疊的顯著關聯位點（表1）。其中，在b45min中發現了8個顯著CpG位點和59個顯著SNP位點位置接近；在DL發現1個CpG位點和1個SNP位點位置接近；在C22:6n-3中發現了3個CpG位點與146個SNP位點位置接近。這些位點被確定為重疊的顯著關聯位點。

表1 GWAS和EWAS重疊的顯著關聯信號統計

2.2 條件關聯分析

為了確定GWAS與EWAS重疊的顯著關聯位點是由SNPs和CpG共同導致或者是由同一位置的SNPs或CpG的獨立信號導致，使用 SNP 作為協變量構建混合線性模型對其中的 CpG 重新進行條件關聯分析。條件關聯分析結果顯示：對于 b45min性狀，在8個顯著的CpG位點中，有7個CpG位點在經過條件關聯分析后仍保持顯著（＜6.25×10?4），表明這7個CpG位點的顯著不受附近的 GWAS 顯著SNP 的影響；對于DL性狀，1個CpG位點在經過條件關聯分析后仍達到顯著閾值（＜0.05），表明該位點的顯著關聯不受附近的顯著SNP影響；對于C22:6n-3性狀，3個CpG位點在經過條件關聯分析后未能達到顯著閾值（＞1.67×10?2），表明這些CpG位點的顯著可能是受到GWAS顯著SNP影響導致的假陽性（表2）。

2.3 鑒定甲基化數量性狀位點

隨后使用條件關聯分析中仍保持顯著的CpG位點的甲基化水平作為因變量，以SNP作為自變量進行全基因組關聯分析鑒定meQTL，即使用b45min性狀的7個CpG位點（SSC1:4 667 485 bp、SSC7:10 489 621 bp、SSC10:36 329 538 bp、SSC12:44 254 675 bp、SSC14: 17 859 513 bp、SSC18:48 429 425 bp和SSC18:50 979 305 bp）和DL性狀的1個CpG位點（SSC7:120 954 182 bp）的甲基化水平作為因變量進行關聯分析從而鑒定meQTL（表3）。

在b45min性狀的7個CpG位點中，在其中3個位點中（SSC1:4 667 485 bp位點（圖1-A）、SSC10: 36 329 538 bp（圖1-C）和SSC18:48 429 425 bp（圖1-F））沒有識別到顯著關聯。對于SSC7:10 489 621 bp（圖1-B）位點，鑒定了2個meQTL，均為trans-meQTL，對應9個顯著SNP位點。對于SSC12:44 254 675 bp（圖1-D）位點，鑒定了2個meQTL，包括1個cis-meQTL和1個trans-meQTL，對應11個顯著SNP位點。對于SSC14:17 859 513 bp（圖1-E）位點，鑒定了2個meQTL，2個均為trans-meQTL，對應46個顯著SNP位點。對于SSC18:50 979 305 bp（圖1-G）位點，GWAS鑒定了3個meQTL， 3個均為trans-meQTL，對應345個顯著SNP位點。對于DL性狀的SSC7: 120 954 182 bp（圖1-H）位點中，鑒定了1個trans- meQTL ，對應4個顯著SNP位點。在鑒定的meQTL中，只有一個cis-meQTL，其SNP位點（SSC12: 43 262 418 bp）位于CpG位點（SSC12:44 254 675 bp）上游約1 Mb的區域，表明該位點可能受到近距離 SNP 調控。該cis-meQTL的SNP基因型與CpG位點甲基化水平如圖2-A。

表2 對GWAS與EWAS共同的顯著信號條件關聯分析結果

value (after correct)是指條件關聯分析的值

value (after correct) refers to thevalue of conditional association analysis

表3 甲基化數量性狀位點鑒定結果

2.4 孟德爾隨機化分析

為了進一步估計DNA甲基化在SNP和表型性狀關聯的中介作用，我們對鑒定的cis-meQTL的SNP位點（SSC12:43 262 418 bp）和對應CpG位點（SSC12:44 254 675 bp）進行孟德爾隨機化分析。使用cis-meQTL作為工具變量，對應CpG位點（SSC12:44 254 675 bp）作為介質變量，表型性狀（b45min）作為結果變量進行孟德爾隨機化分析，發現SNP對CpG位點甲基化的影響與SNP對表型性狀的影響存在一個線性關系（圖2-B），表明該CpG位點與屠宰45 min黃度值（b45min）存在一定的因果關聯。

隨后對該位點注釋發現，距離 CpG 位點（SSC12: 44 254 675 bp）和其cis-meQTL 最近的基因是 NOS2，且CpG位點位于 NOS2 基因內。綜合DNA甲基化組和基因組數據聯合分析結果，本研究推測NOS2 基因可能是肉色性狀關鍵候選基因，其 DNA 甲基化、SNP 共同作用調控基因表達，進而影響肉色性狀相關基因表達。

圖1 CpG位點與基因組SNP關聯分析的曼哈頓圖

3 討論

本研究通過整合DNA甲基化組和全基因組數據對肉質性狀進行EWAS分析和GWAS分析。EWAS鑒定了肉質性狀與DNA甲基化之間的3 107個顯著關聯，GWAS鑒定了肉質性狀與SNP之間的514個顯著關聯。將EWAS結果與GWAS結果聯合分析發現，在屠宰45 min黃度值（b45min）、滴水損失（DL）和二十二碳六烯酸（C22:6n-3）上EWAS和GWAS都在同一基因組區域鑒定到顯著關聯，即EWAS鑒定的顯著CpG位點和GWAS鑒定的顯著SNP位點在基因組的位置彼此接近。因此本研究通過篩選GWAS的顯著SNP位點與EWAS顯著CpG位點距離在2 Mb范圍內重疊的顯著關聯位點進行進一步分析。經過重疊信號的比較分析后，在b45min中發現了8個顯著CpG位點和59個顯著SNP位點位置接近；在DL發現1個CpG位點和1個SNP位點位置接近；在C22:6n-3中發現了3個CpG位點與146個SNP位點位置接近。這些位點被確定為重疊的顯著關聯位點，也是本研究進行下一步甲基化數量性狀位點鑒定的基礎。

A. CpG位點SSC12:44 254 675 bp的cis-meQTL的箱線圖，橫坐標表示SNP基因型，縱坐標表示CpG甲基化水平；B. cis-meQTL的孟德爾隨機化分析的散點圖，橫坐標表示SNP位點（SSC12:43 262 418 bp）對CpG 位點（SSC12:44 254 675 bp）甲基化的影響，縱坐標表示SNP位點（SSC12: 43 262 418 bp）對表型性狀（b45min）的影響

3.1 條件關聯分析

基因表達受多種調控方式影響，DNA甲基化和SNP都可以影響基因的表達，有許多研究報道了SNP、甲基化和基因表達的關系[25-42]。因此，一些EWAS的顯著信號可能是由附近的SNP導致的，由SNP或者其他因素影響引起的這種混淆可能是EWAS鑒定到大量與性狀相關的顯著DNA甲基化的主要驅動因素。這一結果也在許多人類研究中被報道。例如，SHEN等[38]通過對人類疾病蛋白質生物標志物進行GWAS和EWAS并綜合比較分析后，發現與GWAS重疊的EWAS顯著信號是由潛在的遺傳變異以及環境因素引起的。因此，本文通過以SNP為協變量對這些EWAS與GWAS重疊的顯著關聯位點進行條件關聯分析發現，在b45min、DL和C22:6n-3這3個性狀中，b45min的7個CpG位點在條件關聯分析后仍保持顯著，DL的1個CpG位點在條件關聯分析后仍保持顯著，表明這些CpG位點關聯的顯著性是獨立的，不受GWAS顯著SNP的影響；而C22:6n-3的3個位點在使用SNP作為協變量分析后不再顯著，表明這些CpG位點的關聯顯著性可能受到附近SNP的影響而產生的假陽性。

3.2 孟德爾隨機化分析

EWAS旨在系統地將整個表觀基因組中的DNA甲基化水平的變異與表型變異關聯起來。多項研究顯示遺傳變異（通常為SNPs）可以對DNA甲基化產生重大影響[25]，這些位點即為甲基化數量性狀位點（methylation quantitative trait loci，MeQTL）。MeQTL可以根據遺傳變異與CpG位點的距離而分為順式meQTL（cis-meQTL）和反式meQTL（trans- meQTL）。本文對b45min的7個CpG位點以及DL的1個CpG位點，以CpG位點的甲基化水平為因變量，以SNP為自變量進行關聯分析，從而鑒定meQTL。結果發現了10個meQTL，但絕大多數是trans-meQTL，只有一個CpG位點（SSC12:44 254 675 bp）鑒定到一個cis-meQTL，表明該位點可能受到近距離SNP調控。由于meQTL潛在的混雜作用，可能會造成DNA甲基化水平與表型之間出現顯著關聯，并且大多數研究中并沒有考慮meQTL。因此，有研究者提出在進行EWAS之前利用meQTL效應來校正DNA甲基化值以解決這個問題[43-44]，從而將meQTL整合到EWAS中。在EWAS中整合meQTL的另一種應用是通過孟德爾隨機化深入分析DNA甲基化信號與相關表型之間關聯的假定因果方向[45]。在表觀遺傳分析的背景下，將meQTL作為工具變量，DNA甲基化水平作為暴露因素，而表型性狀作為結果[46]。已經有多項研究應用孟德爾隨機化將meQTL整合到EWAS分析中[39,47-48]。此外，EWAS和meQTL信號的組合也可以用于探索遺傳-環境的相互作用[49]。因此，通過EWAS鑒定的與性狀相關的CpG位點有可能反映性狀對DNA甲基化的影響，而不是DNA甲基化對性狀的作用[48-50]。基于此，本文使用cis-meQTL作為工具變量進行孟德爾隨機化分析，結果表明了DNA甲基化對性狀的影響。

3.3 NOS2基因

距離CpG位點（SSC12:44 254 675 bp）和其cis- meQTL最近的基因是NOS2，且CpG位點位于NOS2基因內。在該基因組區域，Choi等[51]也鑒定到與肉色b值相關的QTL區域，此外，Kim等[52]也在該區域鑒定到與肉色L值相關的QTL區域。豬肌肉色素主要由肌紅蛋白、血紅蛋白以及微量有色代謝物組成。在放血充分的情況下，對肉色起主要作用的是肌紅蛋白。NOS2基因屬于一氧化氮合酶家族，一氧化氮合酶催化L-精氨酸產生NO和L-瓜氨酸[53]。Hattori等[54]研究表明在脂多糖和細胞因子刺激的血管平滑肌細胞中，PI3K-Akt通路可以通過NF-kappa B（核因子kappa B）調控iNOS的表達。Silva等[55]在小鼠上的研究發現ATP通過刺激PI3激酶和Akt來增強一氧化氮酶的活性，進而增加NO的產生。NO是一種信使分子，在許多組織、細胞和器官中的核心生物過程中作為關鍵的調節分子。NO也可以與OxyMb反應形成MetMb和硝酸鹽[56]，NO對線粒體呼吸也有影響，NO可以抑制線粒體呼吸鏈[57]，導致ATP產生減少，氧化劑產生增加，從而減少線粒體對氧氣的消耗。有研究表明，肌紅蛋白也可以通過氧化、配體結合以及形成生物活性S-亞硝基肌紅蛋白的方式調節NO平衡[58]。因此NOS2基因可能通過影響細胞內NO的濃度進而影響肉色性狀。

4 結論

利用140頭大白豬背最長肌的28個肉質性狀EWAS分析和GWAS分析的結果，篩選了GWAS的顯著SNP位點與EWAS顯著CpG位點重疊的顯著關聯位點，通過條件關聯分析篩選了具有獨立效應的顯著位點，隨后進行了meQTL鑒定和孟德爾隨機化分析。結果表明CpG位點（SSC12:44 254 675 bp）鑒定的cis-meQTL與屠宰45 min黃度值（b45min）存在一定的因果關聯，且該位點位于NOS2基因內。因此，本論文推測NOS2基因可能是肉色性狀關鍵候選基因，其DNA甲基化、SNP共同作用調控基因表達，進而影響肉色性狀相關基因表達，但其在肉色性狀調控的具體分子機制需進一步試驗研究。

[1] 孫建廣, 張石蕊. 豬肉品質研究進展. 中國豬業, 2008, 3(7): 58-59.

SUN J G, ZHANG S R. Research progress of pork quality. China Swine Industry, 2008, 3(7): 58-59. (in Chinese)

[2] 宋志芳, 解幼志, 刑荷巖, 蘆春蓮, 曹洪戰. 影響豬肉質性狀相關基因的研究進展. 中國豬業, 2017, 12(7): 48-50.

SONG Z F, XIE Y Z, XING H Y, LU C L, CAO H Z. Research progress on genes related to pork quality traits. China Swine Industry, 2017, 12(7): 48-50. (in Chinese)

[3] SHEN Q W, UNDERWOOD K R, MEANS W J, MCCORMICK R J, DU M. The halothane gene, energy metabolism, adenosine monophosphate-activated protein kinase, and glycolysis in postmortem pig longissimus dorsi muscle. Journal of Animal Science, 2007, 85(4): 1054-1061.

[4] SALAS R C D, MINGALA C N. Genetic factors affecting pork quality: halothane and rendement napole genes. Animal Biotechnology, 2017, 28(2): 148-155.

[5] 李夢云, 鄭萍, 李婉濤, 陳代文. PRKAG3基因在不同品種豬不同生長階段骨骼肌中的表達差異及其表達量與肉質的關系. 動物營養學報, 2017, 29(5): 1661-1669.

LI M Y, ZHENG P, LI W T, CHEN D W. The expression difference ofgene in skeletal muscle at different growth stages for different breeds of pigs and the relationship betweengene expression level and meat quality. Chinese Journal of Animal Nutrition, 2017, 29(5): 1661-1669. (in Chinese)

[6] MüNZBERG H. Differential leptin access into the brain-a hierarchical organization of hypothalamic leptin target sites? Physiology & Behavior, 2008, 94(5): 664-669.

[7] 于太永. 脂肪細胞因子leptin和TNF: α對豬骨骼肌成肌細胞增殖分化的影響及其機制[D]. 楊凌: 西北農林科技大學, 2007.

YU T Y. Effects of adipocytokines leptin and TNF-α on the proliferation and differentiation of porcine skeletal myoblasts and its mechanism[D]. Yangling: Northwest A & F University, 2007. (in Chinese)

[8] 李志娟. 豬心臟型脂肪酸結合蛋白研究進展. 養豬, 2021(2): 55-56.

LI Z J. Research progress of heart fatty acid-binding protein. Swine Production, 2021(2): 55-56. (in Chinese)

[9] GRZES M, SADKOWSKI S, RZEWUSKA K, SZYDLOWSKI M, SWITONSKI M. Pig fatness in relation to FASN andgenes polymorphism and their transcript level. Molecular Biology Reports, 2016, 43(5): 381-389.

[10] ZAPPATERRA M, LUISE D, ZAMBONELLI P, MELE M, SERRA A, COSTA L N, DAVOLI R. Association study between backfat fatty acid composition and SNPs in candidate genes highlights the effect ofpolymorphism in large white pigs. Meat Science, 2019, 156: 75-84.

[11] VERNER J, HUMPOLí?EK P, KNOLL A. Impact of MYOD family genes on pork traits in Large White and Landrace pigs. Journal of Animal Breeding and Genetics, 2007, 124(2): 81-85.

[12] KIM J M, CHOI B D, KIM B C, PARK S S, HONG K C. Associations of the variation in the porcine myogenin gene with muscle fibre characteristics, lean meat production and meat quality traits. Journal of Animal Breeding and Genetics = Zeitschrift Fur Tierzuchtung Und Zuchtungsbiologie, 2009, 126(2): 134-141.

[13] LEE E A, KIM J M, LIM K S, RYU Y C, JEON W M, HONG K C. Effects of variation in porcinegene on muscle fiber characteristics, lean meat production, and meat quality traits. Meat Science, 2012, 92(1): 36-43.

[14] 周李生, 楊杰, 劉先先, 張志燕, 楊斌, 麻駿武. 蘇太豬宰后72 h pH和肉色性狀的全基因組關聯分析. 中國農業科學, 2014, 47(3): 564-573. doi: 10.3864/j.issn.0578-1752.2014.03.016.

ZHOU L S, YANG J, LIU X X, ZHANG Z Y, YANG B, MA J W. Genome-wide association analyses for musle pH 72 h value and meat color traits in sutai pigs. Scientia Agricultura Sinica, 2014, 47(3): 564-573. doi: 10.3864/j.issn.0578-1752.2014.03.016.(in Chinese)

[15] GAO G X, GAO N, LI S C, KUANG W J, ZHU L, JIANG W, YU W W, GUO J B, LI Z L, YANG C Z, ZHAO Y X. Genome-wide association study of meat quality traits in a three-way crossbred commercial pig population. Frontiers in Genetics, 2021, 12: 614087.

[16] LIU X X, ZHANG J J, XIONG X W, CHEN C Y, XING Y Y, DUAN Y Y, XIAO S J, YANG B, MA J W. An integrative analysis of transcriptome and GWAS data to identify potential candidate genes influencing meat quality traits in pigs. Frontiers in Genetics, 2021, 12: 748070.

[17] YANG Y L, FAN X H, YAN J Y, CHEN M Y, ZHU M, TANG Y J, LIU S Y, TANG Z L. A comprehensive epigenome atlas reveals DNA methylation regulating skeletal muscle development. Nucleic Acids Research, 2021, 49(3): 1313-1329.

[18] HALU?KOVá J, HOLE?KOVá B, STANI?OVá J. DNA methylation studies in cattle. Journal of Applied Genetics, 2021, 62(1): 121-136.

[19] RAN J S, LI J J, YIN L Q, ZHANG D H, YU C L, DU H R, JIANG X S, YANG C W, LIU Y P. Comparative analysis of skeletal muscle DNA methylation and transcriptome of the chicken embryo at different developmental stages. Frontiers in Physiology, 2021, 12: 697121.

[20] FAN Y X, LIANG Y X, DENG K P, ZHANG Z, ZHANG G M, ZHANG Y L, WANG F. Analysis of DNA methylation profiles during sheep skeletal muscle development using whole-genome bisulfite sequencing. BMC Genomics, 2020, 21(1): 327.

[21] KANG B N, WANG J J, ZHANG H F, SHEN W, EL-MAHDY OTHMAN O, ZHAO Y, MIN L J. Genome-wide profile in DNA methylation in goat ovaries of two different litter size populations. Journal of Animal Physiology and Animal Nutrition, 2022, 106(2): 239-249.

[22] 唐韶青, 張沅, 徐青, 孫東曉, 俞英. 不同動物部分組織基因組甲基化程度的差異分析. 農業生物技術學報, 2006, 14(4): 507-510.

TANG S Q, ZHANG Y, XU Q, SUN D X, YU Y. Analysis of methylation level of genome in various tissues of different animal species. Journal of Agricultural Biotechnology, 2006, 14(4): 507-510. (in Chinese)

[23] LI M Z, WU H L, LUO Z G, XIA Y D, GUAN J Q, WANG T, GU Y R, CHEN L, ZHANG K, MA J D, et al. An atlas of DNA methylomes in porcine adipose and muscle tissues. Nature Communications, 2012, 3: 850.

[24] 張小麗. 豬背部淺層和背部深層脂肪組織全基因組甲基化研究[D]. 雅安: 四川農業大學, 2013.

ZHANG X L. Genome-wide DNA methylation changes between the superficial and deep backfat tissues of the pig[D]. Yaan: Sichuan Agricultural University, 2013. (in Chinese)

[25] VILLICA?A S, BELL J T. Genetic impacts on DNA methylation: research findings and future perspectives. Genome Biology, 2021, 22(1): 127.

[26] 李娜. 基于GWAS的豬肉品質性狀候選基因研究[D]. 北京: 中國農業大學, 2016.

LI N. Genome-wide association studies for pig meat traits and exploration of major genes[D]. Beijing: China Agricultural University, 2016. (in Chinese)

[27] 秦一禾. 豬肌內脂肪含量的全基因組關聯研究[D]. 昆明: 云南大學, 2017.

QIN Y H. Genome-wide association studies for intramuscular fat content in swine[D]. Kunming: Yunnan University, 2017. (in Chinese)

[28] 張俊杰. 利用多品種大樣本群體精細定位影響豬背最長肌脂肪酸組成的QTL[D]. 南昌: 江西農業大學, 2018.

ZHANG J J. Mapping and fine-mapping QTL affecting longissimus muscle fatty acid composition in large scale multiple pig populations[D]. Nanchang: Jiangxi Agricultural University, 2018. (in Chinese)

[29] BOLGER A M, LOHSE M, USADEL B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics, 2014, 30(15): 2114-2120.

[30] KRUEGER F, ANDREWS S R. Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics, 2011, 27(11): 1571-1572.

[31] LANGMEAD B, SALZBERG S L. Fast gapped-read alignment with Bowtie 2. Nature Methods, 2012, 9: 357-359.

[32] LI H, DURBIN R. Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics, 2009, 25(14): 1754-1760.

[33] LI H, HANDSAKER B, WYSOKER A, FENNELL T, RUAN J, HOMER N, MARTH G, ABECASIS G, DURBIN R. 1 0 0 0 Genome Project Data Processing Subgroup. The Sequence Alignment/Map format and SAMtools. Bioinformatics, 2009, 25(16): 2078-2079.

[34] DEPRISTO M A, BANKS E, POPLIN R, GARIMELLA K V, MAGUIRE J R, HARTL C, PHILIPPAKIS A A, DEL ANGEL G, RIVAS M A, HANNA M, et al. A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nature Genetics, 2011, 43: 491-498.

[35] CHANG C C, CHOW C C, TELLIER L C, VATTIKUTI S, PURCELL S M, LEE J J. Second-generation PLINK: rising to the challenge of larger and richer datasets. GigaScience, 2015, 4: 7.

[36] ZHOU X, STEPHENS M. Genome-wide efficient mixed-model analysis for association studies. Nature Genetics, 2012, 44: 821-824.

[37] AULCHENKO Y S, RIPKE S, ISAACS A, VAN DUIJN C M. GenABEL: An R library for genome-wide association analysis. Bioinformatics, 2007, 23(10): 1294-1296.

[38] SHEN Q, WANG K, WANG S J, ZHAO Z J, JI X, CHEN D, YU Y, CUI S D, WANG J G, CHEN Z Y, GU Y R, TANG G Q. Identifying key DNA methylation sites and their cis-methylation quantitative loci for intramuscular fatty acid traits using genome and methylome data in Yorkshire pigs. Italian Journal of Animal Science, 2023, 22(1): 829-840.

[39] AHSAN M, EK W E, RASK-ANDERSEN M, KARLSSON T, LIND-THOMSEN A, ENROTH S, GYLLENSTEN U, JOHANSSON ?. The relative contribution of DNA methylation and genetic variants on protein biomarkers for human diseases. PLoS Genetics, 2017, 13(9): e1007005.

[40] YU F T, QIU C, XU C, TIAN Q, ZHAO L J, WU L, DENG H W, SHEN H. Mendelian randomization identifies CpG methylation sites with mediation effects for genetic influences on BMD in peripheral blood monocytes. Frontiers in Genetics, 2020, 11: 60.

[41] BURGESS S, THOMPSON S G. Use of allele scores as instrumental variables for Mendelian randomization. International Journal of Epidemiology, 2013, 42(4): 1134-1144.

[42] ZHI D G, ASLIBEKYAN S, IRVIN M R, CLAAS S A, BORECKI I B, ORDOVAS J M, ABSHER D M, ARNETT D K. SNPs located at CpG sites modulate genome-epigenome interaction. Epigenetics, 2013, 8(8): 802-806.

[43] BIRNEY E, SMITH G D, GREALLY J M. Epigenome-wide association studies and the interpretation of disease-omics. PLoS Genetics, 2016, 12(6): e1006105.

[44] LAPPALAINEN T, GREALLY J M. Associating cellular epigenetic models with human phenotypes. Nature Reviews Genetics, 2017, 18: 441-451.

[45] DAVEY SMITH G, HEMANI G. Mendelian randomization: genetic anchors for causal inference in epidemiological studies. Human Molecular Genetics, 2014, 23(R1): R89-R98.

[46] RELTON C L, DAVEY SMITH G. Mendelian randomization: applications and limitations in epigenetic studies. Epigenomics, 2015, 7(8): 1239-1243.

[47] BATTRAM T, RICHMOND R C, BAGLIETTO L, HAYCOCK P C, PERDUCA V, BOJESEN S E, GAUNT T R, HEMANI G, GUIDA F, CARRERAS-TORRES R, et al. Appraising the causal relevance of DNA methylation for risk of lung cancer. International Journal of Epidemiology, 2019, 48(5): 1493-1504.

[48] WAHL S, DRONG A, LEHNE B, LOH M, SCOTT W R, KUNZE S, TSAI P C, RIED J S, ZHANG W H, YANG Y W, et al. Epigenome- wide association study of body mass index, and the adverse outcomes of adiposity. Nature, 2017, 541: 81-86.

[49] TSAPROUNI L G, YANG T P, BELL J, DICK K J, KANONI S, NISBET J, VI?UELA A, GRUNDBERG E, NELSON C P, MEDURI E, et al. Cigarette smoking reduces DNA methylation levels at multiple genomic loci but the effect is partially reversible upon cessation. Epigenetics, 2014, 9(10): 1382-1396.

[50] ZWAMBORN R A J, SLIEKER R C, MULDER P C, ZOETEMELK I, VERSCHUREN L, SUCHIMAN H E D, TOET K H, DROOG S, SLAGBOOM P E, KOOISTRA T, et al. Prolonged high-fat diet induces gradual and fat depot-specific DNA methylation changes in adult mice. Scientific Reports, 2017, 7: 43261.

[51] CHOI I, STEIBEL J P, BATES R O, RANEY N E, RUMPH J M, ERNST C W. Identification of carcass and meat quality QTL in an F(2) duroc × pietrain pig resource population using different least-squares analysis models. Frontiers in Genetics, 2011, 2: 18.

[52] KIM J J, ZHAO H H, THOMSEN H, ROTHSCHILD M F, DEKKERS J C M. Combined line-cross and half-sib QTL analysis of crosses between outbred lines. Genetical Research, 2005, 85(3): 235-248.

[53] SESSA W C. The nitric oxide synthase family of proteins. Journal of Vascular Research, 1994, 31(3): 131-143.

[54] HATTORI Y, HATTORI S, KASAI K. Lipopolysaccharide activates Akt in vascular smooth muscle cells resulting in induction of inducible nitric oxide synthase through nuclear factor-kappa B activation. European Journal of Pharmacology, 2003, 481(2/3): 153-158.

[55] SILVA G B, GARVIN J L. Akt1 mediates purinergic-dependent NOS3 activation in thick ascending limbs. American Journal of Physiology Renal Physiology, 2009, 297(3): F646-F652.

[56] KOEBKE K J, PAULY D J, LERNER L, LIU X E, PACHECO A A. Does the oxidation of nitric oxide by oxyMyoglobin share an intermediate with the metMyoglobin-catalyzed isomerization of peroxynitrite? Inorganic Chemistry, 2013, 52(13): 7623-7632.

[57] LOKE K E, LAYCOCK S K, MITAL S, WOLIN M S, BERNSTEIN R, OZ M, ADDONIZIO L, KALEY G, HINTZE T H. Nitric oxide modulates mitochondrial respiration in failing human heart. Circulation, 1999, 100(12): 1291-1297.

[58] RAYNER B S, HUA S S, SABARETNAM T, WITTING P K. Nitric oxide stimulates myoglobin gene and protein expression in vascular smooth muscle. The Biochemical Journal, 2009, 423(2): 169-177.

Integrated Aanalysis of Genome and DNA Methylation for Screening Key Genes Related to Pork Quality Traits

ZHAO ZhenJian, WANG Kai, CHEN Dong, SHEN Qi, YU Yang, CUI ShengDi, WANG JunGe, CHEN ZiYang, YU ShiXin, CHEN JiaMiao, WANG XiangFeng, TANG GuoQing

College of Animal Science and Technology, Sichuan Agricultural University/Key Laboratory of Livestock and Poultry Multi-Omics, Ministry of Agriculture and Rural Affairs/Farm Animal Genetic Resources Exploration and Innovation Key Laboratory of Sichuan Province/State Key Laboratory of Swine and Poultry Breeding Industry, Chengdu 611130

【Background】 The meat quality traits of pigs are important economic traits. Studying the molecular mechanisms that affect various meat quality traits and discovering key genes can guide genetic improvement of pigs and have significant implications for improving pork quality. Currently, those researches on meat-related mechanisms were mainly based on DNA genomics, while studies on DNA methylation related to meat quality traits, as well as integrated analysis of combining genomics and methylation, are scarce. 【Objective】 The potential key genes affecting pork meat quality were screened and identified by combined analysis of genome and DNA methylome, so as to provide a reference for the genetic improvement of pork quality. 【Method】 In this study, 28 meat quality traits of the ongissimus dorsi muscle of 140 Large White pigs were examined. Epigenome-wide association analysis (EWAS) and genome-wide association analysis (GWAS) were performed to identify CpG and SNP sites significantly associated with each trait. Subsequently, the conditional association analysis was performed by using SNPs as covariates on the significantly associated sites overlapping between GWAS and EWAS to further identify CpG sites with independent effects. Association analysis was then conducted to identify methylation quantitative trait loci (meQTL) by using the methylation levels of CpG sites as the dependent variable and SNPs as the independent variable. Finally, the cis-methylation quantitative trait loci (cis-meQTL) were used as instrumental variables for Mendelian randomization analysis to infer the causal relationship between cis-meQTL and phenotypes, while potential key genes at the loci were annotated and identified. 【Result】 (1) The significant associated sites were identified in the same genomic regions for the meat quality traits, namely yellowness value at slaughter 45 minutes postmortem (b45min), drip loss (DL), and docosahexaenoic acid (C22:6n-3), by both EWAS and GWAS. (2) After conditional association analysis, seven CpG sites for b45minand one CpG site for DL remained significant, while the three sites for C22:6n-3 were no longer significant after using SNPs as covariates, indicating that the significant associations of the seven CpG sites for b45minand one CpG site for DL identified by EWAS were not influenced by nearby significant SNPs. (3) A total of ten meQTL were identified for the seven CpG sites for b45minand one CpG site for DL, but the majority were trans-meQTL, with only one CpG site (SSC12:44 254 675 bp) identifying a cis-meQTL, suggesting that this site might be regulated by nearby SNPs. (4) Mendelian randomization analysis showed a causal relationship between the CpG site (SSC12:44 254 675 bp) and the b45minphenotype. (5) annotation of the locus revealed that the nearest gene to the CpG site (SSC12:44 254 675 bp) and its cis-meQTL was NOS2, and the CpG site was located within the NOS2 gene. 【Conclusion】Based on the integrated analysis of DNA methylation and genomics data, this study proposed that the NOS2 gene might be a key candidate gene for meat color traits. DNA methylation and SNP jointly regulated gene expression, thereby affecting the expression of genes related to meat color traits.

epigenome-wide association study (EWAS); genome-wide association study (GWAS); meat quality traits; DNA methylation; meQTL

10.3864/j.issn.0578-1752.2024.07.014

2023-10-18；

2023-12-31

四川省科技廳項目（2020YFN0024, 2021ZDZX0008, 2021YFYZ0030）、四川省豬創新團隊項目（sccxtd-2022-08）

趙真堅，E-mail：530773281@qq.com。通信作者唐國慶，E-mail：tyq003@163.com

（責任編輯 林鑒非）