不同性別大白豬肌肉全基因組高分辨率單堿基甲基化差異分析

郭添福，張志燕，陳 冬，姚天雄，肖石軍，黃路生

(江西農業大學省部共建豬遺傳改良與養殖技術國家重點實驗室，南昌 330045)

DNA甲基化作為一種關鍵的表觀遺傳調控機制，在調節基因表達、基因組印記、細胞分化和胚胎發育等生物過程中發揮重要作用[1]。所有組織或細胞普遍存在DNA甲基化和非甲基化機制，不同性別DNA甲基化水平是否存在差異是一個普遍關注的研究熱點。El-Maarri等[2]通過比較不同性別健康人體全血基因組ALU、LINE-1 重復序列的甲基化水平差異，證實其所有CpG島的甲基化水平男性均顯著高于女性。Fuke等[3]和Shimabukuro等[4]都發現，外周血白細胞男性甲基化水平顯著高于女性。此外，Xiao等[5]發現，人類肝中DNMT3a表達水平在兩性之間相似，而女性的DNMT3b表達水平顯著高于男性。以上結果顯示，性別可能是影響DNA甲基化狀態的重要因素之一，但有研究認為，不同性別的DNA甲基化水平無顯著差異[6-7]。白小青等[8]研究發現，30日 齡榮昌公豬的DNA甲基化含量高于母豬，但在統計上無顯著性差異。隨著新一代高通量測序技術的發展，全基因組甲基化測序在動物研究中越來越受到關注，但到目前為止，還未見大白豬不同性別間全基因組高分辨率單堿基甲基化水平差異分析的報道。

本研究利用全基因組重亞硫酸鹽測序法(whole genome bisulfite sequencing，WGBS)對不同性別大白豬肌肉組織的甲基化水平進行了分析，初步揭示豬性別差異的表觀調控機制，為深入開展不同性別豬肌肉的分子調控研究奠定理論基礎。

1 材料與方法

1.1 試驗材料

本試驗中大白豬來自江西樟樹綠環牧業集團法系育種核心群，屠宰210 日齡健康公母豬各 2頭，采集背最長肌組織，置于-80 ℃冰箱中備用。

1.2 大白豬背最長肌全基因組重亞硫酸鹽測序(WGBS)

1.2.1 DNA提取及樣品檢測 利用標準的苯酚-氯仿法提取基因組DNA，使用分光光度計測量DNA的濃度，瓊脂糖凝膠電泳分析DNA降解程度以及是否有RNA污染，Nanodrop檢測DNA的純度(OD260 nm/OD280 nm比值)，用Qubit對DNA濃度進行精確定量。

1.2.2 文庫構建 檢測合格的基因組DNA樣品，首先用Bioruptor (Diagenode, Belgium) 打斷成平均大小為 200～300 bp的片段，對打斷后的DNA片段進行末端修復、3′端加A堿基，并連接甲基化接頭；隨后進行Bisulfite處理(采用EZ DNA Methylation Gold Kit, Zymo Research)；用 2% 的瓊脂糖凝膠電泳選擇片段；用QIAquick Gel Extraction kit(Qiagen)回收DNA片段。

1.2.3 文庫質量檢測 文庫構建完成后，用Agilent 2100 檢測文庫DNA片段的完整性及插入片段大小，符合預期后，使用Q-PCR方法對文庫的有效濃度進行檢測，以保證文庫質量。

1.2.4 上機測序 DNA文庫檢測合格后，把不同文庫按照有效濃度及目標下機數據量的需求做pooling后進行Hiseq測序，測序策略為雙末端測序。測序地點為北京諾禾致源科技股份有限公司。

1.3 大白豬背最長肌全基因組測序數據的生物信息學分析

