仔豬生長(zhǎng)發(fā)育關(guān)鍵基因篩選和功能分析

宋 陽(yáng) ，蘆春雪，殷淑潔 ，馮 靜 ，劉 歡 ，蘇 辰 ，周百靈，沈叢叢，陳清帥，黃平平*，于 雪*

（1.德州學(xué)院生物物理研究院/山東省生物物理重點(diǎn)實(shí)驗(yàn)室，山東 德州 253023；2.德州學(xué)院生命科學(xué)學(xué)院，山東 德州）

斷奶仔豬是指處在從依靠母豬乳汁到獨(dú)立依靠固體飼料的仔豬[1]。斷奶過(guò)渡期是仔豬生長(zhǎng)過(guò)程中最重要的階段之一，該階段仔豬生長(zhǎng)發(fā)育快，但消化系統(tǒng)尚未發(fā)育完善，因此容易出現(xiàn)采食量差、生長(zhǎng)性能下降和胃腸功能受損等嚴(yán)重后果[2]。過(guò)去研究者們主要通過(guò)對(duì)營(yíng)養(yǎng)、健康和管理等有關(guān)技術(shù)進(jìn)行改進(jìn)，盡量減少斷奶應(yīng)激對(duì)仔豬生長(zhǎng)的不利影響[3]。另外，骨骼肌和脂肪組織不僅是重要的營(yíng)養(yǎng)部位，更是生長(zhǎng)的關(guān)鍵組成部分，且能夠?qū)Ω鞣N脅迫進(jìn)行生理適應(yīng)，包括疾病和壓力等[4]。因此，揭示在仔豬背最長(zhǎng)肌和皮下脂肪組織生長(zhǎng)發(fā)育過(guò)程起關(guān)鍵作用的基因，將有助于改善家畜肉質(zhì)和產(chǎn)量，同時(shí)有利于治療和預(yù)防相關(guān)疾病。

近年來(lái)隨著高通量測(cè)序技術(shù)的發(fā)展，與仔豬生長(zhǎng)發(fā)育相關(guān)基因的挖掘越來(lái)越受到國(guó)內(nèi)外研究者們的關(guān)注[5]。如Pilcher等[5]和Chen等[6]分別基于豬Ensembl release 70和Ensembl release 93版本的基因組和注釋文件分析背最長(zhǎng)肌和皮下脂肪組織發(fā)育相關(guān)的關(guān)鍵基因。根據(jù) Ensembl官方網(wǎng)站，截止 2021年底豬基因組和注釋文件已經(jīng)更新至Ensembl release 104版本。因此本研究在上述工作基礎(chǔ)上從 GEO數(shù)據(jù)庫(kù)下載背最長(zhǎng)肌和皮下脂肪組織各18個(gè)轉(zhuǎn)錄組測(cè)序數(shù)據(jù)，采用Ensembl release 104版本的豬基因組文件和注釋文件進(jìn)行生物信息學(xué)分析，進(jìn)一步挖掘與這兩種組織發(fā)育相關(guān)的關(guān)鍵基因，對(duì)補(bǔ)充認(rèn)識(shí)仔豬生長(zhǎng)發(fā)育相關(guān)的分子機(jī)制具有重要意義。

1 材料與方法

1.1 材料

本實(shí)驗(yàn)原始數(shù)據(jù)來(lái)自美國(guó)國(guó)家生物技術(shù)中心（National Center for Biotechnology Information，NCBI）的基因表達(dá)數(shù)據(jù)庫(kù)（Gene Expression Omnibus, GEO），共36個(gè)樣本，ID為GSE 65983。其中，18個(gè)樣本來(lái)自豬的背最長(zhǎng)?。↙D M），18個(gè)樣本來(lái)自皮下脂肪（BF）組織。豬的基因組和注釋文件自 Ensembl數(shù)據(jù)庫(kù)下載（http://ftp.ensembl.org/pub/release-104/fasta/sus_scrofa/；http://ftp.ensembl.org/pub/release-104/gtf/sus_scrofa/）。

1.2 方法

1.2.1 原始數(shù)據(jù)質(zhì)控 為了獲得高質(zhì)量的 clean data以便增強(qiáng)后續(xù)數(shù)據(jù)分析的準(zhǔn)確性，需要對(duì)獲得的原始數(shù)據(jù)質(zhì)控。采用 Trimmomatic軟件將原始數(shù)據(jù)中包含的測(cè)序接頭、低質(zhì)量 reads以及 N（無(wú)法確定堿基信息的比例大于5 % 的reads）、較短序列進(jìn)行過(guò)濾去除，最終獲得高質(zhì)量 clean data。

