朱雯 湯瑩瑩 孫昕旸 周明 張子軍 陳興勇

（安徽農(nóng)業(yè)大學動物科技學院，合肥 230036）

近年來，畜牧業(yè)的快速發(fā)展增加了畜禽飼料原料的需求量。與之相比，我國糧食產(chǎn)量漲幅相對較低。飼料資源尤其是優(yōu)質蛋白質飼料較為緊缺，如大豆每年進口8.0×107t，中美貿(mào)易戰(zhàn)和新發(fā)展格局下，蛋白資源緊缺威脅到飼料糧安全，非常規(guī)蛋白資源開發(fā)、低蛋白平衡日糧等成為問題解決方向［1］。然而，大部分研究表明，低蛋白飼糧會影響反芻動物生產(chǎn)性能［2-3］。肝臟是調節(jié)糖脂代謝的重要器官，瘤胃發(fā)酵以及飼糧來源的各種營養(yǎng)素均通過門靜脈進入肝臟［4］。由此可見，肝臟功能對反芻動物生產(chǎn)性能的影響起著至關重要的作用。目前，低蛋白日糧在肝臟中代謝調控機制仍不明確，現(xiàn)代組學技術發(fā)展迅速，可有效闡明這一問題。

轉錄水平的調控是動物生命活動調節(jié)的重要調控形式，轉錄組測序是了解動物基因總體表達情況的強大工具。目前高通量轉錄組測序技術已經(jīng)應用在營養(yǎng)和環(huán)境對反芻動物生長性能調控的影響機制研究中［5-6］。一大批涉及營養(yǎng)物質吸收及信號轉導的基因得到鑒定，對于研究動物生長性能調控的分子機制起到很大的推動作用。

團隊前期研究結果發(fā)現(xiàn)，飼糧粗蛋白水平較NRC推薦水平低2.80個百分點時顯著降低安徽白山羊的生長與生產(chǎn)性能［7］。本研究采用轉錄組測序（RNA-Seq）技術，分析低蛋白飼糧影響肝臟基因的表達特性，重點挖掘營養(yǎng)素代謝相關基因，以期揭示低蛋白飼糧影響動物生產(chǎn)性能的分子機制，為低蛋白飼糧的營養(yǎng)調控提供理論依據(jù)。

1 材料與方法

1.1 材料

安徽白山羊購自六安綠潔畜牧有限公司，RNA提取試劑盒、RNA純化試劑盒均購自 QIAGEN公司。

1.2 方法

1.2.1 試驗動物、日糧和樣品采集 試驗前對羊只進行檢疫，預試期進行驅蟲、防疫、定期消毒，保證羊舍干凈。試驗羊單欄飼喂，每天06：30和18：30各飼喂1次，確保槽內約10%的剩料。試驗期內所有羊只自由飲水和采食。試驗設計和動物飼養(yǎng)管理同本課題組前期研究報道［7］。簡述如下：本試驗采用單因子隨機分組設計，選取24只體重接近（16.0±1.35）kg的健康安徽白山羊公羊（3-4月齡），隨機分為2組，分別飼喂正常蛋白質水平（14.8%）飼糧（對照組，CK組）和80%蛋白水平（12.0%）飼糧（低蛋白組，LCP組）。基礎飼糧參照NRC（2007）［8］，體重20 kg、日增重150 g/d的公羊營養(yǎng)需要量自行配制。試驗飼糧為全混合顆粒飼料，直徑為2.5 mm，長度為8 mm，試驗飼糧組成及營養(yǎng)水平見表1。試驗結束后，每組隨機抽取3只羊屠宰，宰前禁食24 h，禁水2 h。屠宰后立即取肝臟組織放入無酶凍存管，保存于液氮備用。本實驗通過安徽農(nóng)業(yè)大學科學倫理委員會批準（批準號SYDWP20190600601）。

表1 試驗飼糧組成及營養(yǎng)水平（干物質基礎）Table 1 Ingredients and chemical composition of the experimental diets（dry matter basis）

1.2.2 RNA提取及轉錄組測序 通過Trizol法進行山羊肝臟組織總RNA提取，分別使用NanoPhotometer spectrophotometer和Qubit 2.0檢 測提 取RNA的純度和濃度，使用Aglient2100檢測RNA的完整性。所有樣品RNA提取并檢測合格后，用干冰保存并送至廣州基迪奧生物公司，應用Illumina Novaseq6000高通量測序平臺進行測序分析。

