大型迪慶藏豬不同生長階段肌肉生長差異基因及其調控通路分析

聶靖茹 張 博 馬 黎 張 浩 李國美 嚴達偉 董新星*

(1.云南農業大學 動物科學技術學院，昆明 650201； 2.中國農業大學 動物科學技術學院，北京 100193； 3.云南農業職業技術學院 畜牧獸醫學院，昆明 650212)

豬肉是我國居民肉類消費的主體，隨著我國社會經濟的發展和消費結構升級，對豬肉的追求已從“有肉吃”向“吃好肉”轉變。我國地方豬種有肉質好的優點也有瘦肉組織生長慢的不足，生長慢、瘦肉率低是制約其大規模應用的瓶頸，保持肉質優良的同時提高產肉量是我國地方豬種產業化急需解決的關鍵問題，篩選地方豬種肌肉生長的關鍵基因并解析其調控機制，對地方豬種的選育利用具有重要意義。

肌肉生長是一復雜的生理過程，一系列起正調控或負調控的基因或因子以多種復雜方式(時空表達、信號級聯、轉錄調控和反饋機制等)精確調控對肌肉生長至關重要。隨著高通量測序技術的發展，通過不同豬種的比較發現了一批影響豬肌肉生長的轉錄因子和調節肌肉生成的基因調控網絡。Shang等對藏豬、烏金豬和大白豬進行轉錄組和蛋白組測序發現，與成肌細胞分化和肌纖維形成有關的

CRYAB

、

FSCN1

和

MAPK12

等20個基因可能在決定豬的出生后生長速度和理論體重中起重要作用。Xu等比較體重達90 kg民豬和長白山野豬，發現4 522個差異基因，富集到豬骨骼肌的肌纖維發育、分化和生長。Zhang等對長白豬和五指山小型豬出生后18、21和28 d的背最長肌進行了亞硫酸氫鹽全基因組測序，發現

Tet1

參與了肌生成素啟動子的去甲基化，促進小鼠成肌細胞系C2C12的永生化和豬胚胎成肌細胞分化。Zhao等比較了通城豬和約克夏豬胚胎30 d到出生后5周的背最長肌，發現

CXCL10

、

EIF2B5

、

PSMA6

、

FBXO32

和

LOC100622249

基因在通城豬的肌肉調節網絡中起著重要作用，而

SGCD

、

ENG

、

THBD

、

AQP4

和

BTG2

基因在大白豬中起主要作用，表現出品種特異性和發育依賴的差異表達模式。但上述研究均是進行品種間的比較，未見利用品種內的差異群體揭示我國地方豬肌肉生長發育調控機制的報道，能成熟應用于提高我國地方豬產肉量的功能基因或分子標記是非常有限的。

迪慶藏豬是我國特有的高原型豬種之一，具有抗逆性強、肉質好等優點但也有生長慢、產仔少的不足。迪慶藏豬可分為大、中、小三種類型，大型豬又稱大架子豬，成年體質量70～150 kg，育肥期日增重200～250 g；中型豬又稱二虎頭豬，成年體質量55～65 kg，育肥期日增重100～200 g；小型豬又稱荷包豬，成年體質量45～55 kg，育肥期日增重100～120 g，目前關于迪慶藏豬的報道主要集中在低氧適應、生長和肉質，未見其肌肉生長發育功能基因的報道。本研究采用RNA-seq技術篩選大型迪慶藏豬不同生長階段(S1：10～40 kg、S2：40～80 kg、S3：80～120 kg)背最長肌肌肉生長相關的差異表達基因并分析其調控網絡，豐富我國地方豬種肌肉生長發育調控基礎數據，為迪慶藏豬的遺傳改良提供參考。

1 材料與方法

1.1 試驗材料

選擇胎次相同，出生日期及體重相近的純種大型迪慶藏豬36頭，隨機分為3組，每組12頭(豬只分組情況見表1)，在云南省香格里拉市綠源生態種養專業合作社藏豬養殖基地進行育肥試驗，分別在體重達40、80和120 kg左右時屠宰，每組選擇接近平均體重的3頭，采集背最長肌組織，液氮速凍運回實驗室，保存于-80 ℃冰箱，用于提取總RNA。

表1 試驗豬只分組
Table 1 Group of experimental pigs