1.2.2 參考序列比對(duì) 目前研究表明，Hisat2（http://daehwankimlab.github.io/hisat2/）相較于TopHat2等對(duì)比軟件具有運(yùn)行速度快、精確度較好等優(yōu)點(diǎn)。因此本研究使用Hisat2軟件進(jìn)行序列比對(duì)，為提高其對(duì)比速度，實(shí)驗(yàn)前需構(gòu)建其基因的索引（index），具體操作如下：（1）下載豬（sus_scrofa）基因組文件。（2）索引構(gòu)建（hisat 2-build–p 16 genome.fa genome）。（3）進(jìn)行 Hisat2比對(duì)。比對(duì)之后利用 samtools軟件進(jìn)行數(shù)據(jù)格式轉(zhuǎn)換，將sam文件轉(zhuǎn)為排好序的bam文件。

1.2.3 轉(zhuǎn)錄本拼接、定量及差異表達(dá)分析 分別采用 Stringtie軟件組裝轉(zhuǎn)錄本、Samtools軟件對(duì)所有轉(zhuǎn)錄本定量，轉(zhuǎn)錄本的 TPM 值至少在一個(gè)樣本中大于 1被認(rèn)為是表達(dá)的。將 36個(gè)樣本根據(jù)組織部位（背最長(zhǎng)肌和皮下脂肪組織）分為 2組，然后采用 DEseq2軟件對(duì)所有轉(zhuǎn)錄本進(jìn)行差異表達(dá)分析。滿足以下條件的轉(zhuǎn)錄本認(rèn)為是差異表達(dá)的：在兩組樣本中∣log2FC∣＞1且 padj.＜0.05。

1.2.4 功能富集分析 為篩選與背最長(zhǎng)肌和皮下脂肪組織發(fā)育相關(guān)的關(guān)鍵基因，本研究利用DAVID對(duì)獲得的差異表達(dá)的 mRNA進(jìn)行基因本體論（Gene Oncology, GO）功能注釋和京都基因與基因組百科全書（Kyoto Encyclopedia of Genes and Genomes, KEGG）通路富集分析。

2 結(jié)果與分析

2.1 轉(zhuǎn)錄本數(shù)據(jù)統(tǒng)計(jì)

轉(zhuǎn)錄本的 TPM 值至少在一個(gè)樣本中大于 1認(rèn)為是表達(dá)的，36個(gè)樣本共獲得 79 800個(gè)轉(zhuǎn)錄本，其中包括 29 118個(gè) mRNA（11 858種基因），2 105個(gè)已知的lncRNA（2種基因，僅6個(gè)轉(zhuǎn)錄本含基因名），以及 47 142個(gè)新的轉(zhuǎn)錄本（表1）。

表1 轉(zhuǎn)錄本總體情況

2.2 差異表達(dá)分析

共獲得 34 873個(gè)差異表達(dá)轉(zhuǎn)錄本，其中9 671個(gè)轉(zhuǎn)錄本在背最長(zhǎng)肌組織上調(diào)；25 202個(gè)轉(zhuǎn)錄本在背最長(zhǎng)肌組織下調(diào)。另外包括12 536個(gè)差異表達(dá)的mRNA，3 352個(gè)mRNA在背最長(zhǎng)肌組織上調(diào)，9 184個(gè) mRNA在背最長(zhǎng)肌組織下調(diào)（表2）。

表2 轉(zhuǎn)錄本差異表達(dá)情況

2.3 差異表達(dá)mRNA GO功能富集分析

將獲得的在兩種組織中差異表達(dá)的mRNA進(jìn)行GO富集分析，結(jié)果分別獲得168個(gè)BP（biological process, BP）通路、139個(gè) CC（cellular component, CC）通路和 122個(gè) MF（molecular function，MF）通路（P-value < 0.05），圖1a、b、c中分別展示位于前 20位的 BP、CC和 MF通路。在 BP通路中，與肌肉和脂肪組織相關(guān)的通路包括肌節(jié)組織（sarcomere organization）、肌動(dòng)蛋白細(xì)胞骨架組織（actin cytoskeleton organization）和細(xì)胞遷移的積極調(diào)解（positive regulation of cell migration）等。61個(gè)mRNA轉(zhuǎn)錄本富集在肌節(jié)組織，如基因MYPN和TNNT1；124個(gè)mRNA轉(zhuǎn)錄本富集在肌動(dòng)蛋白細(xì)胞骨架組織，如基因ARHGAP26和PDLIM3；150個(gè)mRNA轉(zhuǎn)錄本富集在細(xì)胞遷移的積極調(diào)解，如基因PECAM1。在 CC通路中，與肌肉和脂肪組織相關(guān)的通路包括絲狀肌動(dòng)蛋白（filamentous actin）、肌動(dòng)蛋白絲（actin filament）和脂質(zhì)顆粒（lipid particle）等。38個(gè) mRNA轉(zhuǎn)錄本富集在絲狀肌動(dòng)蛋白，如基因PDLIM3；82個(gè) mRNA轉(zhuǎn)錄本富集在肌動(dòng)蛋白絲，如基因TPM4；83個(gè)mRNA轉(zhuǎn)錄本富集在脂質(zhì)顆粒，如基因PLIN1。在 MF通路中，與肌肉和脂肪組織相關(guān)的通路包括磷脂酰肌醇結(jié)合（phosphatidylinositol binding）、脂質(zhì)結(jié)合（lipid binding）和對(duì)旋肌球蛋白結(jié)合（tropomyosin binding）等。150個(gè)mRNA轉(zhuǎn)錄本富集在磷脂酰肌醇結(jié)合，94個(gè) mRNA轉(zhuǎn)錄本富集在脂質(zhì)結(jié)合，如基因APOA2；31個(gè)mRNA轉(zhuǎn)錄本富集在對(duì)旋肌球蛋白結(jié)合，如基因TNNT2和TMOD4。