1.2.3 測序數(shù)據(jù)分析 對測序所獲得的原始數(shù)據(jù)，按照基迪奧生物公司的數(shù)據(jù)處理方法，進行原始數(shù)據(jù)的過濾，獲得高質量的clean reads。然后通過HISAT2軟件進行序列對比，獲得山羊肝臟樣本特異序列信息，并進行基因組定位分析。

1.2.4 差異表達基因（DEGs）的篩選以及功能富集分析 以FDR＜0.05 且 |log2FC|＞1為閾值進行差異表達基因的篩選，其中Fold Changes（FC）代表不同處理之間基因表達的差異。結合GO（gene ontology）和KEGG（kyoto encyclopedia of genes and genomes）數(shù)據(jù)庫對篩選獲得的差異表達基因進行功能注釋、GO及KEGG分析，P＜0.05 的通路為相關基因富集到的差異生物學通路。

1.2.5 基因的蛋白質互作網(wǎng)絡分析 采用STRING數(shù)據(jù)庫的蛋白質互作信息構建低蛋白飼糧調控相關差異表達基因的蛋白質互作網(wǎng)絡。

1.2.6 差異表達基因的驗證 為了對轉錄組測序的結果進行驗證，隨機挑選5個差異表達基因，利用Primer5.0軟件設計擴增引物（表2），采用qRT-PCR的方法，對其表達水平進行驗證。將各處理的樣品總 RNA 反轉錄合成 cDNA，并以此為模板，以 β-actin為內參，進行 qRT-PCR反應。qRT-PCR的反應體系、反應程序以及基因相對表達量的計算均參照Li等［5］方法進行。

表2 qRT-PCR反應所用引物Table 2 Primers used for qRT-PCR

2 結果

2.1 測序質量分析

采用高通量測序平臺，對低蛋白飼糧處理及對照組的肝臟樣品進行了RNA測序。結果顯示，共獲得46 719 458-62 169 938條reads，測序數(shù)據(jù)去除接頭以及低質量數(shù)據(jù)后分別獲得7 007 918 700-9 325 490 700個高質量數(shù)據(jù)（clean reads），所有樣品過濾后高質量序列數(shù)占原始下機序列的比例均大于95%，表明6個樣品構建的文庫質量較好；同時Q30的數(shù)值均大于90%，說明本次測序的質量非?？煽浚ū?）。

表3 各樣品測序質量Table 3 Sequencing data quality of each sample

進一步對樣品重復間數(shù)據(jù)的相關性分析，結果（圖1）表明，重復之間的皮爾森相關系數(shù)均大于0.8，最高則達到了0.965 6，說明數(shù)據(jù)之間的重現(xiàn)性較好，可用于下一步的分析。

圖1 樣品重復之間的相關性分析Fig.1 Pearson correlation analysis between the repeated samples

2.2 差異表達基因統(tǒng)計分析

將低蛋白飼糧組和對照組樣品的轉錄組數(shù)據(jù)進行對比分析，以FDR ＜ 0.05 以及|log2FC| ＞ 1為篩選條件，篩選差異表達基因。結果表明，低蛋白飼糧處理后，差異表達基因的數(shù)量達到了62個，其中，上調49個，下調13個。上調基因包含與蛋白質水解相關基因絲氨酸蛋白酶35（PRSS33）和PCOLCE2；能量代謝相關基因重組人碳酸酐酶14（CA14）和MDF1；脂肪酸代謝相關基因人源全長重組蛋白（ACSBG1）、脂肪酸合成酶（FASN）、?；o酶A硫酯酶（acyl-coenzyme A thioesterase 4）和LOC102184912。下調基因包含與糖代謝相關基因丙酮酸羧化酶（PC）和LOC108635023；胰島素代謝相關基因胰島素受體底物（IRS2）。

2.3 qRT-PCR驗證

為了驗證RNA測序結果的可靠性，隨機選取了PC、FASN、促生長因子（IGF1）、IRS2以及CYP4A18五個基因進行qRT-PCR擴增，并與其測序結果進行對比（圖2），從圖2可以看出，所選擇的5個基因在兩種方法中表達趨勢一致，說明本研究RNA-Seq測序結果較為可靠。

圖2 RNA-Seq數(shù)據(jù)（A）與qRT-PCR（B）基因表達水平的比較Fig.2 Comparison and verification of the quantitative results for selected genes from the RNA-Seq（A）and qRT-PCR analyses（B）

2.4 差異表達基因的GO功能分析