組別Group開始體重/kgInitial weight結束體重/kgFinal weight10～40 kg日增重/gAverage dailygain of S140～80 kg日增重/gAverage dailygain of S280～120 kg日增重/gAverage dailygain of S3試驗1組Group 18.90±1.2537.21±3.45225±18 試驗2組Group 29.00±1.3580.50±6.27238±18 B313±25 A試驗3組Group 39.05±1.02111.50±10.99 235±11 bC 342±45 aA297±56 bA

注：同行數據為相同字母或不標字母表示差異不顯著(>0.05)，標有不同小寫字母表示差異顯著<0.05，大寫字母表示差異極顯著<0.01。

Note: The same letter or no letter indicates that the difference is not significant (>0.05). Different lowercase letters represent significant differences between treatments (<0.05). Different uppercase letters indicate extreme significant differences (<0.01).

1.2 試驗方法

1

.

2

.

1

RNA提取和質量檢測

采用TRIzol結合RNApure Tissue Kit的方法，參照文獻[12]的步驟提取背最長肌總RNA，去除rRNA。檢測RNA質量濃度和純度，檢測合格的總RNA用于轉錄組測序。

1

.

2

.

2

去核糖體鏈特異性文庫構建及測序

從9 個檢測合格的RNA中各取3 μg建庫，將mRNA進行Oligo(dT)磁珠富集，反轉錄合成雙鏈cDNA文庫。對構建好的文庫進行質檢，質檢合格的文庫利用Illumina Hiseq2 000進行上機測序。

1

.

2

.

3

測序數據處理和分析對測序獲得原始數據(Raw reads)進行質量控制，去除污染和低質量片段后獲得過濾后的reads(Clean reads)，采用HISAT2軟件將Clean reads與豬的參考基因組(

Sus

scrofa

11.1)進行比對，對比對到基因組各染色體上的Reads進行統計，以|logFC|>1且

P

value<0.05作為篩選條件，將片段FPKM值歸一化進行差異基因分析。

1

.

2

.

4

差異基因的功能富集分析對篩選出的差異表達基因進行GO和KEGG功能富集分析，篩選與迪慶藏豬肌肉生長相關的差異基因，按富集度Count≥2且矯正

P

<0.05作為顯著富集基因的閾值。

1

.

2

.

5

差異基因的STEM分析由于本試驗中基因表達譜數據只有3 個時間節點，故采用短時間序列表達分析(Short time-series expression miner, STEM)對差異基因進行趨勢分析，數據過濾標準為3個時間節點兩兩比較的比較組內表達差異倍數在2 倍以上且

P

<0.05，對顯著差異表達的模塊進行GO和KEGG功能富集分析。

1

.

2

.

6

差異基因互作網絡分析

利用Cytoscape 3.7.1軟件挖掘差異基因間的互作關系，構建差異基因的互作網絡。

1

.

2

.

7

qPCR定量驗證用Primer Premier 5.0軟件設計qPCR引物，利用反轉錄試劑盒將RNA反轉錄為cDNA，以

GAPDH

作為內參，對篩選到的重要差異基因進行實時熒光定量檢測，篩選差異基因的引物信息見表2。

表2 差異表達基因引物信息
Table 2 Primer information of differential expression genes

序號Serial number名稱Name引物序列(5'-3')Primer sequence (5'-3')退火溫度/℃Tm產物長度/bpLength1ATF3F: TGCTAACCTGACACCCTTTGTR: GCACTCCGTCTTCTCCTTCTT602382CCL2F: TTGAATCCTCATCCTCCAGCAR: GATCTCCTTGCCCGCGATG592423CCN1F: CACCAATGACAACCCCGACTR: ACTTGGGCCGGTACTTCTTC601854EGR1F: AGTTTGCCAGGAGCGATGAAR: AGGCCACACTTTTGTCTGCT60 845FOXO1F: AAGAGCGTGCCCTACTTCAAR: CCTGGGGGATTTCCCACTCT601806NOR-1F: AATGCCCTTGTCCGAGCTTTR: TGGAGGCTGTGAGAAGGTTG601367PDK4F: GCTGGTGACTGGTGTATCCCR: TCACACGCACACATTCAGGA601418PPARDF: AGAACCGCAACAAGTGCCAGR: AGCTTCCTTTTCTCTGCCTCG611099PPARGC-1F: AACCCACAGAGACCCGAAACR: CCCTTGGGGTCATTTGGTGA6014510SFRP1F: GAACAAGGACACGTTGCCAGR: AATGGAGAGTGCCACAAGCA6011911SOX9F: GAAAGTCGGTGAAGAACGGCR: CTTGAAGATGGCGTTGGGAG598012STAT1F: TCCGACACCTGCAACTGAAAR: GGCAGAGAGGTGGTCTCAAG6015113GAPDHF: TCGGAGTGAACGGATTTGR: CCTGGAAGATGGTGATGG60219

