轉(zhuǎn)錄組測序揭示不同生長速度番鴨胸肌組織的基因表達(dá)差異

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

（1. 華南農(nóng)業(yè)大學(xué)動物科學(xué)學(xué)院/農(nóng)業(yè)農(nóng)村部雞遺傳育種與繁殖重點實驗室，廣東 廣州 510642；2. 溫氏食品集團股份有限公司，廣東 云浮 527400）

【研究意義】畜禽的生長速度對于滿足全球日益增長的蛋白質(zhì)需求至關(guān)重要。我國生產(chǎn)的禽肉中，鴨肉產(chǎn)量僅次于雞肉，提高鴨的生長性能可以帶來極高的經(jīng)濟價值。通過基因選擇育種來提高飼料轉(zhuǎn)化效率和肌肉發(fā)育質(zhì)量是一種卓有成效的方法，因此，研究調(diào)控肉鴨生長發(fā)育速度的相關(guān)基因具有很高的研究價值和應(yīng)用價值。番鴨（Muscovy duck）相較于一般的家鴨具有體型大、產(chǎn)肉量大、飼料報酬高、肉質(zhì)性能好和瘦肉率高等特點，是廣泛養(yǎng)殖的肉鴨品種。然而，同一品種番鴨的生長速度可以表現(xiàn)出很大的差異，其中涉及到的遺傳機制目前尚不明確，有待進一步研究。

【前人研究進展】畜禽肉的產(chǎn)量直接關(guān)系到養(yǎng)殖行業(yè)的經(jīng)濟效益，因此肌肉發(fā)育調(diào)控的研究一直是畜牧業(yè)研究的重點。家禽不同品種所表現(xiàn)的生長性能差異為肌肉發(fā)育調(diào)控相關(guān)研究提供了理想的模型。運用轉(zhuǎn)錄組測序技術(shù)（RNA-Seq）對生長性能表現(xiàn)出差異的品種進行測序分析，有利于挖掘與生長性狀相關(guān)的調(diào)控基因，分析其涉及的生物功能和調(diào)控通路[1]。迄今為止，在鴨中已經(jīng)鑒定出許多與生長性狀有關(guān)的候選基因和基因家族，包括肌源性調(diào)節(jié)因子（MRF）家族成員（MYF5、MRF4、MYOD和MYOG）、肌細(xì)胞增強因子2基因家族（MEF2）、脂蛋白脂肪酶（LPL）和胰島素樣生長因子（IGFs）[2]。利用轉(zhuǎn)錄組測序技術(shù)揭示了黑番鴨骨骼肌不同生長發(fā)育階段（胚胎期17、21、31、34 d及6月齡）的差異表達(dá)基因（Differentially expressed genes，DEGs）主要參與細(xì)胞生長、肌肉發(fā)育和細(xì)胞活動（連接、遷移、組裝、分化和增殖）的過程，涉及到粘著斑、MAPK信號通路和肌動蛋白骨架調(diào)節(jié)信號通路，F(xiàn)AF1、RGS8、GRB10、SMYD3和TNNI2基 因是參與這些通路且可能誘導(dǎo)肌肉生長差異的候選基因[3]。

【本研究切入點】持續(xù)而密集的育種選擇使得群體產(chǎn)生表型變化和相應(yīng)的遺傳變化。相似遺傳背景的兩個群體之間的表型差異能夠用基因組變異和轉(zhuǎn)錄組變化來解釋。轉(zhuǎn)錄組測序技術(shù)可以測定給定樣本的所有轉(zhuǎn)錄組集合，從整體水平研究基因的功能和結(jié)構(gòu)，揭示特定生物學(xué)過程的分子機理。轉(zhuǎn)錄組測序技術(shù)已經(jīng)廣泛運用在雞、豬和牛等畜禽動物的生長性能表型差異研究中[4-7]，但采用該技術(shù)以番鴨生長性能差異為對象的研究較少?！緮M解決的關(guān)鍵問題】本研究利用快速生長型番鴨和緩慢生長型番鴨的胸肌組織進行轉(zhuǎn)錄組測序，分析兩種表型番鴨的基因表達(dá)差異，旨在挖掘影響番鴨生長性能的基因，為闡明生長性能相關(guān)差異的遺傳機制提供理論基礎(chǔ)。

1 材料與方法