1.3.1 參考序列比對分析 將獲得的原始測序序列數據(Raw data)進行過濾，包括剔除帶接頭的reads，剔除未能確定堿基比例超過 10% 的reads，剔除低質量的reads(質量值沒有超過20的堿基數占整個read的50%以上的reads)。得到的過濾后數據通過比對軟件BSMAP[9](采用默認參數)比對到豬的參考基因組上(Sscrofa 11.1版本，下載地址：ftp://ftp.ncbi.nlm.nih.gov/genomes/Sus_scrofa/Assembled_chromosomes/seq)，并計算每個樣品的比對率和bisulfite(BS) 轉化率等信息。

目前，微波檢測領域的研究主要集中在微波介質材料的研究上.溶液物質在微波場中的行為與自身的極性有著密切的關系.極化程度可以由介電常數加以表示.溶液物質的介電常數通常與外界因素如濃度，溫度和頻段等存在一個確定的非線性關系[4].可以利用這種關系確定液體的成分.

1.3.2 測序深度統計 每個堿基位點比對到染色體上的讀長數即為該位點的測序深度。計算CG、CHH和CHG 3種序列環境(H代表A、T和C 堿基)下C位點的測序深度。

1.3.3 不同序列環境下胞嘧啶甲基化分析 對包括甲基化位點(mC)在內共 10 bp堿基統計作圖(上游4個堿基，下游5個堿基)，研究不同序列環境下甲基化胞嘧啶上下游的序列特征，可為甲基化預測提供依據。

1.3.4 甲基化組DMR分析 差異甲基化區域(differentially methylated region，DMR)是在多樣本間甲基化程度存在差異的一段基因組區域。DMR是一種重要的表觀遺傳學研究對象，可通過多種途徑調控基因的表達，從而影響生物學功能。本試驗采用SMART軟件(2.2.1版本，采用默認參數)進行DMR分析[10]。

1.3.5 聚類分析 采取層次聚類法分析樣本間DMRs的相互關系及差異趨勢，以探究甲基化差異的生物學意義。

1.3.6 DMRs相關基因的GO和KEGG富集分析 利用BEDTOOLS2軟件(2.25.0版本)進行基因注釋，得到差異甲基化基因(differentially methylated gene，DMG)。利用DAVID數據庫[11]將DMGs進行生物本體分析(GO)和代謝通路(KEGG)富集分析。對DMGs進行富集分析，可以為生物學過程相關研究提供候選基因。

2 結 果

2.1 豬背最長肌全基因組參考序列比對分析

2.1.1 豬背最長肌全基因組Reads與參考基因組比對情況統計 對 4個樣品進行WGBS測序，平均每個樣品產出 150.765 Gb過濾后數據。經計算，平均每個樣品得到 2 374 661 943 bp的唯一匹配堿基，平均唯一匹配率為 96.67%，平均BS轉化率為 99.62%，30×平均覆蓋率為 73.45%，結果見表1。

表1讀長與參考基因組比對情況

Table1Thecomparisonofreadstothereferencegenome

指標Index公豬_1 Boar_1公豬_2Boar_2母豬_1Sow_1母豬_2Sow_2平均Average總堿基數/bp Clean base number2 456 425 2042 456 371 7242 456 394 4932 456 380 0262 456 392 862匹配上堿基數/bp Mapped base number2 439 649 4742 440 580 7932 424 302 8792 439 262 3412 435 948 872匹配率/% Mapping rate99.3299.3698.6999.3099.17唯一匹配上堿基數/bp Unique mapped bases number2 379 110 1622 381 498 5402 359 996 4842 378 042 5862 374 661 943唯一匹配率/% Uniquely mapping rate96.8596.9596.0896.8196.67BS轉化率/% BS conversion rate99.5899.6399.6299.6799.62平均深度(×) Average depth26.5631.7822.9930.7128.0130× 覆蓋率/% Coverage rate73.6073.6773.1473.4073.45

2.1.2 豬背最長肌全基因組測序深度統計 根據單個堿基位點的讀長數計算測序深度及其累積分布。分析顯示，位于測序深度為(10×)～(30×)區段的位點最多，基本呈現正態分布趨勢，見圖1B；超過 80% 的全基因組位點擁有 15× 以上測序深度，見圖1A，說明單位點覆蓋率高，測序范圍廣，結果可信。

圖1 大白豬背最長肌基因組測序深度分布(A)和累積分布(B)
Fig.1 Genome coverage distribution(A) and cumulative distribution(B) of sequencing depth of genome in Large White pigs LD

柱狀圖表示各染色體讀長的平均測序深度，散點圖表示每條染色體讀長的覆蓋率
Histogram represent the average sequencing depth of reads in different chromosomes, scatter diagram represent coverage rate of reads on each chromosome
圖2 覆蓋率和讀長在各染色體上的分布
Fig.2 Distribution of reads and coverage rate on each chromosome

2.2 豬背最長肌全基因組甲基化水平分析

2.2.1 豬背最長肌全基因組甲基化C位點分布分析 在全基因組范圍約有 5% 的胞嘧啶被甲基化(mC/C)。根據胞嘧啶(C)序列特征可以將其分為CG、CHH和CHG(H代表A、T或C 堿基)3種序列環境[12]，CG環境甲基化占總甲基化的比率(mCG/mC)在3種環境中最高；甲基化主要發生在CG環境，mCHG/CHG和mCHH/CHH環境中甲基化所占比例均不超過 1%；母豬的平均mCG/CG水平要高于公豬。大白豬全基因組胞嘧啶甲基化情況見表2。

2.2.2 豬背最長肌全基因組甲基化密度圖譜 不同性別不同序列環境的全基因組甲基化密度見圖3，結果顯示，不同個體部分染色體以及不同序列環境的甲基化密度存在差異。

2.2.3 豬背最長肌全基因組染色體甲基化水平分析 所有個體不同染色體甲基化水平的差異和趨勢見圖4，除母豬_1的CHG與CHH環境的甲基化分布情況外，其余個體不同環境下甲基化分布趨勢基本一致。所有個體的甲基化程度(mC/C)均在12號染色體最高。

2.2.4 豬背最長肌全基因組基因功能區域甲基化水平分析 為更好地理解基因3′UTR、5′UTR、外顯子、內含子和啟動子(轉錄起始位點上游2kb的區域)等不同功能元件的甲基化水平，將每個基因的各個功能區域分別等分成20個單位長度，再對所有基因的功能元件區域相應單位長度的C位點甲基化水平求平均值后進行統計，見圖5，結果顯示，所有個體不同基因功能元件在3種環境的甲基化水平均相似；不同基因功能元件的甲基化水平存在差異，其中，CG環境下內含子和3′UTR的甲基化水平普遍較高，5′UTR的甲基化程度較低。

表2大白豬全基因組胞嘧啶甲基化情況

Table2Genome-widemethylationlevelofcytosineinLargeWhitepigs

組別Group公豬_1Boar_1公豬_2Boar_2平均Average母豬_1Sow_1母豬_2Sow_2平均AveragemC/C/%4.815.054.934.804.904.85mCG數/個mCG number 39 048 95839 519 68039 284 31939 307 46640 134 30839 720 887mCG/mC/%87.6588.488.0286.0888.6987.39mCG/CG /%70.8671.3371.09572.0872.9672.52mCHG數/個mCHG number1 299 8181 243 2571 271 5371 441 3371 215 1851 328 261mCHG/mC/%2.922.782.853.162.692.92mCHG/CHG /%0.630.600.6150.700.590.65mCHH數/個mCHH number4 204 2463 942 1044 073 1754 914 3063 903 4374 408 872mCHH/mC/%9.448.829.1310.768.639.69mCHH/CHH/%0.610.570.590.730.570.65

圈由外向內依次為參考基因組、公豬_1、公豬_2、母豬_1和母豬_2。紅色和綠色分別代表高甲基化水平和低甲基化水平
The circles from the outside to inside represent reference genome, boar_1, boar_2, sow_1 and sow_2, respectively. Circles display high methylation level in red and low methylation level in green
圖3 不同序列環境(mCG、mCHG和mCHH)甲基化密度圖譜
Fig.3 Methylation density maps of mCG, mCHG and mCHH

橫坐標代表染色體號，縱坐標代表不同序列環境下各染色體上甲基化的平均水平
The abscissa represent different chromosomes, and the ordinate represent methylation level of different chromosomes in different contexts
圖4 C位點平均甲基化程度在染色體上的分布
Fig.4 The average level distribution of methylation at C on chromosomes

2.3 豬背最長肌全基因組在不同序列環境下胞嘧啶甲基化分析

通過對不同序列環境下甲基化位點的序列偏好分析發現，mCHG中的第2位H大部分以T居多，A 次之；mCHH第2位H以C 為主，第3位H以T 為主，見圖6。

2.4 豬背最長肌全基因組DMR分析

2.4.1 豬背最長肌全基因組DMR鑒定結果及結構注釋 在不同性別大白豬的背最長肌組織中鑒定出CG、CHH和CHG環境下顯著的DMRs分別為118、0和1511個，共鑒定到1 629個DMRs。DMRs長度主要集中在20～300 bp，最長的出現在第 4號染色體上，達到 1 104 bp；越長的DMRs出現的頻率越低，見圖7A；性別差異分析顯示，母豬的DMRs甲基化水平明顯高于公豬，見圖7B；不同染色體DMRs長度分布呈多態性，見圖7C。

通過對DMRs進行注釋，共鑒定出 841 個DMGs，CG和CHH環境下發生甲基化最多的元件為基因間區，其次是內含子區，具體結果見圖8(圖中不同顏色代表不同功能元件所占的比率)。

2.4.2 豬背最長肌全基因組DMRs聚類分析 對不同樣本之間的DMRs進行聚類分析發現，相同性別的DMRs聚類在一起，不同性別間DMRs聚類差異顯著，見圖9。

2.4.3 豬背最長肌全基因組DMGs相關GO和KEGG富集分析 經過GO富集分析，共檢測到顯著相關的GO條目(P<0.05)171 個，涉及生物學過程、細胞組成以及分子功能3大類，這些DMGs顯著富集在軸突導向(GO: 0007411)、細胞連接(GO: 0030054)、細胞黏附(GO: 0050839)等相關過程。圖10 展示了 30 個最顯著GO條目。對不同性別DMGs進行KEGG富集分析，共得到 10 個相關信號通路(P<0.05)，顯著富集在ECM受體互作(hsa04512)等通路，詳情見圖11。

將每個基因的各個功能區域分別等分成20個單位長度，再計算每個單位長度的甲基化水平求平均值繪制成折線圖。橫坐標代表不同的基因組功能元件，縱坐標為不同序列環境下各功能元件的甲基化水平。不同顏色代表不同序列環境(CG、CHG和CHH)
Each element of each gene was divided into 20 bins and the average of all bins in each gene elements was obtained. The abscissa represent different gene elements, and the ordinate represent methylation levels of elements in different sequence contexts. Different colors represent different sequence contexts (CG, CHG and CHH)
圖5 甲基化水平在不同基因組功能元件上的分布
Fig.5 Distribution of methylation level in different genomic components

橫坐標為堿基位置，mC的位置定義為0，縱坐標為相應位置堿基的富集程度，每個位置的總高度為堿基的序列含量
The abscissa represent bases position, the position of mC is defined as 0, and the ordinate represent the degree of base enrichment, the total height of each position is the sequence content of the base
圖6 不同序列環境下的胞嘧啶甲基化
Fig.6 Cytosine methylation in different sequence contexts

A圖橫坐標代表DMRs長度，縱坐標代表對應DMRs長度的頻率；B圖橫坐標為不同性別，縱坐標代表DMRs甲基化水平分布；C圖橫坐標代表染色體編號，縱坐標代表DMRs長度分布
The abscissa represents DMRs length, the ordinate represents the frequency of different DMRs length in figureA；The abscissa represents gender, the ordinate represents the methylation level of DMRs in figureB；The abscissa represents different chromosome, and the ordinate represents the distribution of DMRs length in figure C
圖7 DMRs長度和水平分布
Fig.7 Distribution of DMRs length and level

圖8 不同序列環境DMRs所在基因功能元件的百分比
Fig.8 The percentage of DMRs in different gene elements in different sequence contexts

為了更好地理解DNA甲基化與不同性別肌肉差異的相關性，本研究設置2個限制因素進行候選基因篩選：1)選擇GO分析中富集于肌肉相關通路的DMGs；2)第一步選擇的DMGs需要顯著富集在KEGG分析中的信號通路(剔除癌癥相關通路)。結果得出，有 2個GO條目與肌肉相關通路有關，即平滑肌細胞遷移負調控(GO:0014912)和肌肉器官發育(GO:0007517)，包含 14 個DMGs，其中有5個DMGs富集在KEGG分析中的信號通路(DMD、ITGA7、LMNA、SEMA6D和SLIT2)，具體信息見表3。

