转录组测序揭示不同生长速度番鸭胸肌组织的基因表达差异

张思雨，朱炜健，季从亮，郭利金，郑 茗，徐海平，聂庆华

（1. 华南农业大学动物科学学院/农业农村部鸡遗传育种与繁殖重点实验室，广东 广州 510642；2. 温氏食品集团股份有限公司，广东 云浮 527400）

【研究意义】畜禽的生长速度对于满足全球日益增长的蛋白质需求至关重要。我国生产的禽肉中，鸭肉产量仅次于鸡肉，提高鸭的生长性能可以带来极高的经济价值。通过基因选择育种来提高饲料转化效率和肌肉发育质量是一种卓有成效的方法，因此，研究调控肉鸭生长发育速度的相关基因具有很高的研究价值和应用价值。番鸭（Muscovy duck）相较于一般的家鸭具有体型大、产肉量大、饲料报酬高、肉质性能好和瘦肉率高等特点，是广泛养殖的肉鸭品种。然而，同一品种番鸭的生长速度可以表现出很大的差异，其中涉及到的遗传机制目前尚不明确，有待进一步研究。

【前人研究进展】畜禽肉的产量直接关系到养殖行业的经济效益，因此肌肉发育调控的研究一直是畜牧业研究的重点。家禽不同品种所表现的生长性能差异为肌肉发育调控相关研究提供了理想的模型。运用转录组测序技术（RNA-Seq）对生长性能表现出差异的品种进行测序分析，有利于挖掘与生长性状相关的调控基因，分析其涉及的生物功能和调控通路[1]。迄今为止，在鸭中已经鉴定出许多与生长性状有关的候选基因和基因家族，包括肌源性调节因子（MRF）家族成员（MYF5、MRF4、MYOD和MYOG）、肌细胞增强因子2基因家族（MEF2）、脂蛋白脂肪酶（LPL）和胰岛素样生长因子（IGFs）[2]。利用转录组测序技术揭示了黑番鸭骨骼肌不同生长发育阶段（胚胎期17、21、31、34 d及6月龄）的差异表达基因（Differentially expressed genes，DEGs）主要参与细胞生长、肌肉发育和细胞活动（连接、迁移、组装、分化和增殖）的过程，涉及到粘着斑、MAPK信号通路和肌动蛋白骨架调节信号通路，FAF1、RGS8、GRB10、SMYD3和TNNI2基 因是参与这些通路且可能诱导肌肉生长差异的候选基因[3]。

【本研究切入点】持续而密集的育种选择使得群体产生表型变化和相应的遗传变化。相似遗传背景的两个群体之间的表型差异能够用基因组变异和转录组变化来解释。转录组测序技术可以测定给定样本的所有转录组集合，从整体水平研究基因的功能和结构，揭示特定生物学过程的分子机理。转录组测序技术已经广泛运用在鸡、猪和牛等畜禽动物的生长性能表型差异研究中[4-7]，但采用该技术以番鸭生长性能差异为对象的研究较少。【拟解决的关键问题】本研究利用快速生长型番鸭和缓慢生长型番鸭的胸肌组织进行转录组测序，分析两种表型番鸭的基因表达差异，旨在挖掘影响番鸭生长性能的基因，为阐明生长性能相关差异的遗传机制提供理论基础。

1 材料与方法

1.1 试验番鸭来源及样品采集

供试番鸭12周龄6只，来自广东温氏思劳水禽试验场。试验番鸭群体为采样群体体重的两尾样，快速生长组番鸭体重分别为4.15、4.12、4.11 kg，对应样本编号为F_1、F_2和F_3；缓慢生长组番鸭体重分别为2.52、2.50、2.50 kg，对应样本编号为S_1、S_2和S_3。分别取番鸭左侧胸肌组织，称重后立即放入液氮中速冻，然后置于-80℃冰箱保存。

1.2 总RNA抽提与质检

利用氯仿法对番鸭胸肌组织总RNA进行抽提。总RNA抽提后，首先利用1.0%琼脂糖凝胶电泳分析RNA的降解程度以及是否有DNA污染，然后利用Nano Photometer Spectrophotometer（IMPLEN，美国）检测RNA纯度（OD260/OD280比值），再利用Qubit 2.0 Flurometer（Life Technologies，美国）对RNA浓度进行精确定量，最后利用Agilent Bioanalyzer 2100（Agilent Technologies，美国）精确检测RNA的完整性。总RNA质检保证28S : 18S ＞1.0，浓度≥500 ng/μL，OD260/OD280＞2.0，RNA完整值（RIN）≥8。