1.1 試驗番鴨來源及樣品采集

供試番鴨12周齡6只，來自廣東溫氏思勞水禽試驗場。試驗番鴨群體為采樣群體體重的兩尾樣，快速生長組番鴨體重分別為4.15、4.12、4.11 kg，對應(yīng)樣本編號為F_1、F_2和F_3；緩慢生長組番鴨體重分別為2.52、2.50、2.50 kg，對應(yīng)樣本編號為S_1、S_2和S_3。分別取番鴨左側(cè)胸肌組織，稱重后立即放入液氮中速凍，然后置于-80℃冰箱保存。

1.2 總RNA抽提與質(zhì)檢

利用氯仿法對番鴨胸肌組織總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的完整性?？俁NA質(zhì)檢保證28S : 18S ＞1.0，濃度≥500 ng/μL，OD260/OD280＞2.0，RNA完整值（RIN）≥8。

1.3 文庫構(gòu)建與測序

每個樣品取總量為1.5 μg RNA用作制備測序樣品。用附著有poly-T oligo的磁珠從總RNA中純化mRNA，然后在5× NEBNext第一鏈合成反應(yīng)緩沖液（NEB，美國）中進行裂解。用隨機6聚體引物和M-MuLV逆轉(zhuǎn)錄酶（RNase H）合成第一鏈cDNA，隨后使用DNA聚合酶I和RNase H合成cDNA的第二鏈。利用核酸外切酶/聚合酶活性將雙鏈cDNA轉(zhuǎn)化為平末端。DNA片段3'末端進行腺苷酸化后，將具有發(fā)夾環(huán)結(jié)構(gòu)的NEBNext銜接子進行連接起來以便雜交。用AMPure XP系統(tǒng)（Beckman Coulter，Beverly，美國）純化文庫片段，以選擇長度為150～200 bp的cDNA片段。然后在大小選擇和銜接子連接后的cDNA中加入3 μL USER酶（NEB，美國），37 ℃下反應(yīng)15 min，95 ℃下反應(yīng)5 min。之后用Phusion高保真DNA聚合酶，Universal PCR引物和Index（X）引物進行擴增反應(yīng)。最后，純化PCR產(chǎn)物（AMPure XP系統(tǒng)），并在Agilent Bioanalyzer 2100系統(tǒng)上評估文庫質(zhì)量。

根據(jù)TruSeq PE Cluster Kit v3-cBot-HS（Illumina，美國）說明手冊用cBot Cluster生成系統(tǒng)進行樣本聚類，然后用Illumina Hiseq（Illumina，美國）進行測序。

1.4 測序數(shù)據(jù)分析

使用FASTQC軟件對原始數(shù)據(jù)進行質(zhì)量控制，剔除只含有測序接頭序列的reads、含N比例大于0.1%的reads和低質(zhì)量reads（質(zhì)量值Q≤ 20的堿基數(shù)占整條read 50%以上），由此得到高質(zhì)量序列數(shù)據(jù)（clean reads）。由于目前未見公開的番鴨全基因組數(shù)據(jù)，因此采用無參考基因組分析手段進行分析。本研究使用Trinity軟件對clean reads進行拼接以獲取參考序列[8]，并借助Corst軟件將比對上轉(zhuǎn)錄本的reads數(shù)和表達(dá)模式對轉(zhuǎn)錄本進行聚類[9]。將Trinity軟件拼接得到的轉(zhuǎn)錄本序列作為后續(xù)分析的參考序列，使用RSEM軟件將clean reads比對到拼接的參考基因組[10]。使用FPKM（expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced，每百萬堿基對的轉(zhuǎn)錄序列每千堿基的預(yù)期片段數(shù)）法計算基因表達(dá)量[11]，采用 DESeq 篩 選Padj＜ 0.05 和 |log2（fold change）|＞ 2的DEGs[12]，篩選過程中剔除組內(nèi)表達(dá)read count值為0的DEGs。對Corset層次聚類后得到的最長Cluster序列在NR、KOG、GO、Pfam和NT、KEGG、Swiss-Prot數(shù)據(jù)庫進行基因功能注釋，進一步進行差異表達(dá)分析。

2 結(jié)果與分析

2.1 胸肌組織轉(zhuǎn)錄組測序分析