為了進一步了解差異表達基因的功能，通過GO數(shù)據(jù)庫對篩選到的差異表達基因進行基因功能富集分類，包括細胞組分（cellular component，CC）、生物過程（biological process，BP）以及分子功能（molecular function，MF）3個方面（圖3）。在BP過程中，差異表達基因所富集的GO的數(shù)目為21；MF及CC中分別富集到8和13個GO條目。在BP功能中，兩組間的差異表達基因在細胞過程（cellular process）和單機體進程（single-organism process）富集到的差異基因數(shù)目最多。在CC組分中，細胞器（organelle）和細胞結構（cell part）中所占比例最大。在MF功能中，差異基因在結合功能（binding）中所占的比例最高，催化活性（catalytic activity）次之。

圖3 差異表達基因的GO分析Fig.3 GO analysis of differentially expressed gene

在分子功能方面，上調的差異表達基因富集到3個GO條目（表4），其中，結合（binding）基因數(shù)目最多（20，40.8%），其次為催化活性（catalytic activity）基因（10，20.4%）；在生物過程方面，上調的差異表達基因富集到13個GO條目（表5），其中細胞進程（cellular process）基因數(shù)目最多（25，51.0%），其次為單體進程（single-organism process）基因（22，44.9%）和催化進程基因（22，44.9%）。

表4 上調差異表達基因的GO功能注釋分析（分子功能）Table 4 GO functional annotation of up-regulated expressed genes（molecular function）

表5 上調差異表達基因的 GO 功能注釋分析（生物過程）Table 5 GO functional annotation of up-regulated expressed genes（biological process）

在分子功能方面，下調的差異表達基因富集到3個GO條目（表6），其中，結合（binding）基因數(shù)目最多（6，46.2%），其次為催化活性（catalytic activity）基因（2，15.4%）；在生物過程方面，下調的差異表達基因富集到13個GO條目（表7），其中單機體進程（single-organism process）基因數(shù)目最多（6，46.2%），其次為生物調節(jié)（biological regulation）基因（5，38.5%）和細胞進程（cellular process）基因（5，38.5%）。

表6 下調差異表達基因的GO功能注釋分析（分子功能）Table 6 GO functional annotation of down-regulated expressed genes（molecular function）

表7 下調差異表達基因的 GO 功能注釋分析（生物過程）Table 7 GO functional annotation of down-regulated expressed genes（biological process）

2.5 差異基因KEGG通路注釋

通過對差異表達基因進行KEGG通路分析（P≤ 0.05），共獲得23條差異信號通路（圖4），文中列舉富集后最顯著的前5個通路以做示例說明（表8），包括：富集因子最高的為味覺傳導（taste transduction），脂肪酸合成（fatty acid biosynthesis），胰島素信號通路（insulin signaling pathway），視黃醇代謝（retinol metabolism）和癌癥轉錄失調（transcriptional misregulation in cancers）等。部分差異基因富集到與營養(yǎng)素代謝相關的脂肪酸合成以及胰島素信號通路中。如上調基因acyl-coenzyme A thioesterase 4和FASN富集到脂肪酸合成通路；上調基因ACSBG1和細胞色素（CYP4A1）富集到PPAR信號轉導通路；上調基因FASN和蛋白磷酸酶1調節(jié)亞基3B（PPP1R3B）以及下調基因IRS2富集到胰島素信號通路。

表8 低蛋白日糧飼喂條件下山羊肝臟組織富集顯著性差異的前5條通路Table 8 Top 5 pathways enriched from the differential expressed liver genes of low crude protein fed goat

圖4 差異表達基因的KEGG分析Fig. 4 KEGG analysis of differentially expressed genes

2.6 基因的蛋白互作網(wǎng)絡分析

蛋白質互作網(wǎng)絡分析結果發(fā)現(xiàn)（圖5），在62個差異表達基因中，有29個基因之間發(fā)生互作關系。越靠近互作圖中心的位置，表示該基因可能發(fā)揮的作用越強。由圖中可知，IGF1，胰島素樣生長因子結合蛋白1（IGFBP1），胰島素樣生長因子酸不穩(wěn)定亞 基（IGFALS），LOC108634620，LOC108634720，LOC108634721，C102182562，細胞因子信號轉導抑制因子（SOCS2），以及組蛋白去甲基化酶（KDM6A）與其它基因有較強的互作關系。

圖5 差異表達基因的PPI分析Fig.5 PPI analysis of differentially expressed genes

3 討論