1.3 文库构建与测序

每个样品取总量为1.5 μg RNA用作制备测序样品。用附着有poly-T oligo的磁珠从总RNA中纯化mRNA，然后在5× NEBNext第一链合成反应缓冲液（NEB，美国）中进行裂解。用随机6聚体引物和M-MuLV逆转录酶（RNase H）合成第一链cDNA，随后使用DNA聚合酶I和RNase H合成cDNA的第二链。利用核酸外切酶/聚合酶活性将双链cDNA转化为平末端。DNA片段3'末端进行腺苷酸化后，将具有发夹环结构的NEBNext衔接子进行连接起来以便杂交。用AMPure XP系统（Beckman Coulter，Beverly，美国）纯化文库片段，以选择长度为150～200 bp的cDNA片段。然后在大小选择和衔接子连接后的cDNA中加入3 μL USER酶（NEB，美国），37 ℃下反应15 min，95 ℃下反应5 min。之后用Phusion高保真DNA聚合酶，Universal PCR引物和Index（X）引物进行扩增反应。最后，纯化PCR产物（AMPure XP系统），并在Agilent Bioanalyzer 2100系统上评估文库质量。

根据TruSeq PE Cluster Kit v3-cBot-HS（Illumina，美国）说明手册用cBot Cluster生成系统进行样本聚类，然后用Illumina Hiseq（Illumina，美国）进行测序。

1.4 测序数据分析

使用FASTQC软件对原始数据进行质量控制，剔除只含有测序接头序列的reads、含N比例大于0.1%的reads和低质量reads（质量值Q≤ 20的碱基数占整条read 50%以上），由此得到高质量序列数据（clean reads）。由于目前未见公开的番鸭全基因组数据，因此采用无参考基因组分析手段进行分析。本研究使用Trinity软件对clean reads进行拼接以获取参考序列[8]，并借助Corst软件将比对上转录本的reads数和表达模式对转录本进行聚类[9]。将Trinity软件拼接得到的转录本序列作为后续分析的参考序列，使用RSEM软件将clean reads比对到拼接的参考基因组[10]。使用FPKM（expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced，每百万碱基对的转录序列每千碱基的预期片段数）法计算基因表达量[11]，采用 DESeq 筛 选Padj＜ 0.05 和 |log2（fold change）|＞ 2的DEGs[12]，筛选过程中剔除组内表达read count值为0的DEGs。对Corset层次聚类后得到的最长Cluster序列在NR、KOG、GO、Pfam和NT、KEGG、Swiss-Prot数据库进行基因功能注释，进一步进行差异表达分析。

2 结果与分析

2.1 胸肌组织转录组测序分析

利用Trinity软件对clean reads进行拼接，其中200～500 bp的转录本数量最多，与该长度基因数差异最大（表1），不小于总长50%（N50）、90%（N90）的拼接转录本分别有2 041和297个，对应基因有2 846和625个（表2）。所有检测样品测序质量良好，Q20在96.83%～96.95%之间，Q30在92.05%～92.31%之间，GC含量为50.06%～50.88%，样本间基因表达相关性均达到0.69以上（表3），满足分析要求。将Trinity拼接得到的转录组作为参考序列，采用RSEM软件将每个样品的clean reads往参考序列上做比对，平均比对上的reads有42 126 585个，比对率为84.35%（表4）。

表1 拼接转录本频数分布情况Table 1 Splicing transcript frequency distribution

表2 拼接转录本长度分布情况Table 2 Splicing transcript length distribution

表3 样品间基因表达量相关性（R2）分析Table 3 Correlation analysis of gene expression between samples（R2）

表4 测序样本序列质量和比对信息统计Table 4 Statistical analysis for sequence quality and alignment information of RNA-seq libraries

2.2 基因功能注释

为了获取全面的基因功能信息，对Corset层次聚类后得到的最长Cluster序列在七大数据库进行基因功能注释。其中NT数据库注释成功率最高（58.93%），注释基因数为89 860个，NR数据库注释基因数48 048个、占31.51%，KO数据库注释基因数22 010个、占14.43%，SwissProt数据库注释基因数39 002个、占25.57%，Pfam数据库注释基因数45 199个、占29.64%，GO数据库注释基因数45 308个、占29.71%，KOG数据库注释基因数25 406个、占16.66%。NR、KOG、GO、Pfam和NT数据库共同注释到22 408个基因（图1）。

