轉錄組測序揭示不同生長速度番鴨胸肌組織的基因表達差異

張思雨，朱煒健，季從亮，郭利金，鄭 茗，徐海平，聶慶華

（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]，但采用該技術以番鴨生長性能差異為對象的研究較少?！緮M解決的關鍵問題】本研究利用快速生長型番鴨和緩慢生長型番鴨的胸肌組織進行轉錄組測序，分析兩種表型番鴨的基因表達差異，旨在挖掘影響番鴨生長性能的基因，為闡明生長性能相關差異的遺傳機制提供理論基礎。

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進行抽提?？俁NA抽提后，首先利用1.0%瓊脂糖凝膠電泳分析RNA的降解程度以及是否有DNA污染，然后利用Nano Photometer Spectrophotometer（IMPLEN，美國）檢測RNA純度（OD260/OD280比值），再利用Qubit 2.0 Flurometer（Life Technologies，美國）對RNA濃度進行精確定量，最后利用Agilent Bioanalyzer 2100（Agilent Technologies，美國）精確檢測RNA的完整性?？俁NA質檢保證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條信號通路協調配合調控番鴨肌纖維的生長發育，導致同種番鴨表現出不同的生長速度表型。