2 結果與分析

2.1 迪慶藏豬背最長肌RNA提取質量及測序數據統計

用1%瓊脂糖凝膠電泳檢測提取的總RNA質量，由圖1可見，樣品總RNA的28 S和18 S條帶清晰可見，質量達到RNA-seq上機測序的要求。RNA-seq測序數據顯示(表3)，每個樣品的Raw reads數都在48 M以上，Clean reads比率都超過了88%，Clean Q30 bases比率均在94%以上，3組樣品的Mapping率均在95%以上，符合生物信息學分析要求。

圖1 迪慶藏豬背最長肌總RNA瓊脂糖凝膠電泳圖Fig.1 RNA agarose gel electrophoresis in longissimus dorsi of large Diqing Tibetan pigs

2.2 迪慶藏豬背最長肌基因表達水平分析

利用篩選到的全部差異基因的FPKM值繪制小提琴圖(圖2)并對3 個試驗組的基因表達水平進行比較，組內基因表達水平基本相似，說明樣品表達量穩定，組間存在差異。

2.3 迪慶藏豬背最長肌mRNA差異表達分析

10～40 kg(S1)、40～80 kg(S2)和80～120 kg(S3)體重階段共得到19 490個差異表達基因，741個差異基因在S1階段特異表達，443個差異基因在S2階段特異表達，951個差異基因在S3階段特異表達(圖3)。

表3 RNA-seq測序數據統計
Table 3 RNA-seq data statistics

樣品Sample原始reads數Raw readsnumber有效reads數Clean readsnumber有效reads比例/%Clean readsrate有效Q30比例/%Clean Q30bases rate比對率/%Mapping ratio40 kg-1背最長肌S1-1LD49 985 88646 915 29693.8694.5996.1740 kg-2背最長肌S1-2LD50 819 64445 904 62890.3394.3795.9440 kg-3背最長肌S1-3LD48 429 06444 220 50091.3194.5195.7680 kg-1背最長肌S2-1LD49 312 96846 160 38693.6194.6996.1780 kg-2背最長肌S2-2LD53 814 98847 796 56088.8294.8896.1780 kg-3背最長肌S2-3LD51 209 20047 989 30493.7194.7996.36120 kg-1背最長肌S3-1LD51 524 07647 920 55693.0194.5596.06120 kg-2背最長肌S3-2LD51 219 00646 874 73691.5294.5095.50120 kg-3背最長肌S3-3LD49 398 23645 979 20693.0894.4696.06

LD表示背最長肌，黃色表示S1階段背最長肌表達量，綠色表示S2階段背最長肌表達量，粉色表示S3階段背最長肌表達量。 LD represents the longissimus dorsi muscle. Yellow represents the expression of LD in S1 stage, green represents the expression of LD in S2 stage, and pink represents the expression of LD in S3 stage.圖2 迪慶藏豬背最長肌差異基因表達水平分析圖Fig.2 Differential gene expression level analysis in longissimus dorsi of large Diqing Tibetan pigs

2.4 不同體重階段迪慶藏豬背最長肌顯著差異基因篩選

以|logFC|>1且

P

<0.05為篩選顯著差異基因的閾值，在S1與S2階段共檢測到730個基因顯著差異表達，其中，341個基因上調表達，389個基因下調表達(圖4(a))；在S2與S3階段中共檢測到981個基因顯著差異表達，其中，530個基因上調表達，451個基因下調表達(圖4(b))。

圖3 迪慶藏豬背最長肌三階段樣品間基因表達韋恩圖Fig.3 Venn diagram of gene expression between three-stage samples in longissimus dorsi of large Diqing Tibetan pigs

2.5 不同體重階段迪慶藏豬背最長肌顯著差異基因功能富集分析

