大型迪庆藏猪不同生长阶段肌肉生长差异基因及其调控通路分析

聂靖茹 张 博 马 黎 张 浩 李国美 严达伟 董新星*

(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

的表达下调，肌肉丙酮酸脱氢酶复合物活性增强，对碳水化合物氧化和葡萄糖摄取增加，促进肌纤维变大。