圖1 差異表達(dá)mRNA功能富集分析結(jié)果

2.4 差異表達(dá)mRNA KEGG功能富集分析

通過(guò)KEGG信號(hào)通路分析共獲得187個(gè)通路（P-value < 0.05），圖1d僅展示位居前20位的通路。其中肌動(dòng)蛋白細(xì)胞骨架調(diào)節(jié)（Regulation of actin cytoskeleton）與肌肉組織相關(guān)。203個(gè)mRNA轉(zhuǎn)錄本富集在肌動(dòng)蛋白細(xì)胞骨架調(diào)節(jié)，如基因CFL和FGFR4。另外，脂肪細(xì)胞因子信號(hào)通路（Adipocytokine signaling pathway）、鞘脂信號(hào)通路（Sphingolipid signaling pathway）、脂肪酸代謝（Fatty acid metabolism）等通路與脂肪組織發(fā)育相關(guān)。73個(gè) mRNA轉(zhuǎn)錄本富集在脂肪細(xì)胞因子信號(hào)通路，如基因AKT2和PCK1；113個(gè)mRNA轉(zhuǎn)錄本富集在鞘脂信號(hào)通路，另外，63個(gè)mRNA轉(zhuǎn)錄本富集在脂肪酸代謝結(jié)合，如基因SCD5和FASN。

3 討論

在豬生長(zhǎng)發(fā)育的一系列階段中，斷奶會(huì)對(duì)豬的生長(zhǎng)發(fā)育產(chǎn)生不良影響。骨骼肌和脂肪組織是豬的營(yíng)養(yǎng)部位，在豬的發(fā)育中發(fā)揮重要的作用，同時(shí)也是仔豬生長(zhǎng)的關(guān)鍵部分，骨骼肌與皮下脂肪組織能夠使仔豬對(duì)各種疾病和壓力進(jìn)行生理適應(yīng)。皮下脂肪組織主要與能量代謝途徑有關(guān)，在豬不同的生長(zhǎng)階段，皮下脂肪組織的生長(zhǎng)狀況亦不同[7]。近年來(lái)與豬肌肉和脂肪組織發(fā)育相關(guān)的基因及其功能逐漸被揭示，如Wang等（2015）[8]對(duì)不同品種豬的背最長(zhǎng)肌組織進(jìn)行高通量測(cè)序和生物信息學(xué)分析，發(fā)現(xiàn)CAV2、MYOZ2和FRZB與肌肉生長(zhǎng)緊密相關(guān)，而FASN、SCD和ADORA1與脂質(zhì)沉積緊密相關(guān)。Li等（2016）應(yīng)用 RNA測(cè)序鑒定皖南花豬與約克夏豬背最長(zhǎng)肌差異表達(dá)基因并進(jìn)行功能富集分析，發(fā)現(xiàn)差異表達(dá)的基因主要富集在肌肉發(fā)育的生物學(xué)過(guò)程和脂肪酸代謝等通路[9]。但斷奶仔豬背最長(zhǎng)肌和皮下脂肪組織發(fā)育相關(guān)的基因及其功能的研究仍處于初級(jí)階段。

為了深入研究在仔豬背最長(zhǎng)肌和皮下脂肪組織生長(zhǎng)發(fā)育過(guò)程起關(guān)鍵作用的基因，本研究采用數(shù)據(jù)庫(kù)數(shù)據(jù)進(jìn)行生物信息學(xué)分析。本研究詳細(xì)介紹了可能影響仔豬生長(zhǎng)性能的基因及通路，為仔豬生長(zhǎng)發(fā)育相關(guān)分子生物學(xué)的研究提供理論支撐。在未來(lái)研究中，將進(jìn)一步結(jié)合分子生物學(xué)和細(xì)胞生物學(xué)實(shí)驗(yàn)，揭示關(guān)鍵通路中基因的功能以及發(fā)揮功能的分子機(jī)制，為減少仔豬斷奶應(yīng)激、改善生長(zhǎng)發(fā)育提供理論依據(jù)，對(duì)生豬產(chǎn)業(yè)具有重要意義。