S1vsS2階段：GO分析結果顯示，差異基因主要富集在骨骼肌細胞分化、典型Wnt信號通路、破骨細胞分化等(圖5(a))。KEGG分析顯示，差異基因主要富集在FOXO信號通路、TNF信號通路等(圖5(b))，其中

SFRP1

、

BCL3

、

ACTA2

和

CCN1

等基因在以上通路中富集，均下調表達，

ACVR2B

和

ATF3

等基因上調表達，這些基因與肌肉生長性狀相關。

(a)S1 vs S2階段差異基因火山圖。(b)S2 vs S3階段差異基因火山圖。

S2vsS3階段：GO分析結果顯示，差異基因主要富集在典型Wnt信號通路、成纖維細胞遷移的正調控、JUN激酶活性的負調控等(圖5(c))。KEGG分析結果顯示，差異基因主要富集在ECM受體相互作用、Wnt信號通路、PI3K-Akt信號通路等(圖5(d))，其中

PDK4

、

MYF6

和

ATF3

等基因在以上通路中富集，下調表達；

NFAT5

、

SFRP1

和

SCD

等基因在S2vsS3階段均上調表達，參與肌肉生長過程。

(a)S1 vs S2階段差異表達基因GO富集分析氣泡圖;(b)S1 vs S2階段差異表達基因KEGG富集分析氣泡圖;(c)S2 vs S3階段差異表達基因GO富集分析氣泡圖;(d)S2 vs S3階段差異表達基因KEGG富集分析氣泡圖。

2.6 迪慶藏豬背最長肌差異基因短時間序列表達分析(STEM)

對篩選到的顯著差異基因進行趨勢分析，由圖6 可見，共有4 個差異顯著的模塊，其中模塊11和14的差異基因總體聚為一類，為先上調表達后呈現輕微下調表達或不變；模塊9和10的差異基因聚為一類，為先上調表達后下調表達；其余模塊差異不顯著。

圖6 迪慶藏豬背最長肌STEM分析Fig.6 STEM analysis in longissimus dorsi of large Diqing Tibetan pigs

將模塊11和14的基因集合并，進行GO和KEGG富集分析，GO結果顯示(圖7(a))，差異基因主要在調控成纖維細胞生長因子受體信號通路、MAPK活性激活等過程中富集；KEGG結果顯示(圖7(b))，差異基因主要富集到PI3K-Akt信號通路、ECM受體相互作用等過程。將模塊9和10的基因集合并，對差異基因進行GO和KEGG富集分析，GO結果顯示(圖7(c))，差異基因主要富集到Wnt信號通路的調控、調節肌動蛋白細胞骨架組織；KEGG結果顯示(圖7(d))，差異基因主要富集到T細胞受體信號通路等過程。

2.7 不同體重階段迪慶藏豬背最長肌顯著差異基因互作網絡

2

.

7

.

1

S1vsS2差異基因互作網絡選擇功能富集分析中顯著富集通路中的基因，繪制基因網絡圖(圖8)，可以看出：

FOS

基因位于調控網絡的中心，與

NOR-1

、

EGR2

、

EGR1

、

NFAT5

、

CCL2

、

ARNTL

和

GHSR

等基因有互作關系。

NOR-1

基因上調表達，

FOS

、

EGR2

、

EGR1

和

CCL2

基因下調表達。

2

.

7

.

2

S2vsS3差異基因互作網絡S2vsS3階段迪慶藏豬背最長肌差異基因互作網絡見圖9，由圖可見：

MYC

、

SOX9

、

CCN2

、

FOXO1

、

PDK4

、

PPARGC-1

、

PPARD

、

ATF3

和

FZD7

基因位于網絡中心位置，核心位置的基因除

CCN2

上調表達外，其余基因均下調表達。

2.8 迪慶藏豬背最長肌肌肉生長相關差異基因定量驗證

隨機選擇12個差異基因進行qPCR驗證，

EGR1

、

SFRP1

、

FOXO1

、

CCL2

和

NOR-1

等12個基因在背最長肌組織中的表達量與轉錄組測序結果一致(圖10)。

3 討 論

3.1 S1vsS2階段肌肉生長差異基因

S1vsS2階段，

NOR-1

基因上調，

FOS

、

EGR1

、