利用Trinity軟件對clean reads進行拼接，其中200～500 bp的轉(zhuǎn)錄本數(shù)量最多，與該長度基因數(shù)差異最大（表1），不小于總長50%（N50）、90%（N90）的拼接轉(zhuǎn)錄本分別有2 041和297個，對應(yīng)基因有2 846和625個（表2）。所有檢測樣品測序質(zhì)量良好，Q20在96.83%～96.95%之間，Q30在92.05%～92.31%之間，GC含量為50.06%～50.88%，樣本間基因表達(dá)相關(guān)性均達(dá)到0.69以上（表3），滿足分析要求。將Trinity拼接得到的轉(zhuǎn)錄組作為參考序列，采用RSEM軟件將每個樣品的clean reads往參考序列上做比對，平均比對上的reads有42 126 585個，比對率為84.35%（表4）。

表1 拼接轉(zhuǎn)錄本頻數(shù)分布情況Table 1 Splicing transcript frequency distribution

表2 拼接轉(zhuǎn)錄本長度分布情況Table 2 Splicing transcript length distribution

表3 樣品間基因表達(dá)量相關(guān)性（R2）分析Table 3 Correlation analysis of gene expression between samples（R2）

表4 測序樣本序列質(zhì)量和比對信息統(tǒng)計Table 4 Statistical analysis for sequence quality and alignment information of RNA-seq libraries

2.2 基因功能注釋

為了獲取全面的基因功能信息，對Corset層次聚類后得到的最長Cluster序列在七大數(shù)據(jù)庫進行基因功能注釋。其中NT數(shù)據(jù)庫注釋成功率最高（58.93%），注釋基因數(shù)為89 860個，NR數(shù)據(jù)庫注釋基因數(shù)48 048個、占31.51%，KO數(shù)據(jù)庫注釋基因數(shù)22 010個、占14.43%，SwissProt數(shù)據(jù)庫注釋基因數(shù)39 002個、占25.57%，Pfam數(shù)據(jù)庫注釋基因數(shù)45 199個、占29.64%，GO數(shù)據(jù)庫注釋基因數(shù)45 308個、占29.71%，KOG數(shù)據(jù)庫注釋基因數(shù)25 406個、占16.66%。NR、KOG、GO、Pfam和NT數(shù)據(jù)庫共同注釋到22 408個基因（圖1）。

圖1 NR、NT、Pfam、GO和KOG數(shù)據(jù)庫的基因功能注釋結(jié)果比較Fig. 1 Comparison of gene function annotation results of databases NR, NT, Pfam, GO and KOG

2.3 近緣物種相似性

通過與NR數(shù)據(jù)庫進行比對注釋，可以獲取供試樣本基因序列與近緣物種基因序列的相似性以及試驗樣本的基因功能信息。結(jié)果顯示，供試番鴨與綠頭鴨的基因組親緣關(guān)系最近，相似性達(dá)46.7%，其次是雞（6.4%）和金雕（5.9%），與野生火雞和白頭海雕的相似性較低，分別為3.0%和2.7%。

2.4 基因差異表達(dá)分析

使用FPKM法計算基因表達(dá)量，然后進行DEGs篩選（差異顯著標(biāo)準(zhǔn)為Padj＜ 0.05）。本研究共得到186個DEGs，相比于生長緩慢組番鴨，快速生長組番鴨共有49個DEGs表達(dá)量出現(xiàn)上調(diào)，137個DEGs表達(dá)量出現(xiàn)下調(diào)（圖2）。在這些 DEGs中，MAPK8IP3、PRKAG2和PPP2R3C在快速生長組番鴨胸肌組織中顯著下調(diào)（表5），可能是影響番鴨生長速度的候選基因。

圖2 差異表達(dá)基因火山圖Fig. 2 Volcano diagram of DEGs

2.5 差異表達(dá)基因功能分析

基于NR和Pfam數(shù)據(jù)庫對得到的186個DEGs進行功能注釋，成功注釋了142個。使用Blast2 GO軟件對142個DEGs進行GO條目的功能富集分析[13]，DEGs富集到1 531個GO條目，集中分布在生物學(xué)過程、細(xì)胞成分和分子功能三大類。選取P值最小的前30個GO條目進行分析，發(fā)現(xiàn)DEGs主要富集共生體與宿主的粘附相關(guān)生物過程。此外，這些DEGs還富集到L-半胱氨酸代謝過程、蛋白質(zhì)加工與成熟、碳水化合物運輸和鈣離子轉(zhuǎn)運的正調(diào)控等條目，富集到這些條目的DEGs可能涉及到番鴨機體在生長過程中的能量代謝與相關(guān)蛋白的合成與運輸?shù)龋▓D3）。