行代表個體號,列代表個體的DMR
Each row represent an individual and each column represent an individual DMR
圖9 大白豬不同樣本DMRs聚類分析
Fig.9 DMRs cluster analysis in different samples of Large White pigs

縱坐標代表GO條目；橫坐標代表-lg(P-value)
The ordinate represents the GO terms, the abscissa represents the -lg(P-value)
圖10 最顯著30個GO條目的柱狀圖
Fig.10 The histogram of top 30 GO terms

縱坐標代表KEGG條目；橫坐標代表富集因子。圈的大小代表通路的基因數，圈越大，基因數越多。不同顏色代表-lg(P-value),顏色越紅，-lg(P-value)越大
The ordinate represents the KEGG terms and the abscissa represents the rich factor. The size of circles represent the number of genes contained in the particular class, the larger the circle indicate the more genes. The different colors represent the -lg(P-value), the redder the circles indicate the higher -lg(P-value)
圖11 KEGG條目氣泡圖
Fig.11 The bubble chart of KEGG terms

3 討 論

豬是人類最主要的動物蛋白來源，肉質作為養豬生產中的一個重要經濟性狀，已經越來越受到研究者、生產者和廣大消費者的重視。

DNA甲基化作為重要的表觀遺傳學現象之一，能引起染色質結構、DNA構象和轉錄效率等多方面的改變，從而調控基因表達[13]。近年來，豬的全基因組甲基化方面雖有研究，但主要采用的是DNA甲基化免疫共沉淀測序法[14-23]和簡化代表性重亞硫酸鹽測序法[23-25]，雖然有極少數研究利用“金標準”重亞硫酸鹽進行甲基化測序法[26-28]，但覆蓋度太低；主要集中在不同品種、不同組織或不同處理方式之間的甲基化差異分析，尚未見不同性別間肌肉全基因組甲基化差異的報道。本試驗利用WGBS研究不同性別大白豬的背最長肌全基因組高分辨率單堿基甲基化差異，系統地比較了不同性別大白豬背最長肌全基因組甲基化的差異。