EGR2

和

CCL2

基因下調。

NOR-1

基因(即

NR4A3

)是

NR4A

亞家族的一員，在骨骼肌和其他腎上腺素靶組織(包括脂肪和心臟)中表達，在C2C12骨骼肌細胞系中，

NOR-1

降低會導致肌生長抑制素的表達增加，肌生長抑制素是肌肉質量的負調節因子。Mey等研究表明，

NOR-1

參與介導肌肉的合成代謝，通過調節骨骼肌β-腎上腺素受體激動劑的表達來調節肌生長抑制素，與β-腎上腺素能誘導

NOR-1

表達及骨骼肌肥大的結果一致。肌衛星細胞激活時

NOR-1

基因表達上調，表明該基因與衛星細胞和肌肉生長有關，

NOR-1

在肌肉細胞分化中發揮作用，可作用于肌肉損傷修復。本研究中，

NOR-1

在S1vsS2階段表達上調，可能使肌生長抑制素基因表達下降，導致肌肉肥大和肌生成力增加，生長速度加快，與S1階段到S2階段日增重極顯著增加的結果一致。

FOS

、

EGR1

和

EGR2

基因均屬于立早基因(IEGs)，參與細胞生長和分化。

FOS

基因參與肌生成，據Zhao等研究，分化C2C12成肌細胞時，

FOS

基因下調會完全抑制肌管形成，抑制肌源性分化。

FOS

基因能夠調節

TNF

-

α

的表達，而

TNF

-

α

又通過

NF

-

κB

的激活和

IGF-1

信號通路的損傷來抑制肌源性分化，在分化成肌細胞時，

TNF

-

α

誘導的

NF

-

κB

激活導致分化標志肌原蛋白和肌球蛋白重鏈的表達降低，在分化肌管時，總蛋白和快速型肌球蛋白重鏈水平降低。

FOS

基因激活

IGF-1

信號通路的下游酶

Akt

，調控肌生成，通過抑制酪氨酸磷酸化下調肌生成素表達。

EGR1

與

EGR2

基因均屬于早期生長因子家族，在脊椎動物肌腱形成中起作用，Brent等研究表明，

EGRs

和肌腱膠原的表達受

FGF

的正向調節，雞肌腱細胞的分化又依賴于

FGF

肌肉源性信號。Canty等研究發現，

EGR1

在

EGR2

-/-突變小鼠中的表達降低，與過表達

EGR

組肌腱相比，缺乏

EGR

會導致膠原纖維束減少。

CCL2

可通過肌內巨噬細胞上調

IGF-1

的表達促進急性骨骼肌損傷修復，

CCL2-/-

小鼠對缺血或肌毒因子誘導的急性損傷的反應顯示巨噬細胞浸潤明顯減少，炎癥反應減少伴隨著肌肉再生不良，

CCL2

缺乏可導致急性損傷后肌肉炎癥明顯減輕，肌肉再生受損。本研究在S1 vs S2階段

FOS

、

EGR1

、

EGR2

和

CCL2

基因下調表達，可能抑制成肌細胞肌源性分化、肌生成力減弱、肌腱膠原纖維減少、肌肉再生能力降低，肌肉生長潛力下降。

(a)模塊11及14差異表達基因GO富集分析氣泡圖；(b)模塊11及14差異表達基因KEGG富集分析氣泡圖；(c)模塊9及10差異表達基因GO富集分析氣泡圖；(d)模塊9及10差異表達基因KEGG富集分析氣泡圖。

黃色表示位于網絡中心，橙色表示位于網絡外周。下同。 Yellow means it is located in the center of the network, orange means it is located on the periphery of the network. The same as below.圖8 迪慶藏豬背最長肌S1 vs S2階段差異基因互作網絡圖Fig.8 Interaction network diagram of differential expression genes at S1 vs S2 stage in longissimus dorsi of large Diqing Tibetan pigs

圖9 迪慶藏豬背最長肌S2 vs S3階段 差異基因互作網絡圖Fig.9 Interaction network diagram of differential expression genes at S2 vs S3 stage in longissimus dorsi of large Diqing Tibetan pigs

圖10 迪慶藏豬背最長肌差異表達基因qPCR定量驗證結果Fig.10 QPCR results of differential expression genes in longissimus dorsi of large Diqing Tibetan pigs