表5 不同生長速度番鴨的差異表達(dá)基因Table 5 DEGs in Muscovy ducks with different growth rates

2.6 差異表達(dá)基因通路分析

利用KAAS平臺的KEGG自動注釋服務(wù)對DEGs進行通路分析，共富集到140個通路（q＜1）。圖4展示了DEGs富集程度前20的信號通路。DEGs可以被富集在磷脂酰肌醇三激酶-絲氨酸/蘇氨酸激酶（PI3K-AKT）信號通路、絲裂原活化蛋白激酶（MAPK）信號通路和腺苷酸激活蛋白激酶（AMPK）信號通路，意味著這些通路可能與番鴨的生長相關(guān)。

3 討論

快速生長型番鴨和慢速生長型番鴨具有不同的肌肉生長和發(fā)育表現(xiàn)，為了探究這些不同表現(xiàn)所涉及到的基因和信號通路，本研究對兩種生長類型番鴨的胸肌組織進行了轉(zhuǎn)錄組測序分析。結(jié)果發(fā)現(xiàn)，快速生長型番鴨和慢速生長型番鴨相比存在186個顯著DEGs，其中上調(diào)49個，下調(diào)137個。通過GO富集發(fā)現(xiàn)，這些DEGs主要富集在共生體與宿主的粘附相關(guān)生物學(xué)過程和一些在生長過程中與能量代謝以及相關(guān)蛋白合成與運輸有關(guān)的條目。對DEGs進行KEGG分析，發(fā)現(xiàn)DEGs主要集中在磷脂酰肌醇三激酶-絲氨酸/蘇氨酸激酶（PI3K-AKT）信號通路、絲裂原活化蛋白激酶（MAPK）信號通路和腺苷酸激活蛋白激酶（AMPK）信號通路。

圖3 差異表達(dá)基因GO富集Fig. 3 GO enrichment analysis of DEGs

圖4 差異表達(dá)基因KEGG富集Fig. 4 KEGG pathway enrichment analysis of DEGs

成年動物的肌肉發(fā)育主要靠肌纖維的質(zhì)量增加而不是數(shù)量增加[14]，肌纖維通過兩個途徑來生長發(fā)育：第一，肌纖維在各種調(diào)控因子作用下通過增加蛋白質(zhì)表達(dá)量來增加肌纖維的體積質(zhì)量；第二，衛(wèi)星細(xì)胞收到相關(guān)信號后可以分化為成肌細(xì)胞，融合到肌纖維中增加肌纖維的體積和質(zhì)量。MAPK信號通路、PI3K-AKT信號通路和AMPK信號通路與肌肉的發(fā)育存在不可分割的聯(lián)系。

MAPK是絲氨酸/蘇氨酸蛋白激酶，可磷酸化底物的絲氨酸/蘇氨酸殘基。MAPK通路是聯(lián)系細(xì)胞外刺激（包括生長因子和激素）與細(xì)胞內(nèi)基因表達(dá)調(diào)控的橋梁，可調(diào)節(jié)多種細(xì)胞活動包括增殖、分化、凋亡、存活、炎癥和先天免疫[15]。MAPK通路扮演著肌肉干細(xì)胞命運調(diào)控者的重要角色。p38 MAPK是MAPK的一個亞類，在肌肉分化中發(fā)揮至關(guān)重要的作用[16]。在成肌細(xì)胞分化過程中，p38會被持續(xù)激活，調(diào)節(jié)參與肌肉生成轉(zhuǎn)錄和表觀遺傳調(diào)控的許多參與因子的表達(dá)或活性[17]。p38α/βMAPKs通過直接將 MEF2 和E47磷酸化來促進MEF2轉(zhuǎn)錄活性和MyoD / E47異二聚體形成[18-19]，增強肌肉生成素基因的表達(dá)。蛋白精氨酸甲基轉(zhuǎn)移酶（Prmt7）能夠通過精氨酸殘基70處p38αMAPK的甲基化促進MyoD介導(dǎo)的成肌細(xì)胞分化[20]。NF-κB信號通路也是完成成肌程序的必需通路，p38 MAPK通路和NF-κB通路存在串?dāng)_，NF-κB充當(dāng)p38 MAPK的下游效應(yīng)物誘導(dǎo)成肌細(xì)胞分化[21]。