本研究發現，基因組中大約有 5% 的C位點被甲基化，不同性別間胞嘧啶甲基化水平差異不顯著，與此前相關報道類似[8]。CG環境甲基化比率在3種環境中最高，與此前報道的結果一致[20-22,26-27]。不同性別功能元件區的全基因組甲基化模式相似，CG、CHG和CHH 3種序列環境的甲基化差異可能與不同遺傳元件差異有關。

大部分DMRs為20～300 bp，說明小片段區域C位點甲基化對基因表達的調節起重要作用。不同序列環境下DMGs的功能元件區分布一致，CG和CHH環境主要集中在基因間區和內含子區，其他功能元件區分布均較少。出現這種現象的原因可能是，基因間區和基因的內含子區約占整個基因組的99%，所以大多數的甲基化位點出現在基因間區和內含子區域；另一方面，啟動子等功能區域的甲基化會阻抑基因的表達，尤其是眾多的看家基因以及早期就需要表達的基因，這些基因的功能區域都基本沒有發生甲基化，而內含子和基因間區的甲基化所涉及的調控功能較少，故非功能區間的甲基化位點更多。

在GO分析中，841 個DMGs中有 757 個富集在生物學過程、細胞組成和分子功能。通過2個嚴格的因素(一是選擇GO分析中富集于肌肉相關通路的DMGs；二是第一步選擇的DMGs需要顯著富集在KEGG分析中的信號通路(剔除癌癥相關通路))篩選出5個可能與不同性別肌肉調節相關的DMGs：DMD、ITGA7、LMNA、SEMA6D 和SLIT2。DMD基因(抗肌肉萎縮蛋白基因)發生無義突變或框移突變而引起進行性肌營養不良癥，該癥是一種常見的X連鎖隱性遺傳的神經肌肉疾病，包括Duchenne型肌營養不良和Becker型肌營養不良，主要表現為容易跌倒，蹲站困難，喪失獨立行走能力等[29-33]。ITGA7基因(整合素α7基因)是整合素家族的重要成員之一，編碼的整合蛋白α7為心肌和骨骼肌層黏蛋白受體，與肌萎縮蛋白糖蛋白復合物類似，連接細胞外基質到肌動蛋白細胞骨架內部，該蛋白具有治療DMD的潛質[34]。LMNA基因(核纖層蛋白基因)編碼的核纖層蛋白對核穩定、染色質結構和基因表達有重要作用，為骨骼肌衛星細胞增殖所必需[35-38]，協助維持骨骼肌的體積和強度[39]，該基因突變是引起阻滯型家族性擴張型心肌病最常見的致病基因[40]。SLIT2 基因為Slit家族成員之一，編碼的Slit2蛋白能阻止血小板源性生長因子誘導主動脈平滑肌細胞遷移[41]。SEMA6D 基因是Semaphorin家族成員之一，在肌肉、胚胎和腦等組織中均有表達，該基因與神經系統的發生及軸突的形成有密切關系[42]；在胚胎的發育過程中，該基因通過與不同受體結合，能抑制或誘導心不同部位上皮細胞的遷移，從而在胚胎心血管的形態發生中起重要作用[43]，未見該基因參與肌肉發育相關過程的報道，其生物學功能有待進一步探討。

4 結 論

本研究全面構建了大白豬背最長肌高分辨率單堿基的全基因組甲基化圖譜，并解析了性別差異的甲基化區域和基因，共檢測到 1 629 個DMRs和 841 個DMGs，雖然不同性別總體甲基化水平相似，但母豬的DMRs平均甲基化水平明顯高于公豬。通過富集分析，獲得 171 個GO條目和 10 個KEGG條目，篩選出 5個與不同性別肌肉相關的候選基因，可為未來從遺傳學和表觀遺傳學角度研究豬肉質性狀及人類試驗動物模型的開發提供重要參考。