RNA-Seq技術，主要研究某一時間段內特定細胞在一種功能狀態(tài)下轉錄出的所有RNA的總和，是研究基因表達圖譜以及蛋白質生物信息功能的重要方法［9］。本研究通過RNA-Seq技術對低蛋白飼糧和正常組飼糧飼喂山羊肝臟組織進行了差異表達基因篩選，并對其中的5個基因進行了qRT-PCR驗證，驗證結果與測序結果一致，表明測序結果可靠。本研究中，低蛋白飼糧顯著降低了PC和IGFBP1表達量，增加IGFALS，IGF1和IRS2基因表達量。丙酮酸羧化酶是調節(jié)肝臟糖異生的限速酶［10］。IGF1可促進組織攝取葡萄糖，調節(jié)糖原異生和糖酵解，促進蛋白質和脂肪合成［11］。IRS2是胰島素信號傳遞的重要介質，在血糖和胰島素敏感性的調控中發(fā)揮著重要作用［12］。IGFALS幾乎由肝臟產(chǎn)生并分泌到血液循環(huán)中，以高濃度存于血清中［13］，其可與IGF-1形成復合物穩(wěn)定血清中IGF-1的半衰期［14］。IGFBP1是胰島素樣生長因子結合蛋白家族的成員之一，可負向調控胰島素信號通路［15］。課題組前期研究發(fā)現(xiàn)，在低蛋白飼糧飼喂條件下，肉羊血清葡萄糖濃度顯著下降［7］?？梢姡偷鞍兹占Z可正向調控胰島素信號通路，抑制肝臟糖異生，最終影響肉羊血糖濃度。

本研究結合 GO數(shù)據(jù)庫對差異基因主要參與的生物學過程進行了分類，發(fā)現(xiàn)差異表達基因主要參與了結合和催化活性等生物學過程，表明飼糧蛋白水平對肝臟的生理和代謝有較大的影響。研究發(fā)現(xiàn)，低蛋白飼糧通過提高豬［16］和羊［17］骨骼肌細胞中FASN基因表達，增加骨骼肌肌間脂肪的含量。與前人研究結果相似，本研究結合KEGG富集分析，篩選出23條差異信號通路，其中2個上調表達基因（?；o酶A硫酯酶和FASN）參與脂肪酸合同通路。碳水化合物反應原件（ChoRE）是脂類合成的關鍵轉錄因子，可被高糖激活，誘導下游基因脂肪酸合成酶的轉錄和表達［18］。FASN是脂肪酸從頭合成的關鍵酶，可直接與ChoRE啟動子結合，促進葡萄糖轉化為脂類［19］。因此，在山羊肝臟中FASN表達量的上升提示低蛋白飼糧可增加肝臟脂肪的合成。

過氧化物酶體增殖物活化受體（peroxisome porliferator-activated receptor，PPAR）屬 II 型核受體超家族成員，是一類由配體激活的核轉錄因子，是調節(jié)脂質代謝、脂肪生成、胰島素敏感、炎癥反應、細胞生長和分化的重要因子［20］。KEGG分析得到PPAR信號通路的顯著富集，提示機體能量平衡受影響。PPAR目前發(fā)現(xiàn)有3種亞型，PPARα，PPARβ（δ）和 PPARγ，其在不同組織中表達。PPARα 主要在脂肪酸氧化較高的組織中表達，如肝臟、心臟、肌肉和棕色脂肪組織等，參與脂肪酸代謝、炎癥反應等［21］。上調基因CYP4A1和ACSBG1基因富集于PPAR信號通路。CYP4A1 參與環(huán)境化合物在體內的代謝過程，受PPARα調控，是脂肪酸β酸氧化過程中的關鍵酶。ACSBG1可促進細胞增殖，誘導細胞凋亡［22］。本研究中，CYP4A1和ACSBG1基因表達量在低蛋白飼糧組山羊肝臟組織中均上調，表明低蛋白飼糧可促進山羊肝臟細胞的增殖，誘導肝臟脂肪酸氧化，而過度的脂質氧化會對肝臟造成直接的毒性作用。因此，低蛋白飼糧對肝臟健康不利。通過PPI互作圖發(fā)現(xiàn)KDM6A與其它基因有較強的互作關系。KDM6A在包括轉錄在內的多種生物進程中發(fā)揮非常重要的作用［23］，KDM6A基因敲除會影響細胞分化［24］。本研究中，低蛋白飼糧顯著降低山羊肝臟組織中KDM6A基因的表達，同樣提示低蛋白飼糧對肝臟組織發(fā)育不利。

4 結論

在低蛋白飼糧飼喂條件下，山羊肝臟組織大量的基因發(fā)生了顯著的差異表達，部分基因與蛋白質代謝、能量及糖脂代謝相關。這些差異表達基因的功能主要集中于肝臟胰島素信號通路和脂質代謝調控。低蛋白飼糧促進肝臟脂肪酸的合成和脂質過氧化，抑制肝臟糖異生，不利于肝臟組織發(fā)育。