图1 NR、NT、Pfam、GO和KOG数据库的基因功能注释结果比较Fig. 1 Comparison of gene function annotation results of databases NR, NT, Pfam, GO and KOG

2.3 近缘物种相似性

通过与NR数据库进行比对注释，可以获取供试样本基因序列与近缘物种基因序列的相似性以及试验样本的基因功能信息。结果显示，供试番鸭与绿头鸭的基因组亲缘关系最近，相似性达46.7%，其次是鸡（6.4%）和金雕（5.9%），与野生火鸡和白头海雕的相似性较低，分别为3.0%和2.7%。

2.4 基因差异表达分析

使用FPKM法计算基因表达量，然后进行DEGs筛选（差异显著标准为Padj＜ 0.05）。本研究共得到186个DEGs，相比于生长缓慢组番鸭，快速生长组番鸭共有49个DEGs表达量出现上调，137个DEGs表达量出现下调（图2）。在这些 DEGs中，MAPK8IP3、PRKAG2和PPP2R3C在快速生长组番鸭胸肌组织中显著下调（表5），可能是影响番鸭生长速度的候选基因。

图2 差异表达基因火山图Fig. 2 Volcano diagram of DEGs

2.5 差异表达基因功能分析

基于NR和Pfam数据库对得到的186个DEGs进行功能注释，成功注释了142个。使用Blast2 GO软件对142个DEGs进行GO条目的功能富集分析[13]，DEGs富集到1 531个GO条目，集中分布在生物学过程、细胞成分和分子功能三大类。选取P值最小的前30个GO条目进行分析，发现DEGs主要富集共生体与宿主的粘附相关生物过程。此外，这些DEGs还富集到L-半胱氨酸代谢过程、蛋白质加工与成熟、碳水化合物运输和钙离子转运的正调控等条目，富集到这些条目的DEGs可能涉及到番鸭机体在生长过程中的能量代谢与相关蛋白的合成与运输等（图3）。

表5 不同生长速度番鸭的差异表达基因Table 5 DEGs in Muscovy ducks with different growth rates

2.6 差异表达基因通路分析

利用KAAS平台的KEGG自动注释服务对DEGs进行通路分析，共富集到140个通路（q＜1）。图4展示了DEGs富集程度前20的信号通路。DEGs可以被富集在磷脂酰肌醇三激酶-丝氨酸/苏氨酸激酶（PI3K-AKT）信号通路、丝裂原活化蛋白激酶（MAPK）信号通路和腺苷酸激活蛋白激酶（AMPK）信号通路，意味着这些通路可能与番鸭的生长相关。

3 讨论

快速生长型番鸭和慢速生长型番鸭具有不同的肌肉生长和发育表现，为了探究这些不同表现所涉及到的基因和信号通路，本研究对两种生长类型番鸭的胸肌组织进行了转录组测序分析。结果发现，快速生长型番鸭和慢速生长型番鸭相比存在186个显著DEGs，其中上调49个，下调137个。通过GO富集发现，这些DEGs主要富集在共生体与宿主的粘附相关生物学过程和一些在生长过程中与能量代谢以及相关蛋白合成与运输有关的条目。对DEGs进行KEGG分析，发现DEGs主要集中在磷脂酰肌醇三激酶-丝氨酸/苏氨酸激酶（PI3K-AKT）信号通路、丝裂原活化蛋白激酶（MAPK）信号通路和腺苷酸激活蛋白激酶（AMPK）信号通路。

图3 差异表达基因GO富集Fig. 3 GO enrichment analysis of DEGs

图4 差异表达基因KEGG富集Fig. 4 KEGG pathway enrichment analysis of DEGs

成年动物的肌肉发育主要靠肌纤维的质量增加而不是数量增加[14]，肌纤维通过两个途径来生长发育：第一，肌纤维在各种调控因子作用下通过增加蛋白质表达量来增加肌纤维的体积质量；第二，卫星细胞收到相关信号后可以分化为成肌细胞，融合到肌纤维中增加肌纤维的体积和质量。MAPK信号通路、PI3K-AKT信号通路和AMPK信号通路与肌肉的发育存在不可分割的联系。