3.2 S2vsS3階段肌肉生長差異基因

S2vsS3階段，

SFRP1

和

SFRP2

基因表達上調，

FZD7

、

FOXO1

、

PDK4

和

PPARD

基因表達下調。

SFRP1

和

SFRP2

基因屬于

SFRP

家族，富含半胱氨酸結構域，該結構域與細胞表面卷曲受體的Wnt結合位點同源，可通過競爭性結合調控Wnt信號途徑而調節細胞生長和分化。經典Wnt信號通路在胚胎肌肉發育、成年骨骼肌再生及成肌細胞系體外增殖和分化過程中具有重要作用，

SFRP1

可作為Wnt信號的可溶性調節劑，參與控制骨骼肌發育過程中的細胞增殖和分化，Descamps等研究表明，

SFRP1

和

SFRP2

添加到C2C12或原代衛星細胞中可以抑制肌管形成，其作用是阻止成肌細胞進入終末分化過程。本研究中，

SFRP1

和

SFRP2

在S2 vs S3階段表達上調，可能阻止成肌細胞分化形成肌纖維，導致肌纖維數量減少，肌肉生長潛力下降。

FZD7

基因能激活PI3K-Akt通路，而IGF1-PI3K-Akt/PKB-mTOR通路作為肌肉生長的正調節器，既可刺激衛星干細胞擴增，又刺激肌纖維肥大，本研究從S2到S3階段

FZD7

表達下調，可能減緩肌纖維肥大、肌肉生長速度下降，與S2階段到S3階段日增重顯著降低的結果一致。

FOXO1

基因屬于forkhead轉錄因子家族，主要在骨骼肌、肝臟和脂肪組織表達，Allen等研究表明，

FOXO

轉錄因子可通過增加肌萎縮基因的表達進而刺激蛋白質水解，肌生長抑制素基因是

FOXO1

作用的下游靶點，

FOXO1

可以結合到肌生長抑制素基因上游啟動子區域的

FOXO

位點，增加肌生長抑制素mRNA表達和啟動子活性，抑制合成代謝過程，加速肌肉萎縮。此外，

FOXO1

還能夠上調肌肉特異性泛素連接酶的表達并減小肌肉的大小，本研究中，

FOXO1

基因在S2vsS3階段表達下調，可能減少肌生長抑制素mRNA表達，減緩肌肉萎縮進程。

PDK4

基因在骨骼肌中特異性表達，同時，

FOXO1

可與

PDK4

基因的啟動子區結合直接增加肌肉中

PDK4

mRNA的表達，另外，

PPARD

與

FOXO1

轉錄因子表達相關，可增加肌肉中

PDK4

蛋白含量，從而抑制肌肉丙酮酸脫氫酶復合物活性，進而減弱碳水化合物氧化和葡萄糖攝取，誘發胰島素抵抗。本研究在S2到S3階段

FOXO1

、

PDK4

和

PPARD

均下調表達，可能原因是

PPARD

下調導致

FOXO1

基因表達下調，而

FOXO1

下調導致與

PDK4

基因的啟動子區結合減少，肌肉中

PDK4

的表達下調，肌肉丙酮酸脫氫酶復合物活性增強，對碳水化合物氧化和葡萄糖攝取增加，促進肌纖維肥大。

4 結 論

大型迪慶藏豬從40到80 kg，

NOR-1

基因表達上調抑制了肌生長抑制素的表達，導致肌肉肥大和肌生成力增加，生長速度加快；

FOS

、

EGR1

、

EGR2

和

CCL2

基因下調，抑制成肌細胞肌源性分化、肌腱膠原纖維減少，肌肉生長潛力下降。從80到120 kg，

SFRP1

和

SFRP2

基因上調阻止成肌細胞終末分化，肌纖維數量減少，肌肉生長潛力降低；

FZD7

基因表達下調，減緩肌纖維肥大；

PPARD

下調導致

FOXO1

基因下調，一方面減少肌生長抑制素mRNA表達從而減緩肌肉萎縮進程，另一方面，

FOXO1

下調導致與

PDK4

基因的啟動子區結合減少，肌肉中

PDK4

的表達下調，肌肉丙酮酸脫氫酶復合物活性增強，對碳水化合物氧化和葡萄糖攝取增加，促進肌纖維變大。