PI3Ks是細(xì)胞內(nèi)脂質(zhì)激酶的獨特家族，可磷酸化磷脂酰肌醇，而磷脂酰肌醇是真核細(xì)胞膜的重要組成成分。AKT包含3個結(jié)構(gòu)域：pH結(jié)構(gòu)域、中間激酶和調(diào)節(jié)性羧基末端結(jié)構(gòu)域。PI3K將膜磷脂磷脂酰肌醇-4, 5-二磷酸磷酸化為磷脂酰肌醇-3, 4, 5-三磷酸，為磷脂酸肌醇依賴性激酶1（PDK-1）和AKT磷酸化提供膜脂質(zhì)結(jié)合位點，隨后AKT和PDK-1磷酸化激活[22]。AKT參與多種細(xì)胞過程，包括增殖、代謝和細(xì)胞大小調(diào)節(jié)[23]。PI3Ks/AKT信號通路通過在機體生長和關(guān)鍵細(xì)胞過程（例如葡萄糖穩(wěn)態(tài)、脂質(zhì)代謝、蛋白質(zhì)合成以及細(xì)胞增殖和存活）中介導(dǎo)生長因子信號，在細(xì)胞生理學(xué)中發(fā)揮重要作用。在肌肉生長和再生中，PI3Ks/AKT信號通路不可或缺。胰島素樣生長因子1（IGF-1）可通過PI3Ks/AKT信號通路促進骨骼肌蛋白質(zhì)合成，促進肌肥大，該通路還能抑制FOXO，并且抑制E泛素連接酶的轉(zhuǎn)錄，從而抑制肌肉萎縮[24-25]。

AMPK是細(xì)胞內(nèi)能量狀態(tài)的中央傳感器[26]。運動可以通過AMPK介導(dǎo)的信號通路增強骨骼肌葡萄糖的攝取和脂肪酸的氧化，以便于為肌肉收縮提供充足的能量[27]。在骨骼肌生長發(fā)育中AMPK具有促進蛋白質(zhì)分解和抗蛋白質(zhì)合成的作用。AMPK響應(yīng)食物缺乏、運動和其他壓力源刺激，介導(dǎo)不同類型的自噬途徑，促進肌肉萎縮并清除異常線粒體[28]。AMPK的激活通過抑制雷帕霉素復(fù)合物1（mTORC1）活性從而抑制蛋白質(zhì)合成[29]。

MAPK通路和PI3K/AKT通路在肌肉分化調(diào)控過程中相互依賴[30]，AMPK信號通路與肌細(xì)胞能量代謝和蛋白質(zhì)合成與分解有關(guān)。MAPK8IP3、PRKAG2、PPP2R3C和HSPA2等番鴨生長速度潛在調(diào)控基因可能介導(dǎo)上述3條信號通路，發(fā)揮肌肉發(fā)育調(diào)節(jié)作用。

4 結(jié)論

肌肉大小是由總體肌纖維數(shù)量和單個肌纖維體積共同決定的。本研究利用轉(zhuǎn)錄組測序技術(shù)分析了不同生長速度番鴨的胸肌組織基因表達(dá)差異，篩選出的DEGs富集到絲裂原活化蛋白激酶（MAPK）信號通路、磷脂酰肌醇三激酶-絲氨酸/蘇氨酸激酶（PI3K-Akt）信號通路和腺苷酸激活蛋白激酶（AMPK）信號通路，三者與肌肉生長代謝調(diào)控密切相關(guān)。DEGs（MAPK8IP3、PRKAG2和PPP2R3C等）可能通過MAPK信號通路激活衛(wèi)星細(xì)胞增殖、分化、遷移并融合到肌纖維中促進肌肉生長，通過PI3K-AKT信號通路促進肌肉蛋白質(zhì)合成，調(diào)控AMPK信號通路影響肌肉萎縮，此3條信號通路協(xié)調(diào)配合調(diào)控番鴨肌纖維的生長發(fā)育，導(dǎo)致同種番鴨表現(xiàn)出不同的生長速度表型。