MAPK是丝氨酸/苏氨酸蛋白激酶，可磷酸化底物的丝氨酸/苏氨酸残基。MAPK通路是联系细胞外刺激（包括生长因子和激素）与细胞内基因表达调控的桥梁，可调节多种细胞活动包括增殖、分化、凋亡、存活、炎症和先天免疫[15]。MAPK通路扮演着肌肉干细胞命运调控者的重要角色。p38 MAPK是MAPK的一个亚类，在肌肉分化中发挥至关重要的作用[16]。在成肌细胞分化过程中，p38会被持续激活，调节参与肌肉生成转录和表观遗传调控的许多参与因子的表达或活性[17]。p38α/βMAPKs通过直接将 MEF2 和E47磷酸化来促进MEF2转录活性和MyoD / E47异二聚体形成[18-19]，增强肌肉生成素基因的表达。蛋白精氨酸甲基转移酶（Prmt7）能够通过精氨酸残基70处p38αMAPK的甲基化促进MyoD介导的成肌细胞分化[20]。NF-κB信号通路也是完成成肌程序的必需通路，p38 MAPK通路和NF-κB通路存在串扰，NF-κB充当p38 MAPK的下游效应物诱导成肌细胞分化[21]。

PI3Ks是细胞内脂质激酶的独特家族，可磷酸化磷脂酰肌醇，而磷脂酰肌醇是真核细胞膜的重要组成成分。AKT包含3个结构域：pH结构域、中间激酶和调节性羧基末端结构域。PI3K将膜磷脂磷脂酰肌醇-4, 5-二磷酸磷酸化为磷脂酰肌醇-3, 4, 5-三磷酸，为磷脂酸肌醇依赖性激酶1（PDK-1）和AKT磷酸化提供膜脂质结合位点，随后AKT和PDK-1磷酸化激活[22]。AKT参与多种细胞过程，包括增殖、代谢和细胞大小调节[23]。PI3Ks/AKT信号通路通过在机体生长和关键细胞过程（例如葡萄糖稳态、脂质代谢、蛋白质合成以及细胞增殖和存活）中介导生长因子信号，在细胞生理学中发挥重要作用。在肌肉生长和再生中，PI3Ks/AKT信号通路不可或缺。胰岛素样生长因子1（IGF-1）可通过PI3Ks/AKT信号通路促进骨骼肌蛋白质合成，促进肌肥大，该通路还能抑制FOXO，并且抑制E泛素连接酶的转录，从而抑制肌肉萎缩[24-25]。

AMPK是细胞内能量状态的中央传感器[26]。运动可以通过AMPK介导的信号通路增强骨骼肌葡萄糖的摄取和脂肪酸的氧化，以便于为肌肉收缩提供充足的能量[27]。在骨骼肌生长发育中AMPK具有促进蛋白质分解和抗蛋白质合成的作用。AMPK响应食物缺乏、运动和其他压力源刺激，介导不同类型的自噬途径，促进肌肉萎缩并清除异常线粒体[28]。AMPK的激活通过抑制雷帕霉素复合物1（mTORC1）活性从而抑制蛋白质合成[29]。

MAPK通路和PI3K/AKT通路在肌肉分化调控过程中相互依赖[30]，AMPK信号通路与肌细胞能量代谢和蛋白质合成与分解有关。MAPK8IP3、PRKAG2、PPP2R3C和HSPA2等番鸭生长速度潜在调控基因可能介导上述3条信号通路，发挥肌肉发育调节作用。

4 结论

肌肉大小是由总体肌纤维数量和单个肌纤维体积共同决定的。本研究利用转录组测序技术分析了不同生长速度番鸭的胸肌组织基因表达差异，筛选出的DEGs富集到丝裂原活化蛋白激酶（MAPK）信号通路、磷脂酰肌醇三激酶-丝氨酸/苏氨酸激酶（PI3K-Akt）信号通路和腺苷酸激活蛋白激酶（AMPK）信号通路，三者与肌肉生长代谢调控密切相关。DEGs（MAPK8IP3、PRKAG2和PPP2R3C等）可能通过MAPK信号通路激活卫星细胞增殖、分化、迁移并融合到肌纤维中促进肌肉生长，通过PI3K-AKT信号通路促进肌肉蛋白质合成，调控AMPK信号通路影响肌肉萎缩，此3条信号通路协调配合调控番鸭肌纤维的生长发育，导致同种番鸭表现出不同的生长速度表型。