基于轉錄組測序篩選黃牛低氧適應性相關差異基因

侯孟典,王 會,鐘金城,柴志欣,益西康珠,王吉坤,王嘉博

(1.青藏高原動物遺傳資源保護與利用教育部重點實驗室,西南民族大學,四川 成都 610041； 2.青藏高原動物遺傳資源保護與利用四川省重點實驗室,西南民族大學,四川 成都 610041)

青藏高原作為世界海拔最高的高原,平均海拔為4 000～5 000 m,因此造成了高寒、低氧、紫外線極強等惡劣多變的高原氣候。寒冷和氧氣含量低等原因給恒溫脊椎動物的生存帶來了嚴峻的挑戰(zhàn),由于氧氣分壓的降低減緩了呼吸組織細胞對氧氣的供應,因此限制了能量的代謝而抑制了熱量的產(chǎn)生,研究高原動物如何抵抗低氧與寒冷環(huán)境有助于了解低氧適應性進化進程[1]。藏黃牛是青藏高原地區(qū)的傳統(tǒng)飼養(yǎng)品種,主要分布在海拔3 000～4 500 m的農(nóng)區(qū)、林區(qū)及半農(nóng)半牧區(qū),以放牧飼養(yǎng)為主,可肉、乳、役兼用,是經(jīng)長期自然選擇而逐漸適應高原氣候的小型地方原始品種[2]。由于高海拔地區(qū)獨特的環(huán)境特點,藏黃牛進化出與低氧環(huán)境相適應的生理特征,例如獨特的心臟結構、單位面積內(nèi)肺泡數(shù)增多等[3-9]。藏黃牛作為高原地區(qū)所特有的原始品種,既是促進當?shù)亟?jīng)濟發(fā)展的重要動物,也是作為研究黃牛低氧適應性的理想物種。三江黃牛原產(chǎn)于四川省阿壩藏族自治州汶川縣,具有良好的役用性能、軀干較長、肉質(zhì)好、適應性與抗病性強等特點,是遺傳資源中一個極為珍貴的基因庫[10]。

據(jù)報道,長期生活在低溫、低氧及高輻射的高海拔地區(qū)的動物進化出與低海拔地區(qū)動物不同的呼吸系統(tǒng),例如支氣管管壁變厚、肺泡內(nèi)毛細血管豐富且呈擴張狀態(tài)、心血管系統(tǒng)發(fā)達等[11-12]。齊曉園等[13]對牦牛與藏黃牛血液生理生化指標的研究發(fā)現(xiàn),隨著海拔高度的升高,牦牛通過使血液中紅細胞和血紅蛋白含量不斷增加、增大血細胞比容等方式來適應高原低氧環(huán)境,黃牛通過增加紅細胞的體積、增加血液中血紅蛋白含量與血紅蛋白濃度等方式來適應低氧環(huán)境,從而發(fā)現(xiàn)牦牛與黃牛對于應對低氧適應性產(chǎn)生了不同的機制；王志敏等[14]對藏雞低氧適應性機制研究發(fā)現(xiàn),隨著海拔高度的升高,紅細胞數(shù)、血紅蛋白含量及平均紅細胞體積隨之升高,而平均血紅蛋白量與紅細胞平均血紅蛋白濃度隨著海拔的升高而降低,并且發(fā)現(xiàn)處在同一海拔高度的公母雞之間的血液生理生化指標也存在差異,公雞的各項指標均高于母雞。

心臟是脊椎動物身體中對氧含量變化較敏感的重要器官之一,為血液流動至全身各部位提供動力,使細胞能夠正常的代謝及運作,研究發(fā)現(xiàn)藏黃牛進化出與三江黃牛不同的心血管系統(tǒng),且在低氧情況下能夠通過增加呼吸頻率,加強心泵功能等方式進行自身調(diào)控來適應低氧環(huán)境[15-16],因此，本研究擬通過對藏黃牛和三江黃牛心臟組織進行轉錄組測序分析,以期獲得高原低氧適應性相關的基因及其涉及的信號通路,比較不同海拔高度黃牛心臟組織中差異表達轉錄本,為高海拔環(huán)境脅迫下產(chǎn)生的適應性進化機制研究奠定基礎。

1 材料和方法

1.1 樣品采集

分別選取健康的4.5歲、位于西藏自治區(qū)昌都市類烏齊縣的藏黃牛(海拔高度:3 810 m)與四川省成都市郫都區(qū)的三江黃牛(海拔高度:500 m)各3頭,采集其心臟組織,經(jīng)DEPC水沖洗后,用錫箔紙包裹置于液氮中帶回實驗室保存?zhèn)溆谩?/p>

1.2 總RNA提取、cDNA文庫構建及Illumina測序

取出液氮中放置的各黃牛心臟組織,加入液氮于研缽中研磨,研磨至粉末狀后移入1.5 mL無RNase的離心管,加入1 mL TRIzol試劑(購自Invitrogen公司,美國)提取總RNA,提取步驟參照 TRIzol試劑說明書。提取的總RNA經(jīng)紫外分光光度計測試其濃度與OD260nm/280nm值,再利用1%的瓊脂糖凝膠電泳(BIORAD)檢測其質(zhì)量。將提取的總RNA反轉錄構建cDNA文庫,建好的測序文庫用Illumina HiSeqTM 4000平臺進行測序。

1.3 數(shù)據(jù)處理及分析

為了保證數(shù)據(jù)質(zhì)量,在生物信息學分析前對數(shù)據(jù)進行質(zhì)量控制與過濾,步驟如下:先將初步過濾得到的干凈數(shù)據(jù)序列(Clean reads)進行進一步更嚴格的過濾,去除含有接頭的reads、全部都是A堿基的reads、含N比例大于10%的reads和低質(zhì)量的reads(質(zhì)量值Q≤20的堿基數(shù)占整條reads的50%以上),最終得到高質(zhì)量數(shù)據(jù)用于后續(xù)的信息分析。

使用比對軟件Tophat2(2.1.1)將數(shù)據(jù)對比到黃牛參考基因組(ARS-UCD1.2)進行分析。利用FPKM(Fragments Per Kilobase of transcript per Million mapped reads)法計算轉錄本表達量并進行統(tǒng)計。根據(jù)FPKM定量結果,使用R語言包(http://www.r-project.org/)進行熱圖繪制,計算出所有樣品兩兩之間的相關性。樣品間基因表達水平相關性是檢驗試驗可靠性和樣本選擇是否合理的重要指標,相關系數(shù)越接近1,表明樣品表達模式的相似度越高[13,17]。

使用 Cufflink(http://cufflinks.cbcb.umd.edu/ howitworks.html)軟件根據(jù)比對結果來重構轉錄本,對比重構轉錄本與參考基因組(ARS-UCD1.2)序列(篩選轉錄本長度≥200 bp且外顯子數(shù)目≥2),如果新轉錄本與參考基因組對比上,則得到樣本已知的轉錄本和新的轉錄本。重構轉錄本可能由于延長了轉錄本注釋的5′或者3′端,因此優(yōu)化了轉錄本結構[18]。 利用 rMATS 軟件(http://rnaseq-mats.sourceforge.net/index.html)進行可變剪切事件的分析。rMATS軟件可以將可變剪切事件分成5類:外顯子跳躍(Skipped exon,SE)、可變5′端剪切(Alternative 5′ splice site,A5SS)、可變3′端剪切(Alternative 3′ splice site,A3SS)、外顯子選擇性跳躍(Mutually exclusive exon,MXE)、內(nèi)含子保留(Retained intron,RI)[19-20]。

1.4 差異表達mRNA的 GO富集與KEGG信號通路分析

使用edgeR軟件以RPKM(Reads per kilobase transcriptome per million mapped reads)法計算表達量,對藏黃牛與三江黃牛心臟組織mRNA表達量進行差異分析,FDR(False discovery rate)是多重檢驗校正P值的一種方法,以篩選條件為:FDR<0.05 且 |log2FC|>1(FC為3個樣品表達量的平均值)來篩選差異轉錄本。同時利用皮爾森算法計算樣品間的距離,使用皮爾遜相關系數(shù)r(Pearson′s correlation coefficient)來表示生物學重復相關性,r2越接近1就表示2個樣品間相關性越強,對同一條件下每一對生物學重復樣品的基因表達量的相關性作圖,根據(jù)距離大小建立聚類圖。

對得到的差異轉錄本進行GO功能分析將差異表達轉錄本向GO數(shù)據(jù)庫(http://www.geneontology.org/)的各term(詞條、節(jié)點)映射,并計算每個條目的轉錄本數(shù)[21]。GO功能分析主要包括分子功能(Molecular function)、細胞組分(Cellular component)、生物過程(Biological process)。應用超幾何檢驗,找出與整個轉錄本組背景相比,在差異表達轉錄本中顯著富集的GO條目。P值通過FDR校正之后得到Q值,以Q≤0.05為閾值,滿足此條件的GO條目定義為在差異表達轉錄本中顯著富集的GO條目。

通過KEGG(Kyoto encyclopedia of genes and genomes)數(shù)據(jù)庫進行比對[22],應用超幾何檢驗,找出差異表達轉錄本中顯著性富集的通路,確定差異表達轉錄本參與的主要生化代謝途徑和信號轉導途徑[23]。

1.5 RT-qPCR檢測

為驗證測序結果是否準確,在差異表達基因中隨機選取6個基因,以GADPH為內(nèi)參基因,用SYBR?Premix Ex TaqTMⅡ試劑,以反轉錄的cDNA為模板利用CFX96 實時定量PCR 系統(tǒng)(美國Bio-Rad 公司)進行實時熒光定量PCR(Real-time PCR)驗證。使用Primer Premier 5.0軟件設計熒光定量引物(表1)并由賽默飛世爾(上海)科技有限公司合成。10 μL反應體系:2×SYBR 5 μL；上下游引物各0.4 μL；cDNA模板0.8 μL；無RNAase水3.4 μL。定量反應程序:95 ℃預變性30 s；95 ℃變性5 s,60 ℃退火30 s,39個循環(huán),同時記錄溶解度曲線用來鑒定引物的特異性,同一樣品設置3個技術重復。利用2-ΔΔCT法計算各個樣本的相對表達量,結果以“平均值±標準差”來表示。

表1 熒光定量PCR引物Tab.1 RT-qPCR primers

2 結果與分析

2.1 測序數(shù)據(jù)質(zhì)量評估及分析

本研究共構建藏黃牛與三江黃牛心臟組織轉錄組文庫6個,測序得到的原始數(shù)據(jù)經(jīng)過濾,高質(zhì)量數(shù)據(jù)在測序數(shù)據(jù)中所占比例均高于98%(表2)。經(jīng)過過濾后樣品GC含量均不小于50%,從堿基組成與質(zhì)量分析圖可知,過濾后樣品的堿基分布穩(wěn)定,Q30(數(shù)據(jù)質(zhì)量達到Q30即測序準確率為99.9%的堿基數(shù)目及百分比)值均大于95.9%,證明測序質(zhì)量可靠,可以用于后續(xù)試驗分析[24]。

表2 經(jīng)過過濾數(shù)據(jù)統(tǒng)計Tab.2 Statistics of filtered data

注:HAC1、HAC2、HAC3分別是藏黃牛的3個生物學重復；LAC1、LAC2、LAC3是三江黃牛的3個生物學重復。表3-4同。

Note: HAC1, HAC2, and HAC3 are three biological replicates of Tibetan cattle；LAC1, LAC2, and LAC3 are three biological replicates of Sanjiang cattle. The same as Tab.3-4.

2.2 轉錄組測序數(shù)據(jù)對比分析

比對黃牛參考基因組(ARS-UCD1.2)結果表明,唯一比對率均在86%以上,對比率均高于87%(表3),表明轉錄組測序數(shù)據(jù)質(zhì)量穩(wěn)定[25-26]。結果可知唯一能對比到的reads數(shù)在3個藏黃牛轉錄組中平均有81 930 524條,3個三江黃牛轉錄組中平均有79 135 363條,藏黃牛轉錄本數(shù)高于三江黃牛,而三江黃牛轉錄組對比率(Mapping ratio)高于藏黃牛。

表3 去除核糖體RNA數(shù)據(jù)對比到參考基因組Tab.3 Removal rRNA data to the reference genome

2.3 新轉錄本預測

使用Cufflink軟件根據(jù)比對結果來重構轉錄本,將重構結果與黃牛參考轉錄本(ARS-UCD1.2)序列進行比較,得到已知的轉錄本和新的轉錄本(表4),并對新轉錄本的基因組信息進行功能注釋。結果可知,3個藏黃牛轉錄本中平均預測出新mRNA 5 238個,三江黃牛轉錄本中預測出新mRNA 3 331個。

表4 轉錄本統(tǒng)計匯總Tab.4 Transcript statistics summary

藏黃牛與三江黃牛樣品可變剪切結果可知(圖1),共得到差異可變剪切事件4 742個。在5類可變剪切事件模型中以SE模型為主,占所有差異可變剪切事件的62.5%；RI模型在差異可變剪切中所占比例最小,為3.5%。

圖1 差異可變剪切分類和數(shù)量統(tǒng)計圖Fig.1 Alternative splicing and quantity statistics

2.4 不同海拔高度黃牛心臟mRNA表達水平比較

對所有樣本轉錄本的表達量豐度進行比較,其中橫坐標為log10(FPKM)，數(shù)值越高表示轉錄本表達量越高；縱坐標為轉錄本的豐度，即為對應橫軸表達量的轉錄本數(shù) / 檢測已表達轉錄本的總數(shù)。結果如圖2顯示,三江黃牛3個樣本的轉錄本表達量豐度明顯高于藏黃牛樣本轉錄本的表達量。

圖2 所有樣本表達量豐度分布圖Fig.2 FPKM distribution of all samples

根據(jù)FPKM計算結果對轉錄本表達量進行計算,根據(jù)結果計算出各樣品之間的相關性(圖3)。以三江黃牛樣本1與2組比較為例,r2值為0.973 7,表明其重復性高,數(shù)據(jù)準確且穩(wěn)定。藏黃牛樣本的相關性系數(shù)r2平均值為0.889 0,三江黃牛樣本的相關性系數(shù)r2平均值為0.972 6。

圖3 樣品間表達量相關性熱圖Fig.3 Correlation heat map of gene expression level in samples

2.5 mRNA差異表達分析

以FDR 與 log2FC 來篩選差異轉錄本結果表明,藏黃牛心臟組織轉錄組文庫相比于三江黃牛心臟組織轉錄組文庫分析,共發(fā)現(xiàn)顯著差異表達基因608個,其中上調(diào)基因467個,下調(diào)基因141個(圖4)。在藏黃牛與三江黃牛中表達量最高且表達差異倍數(shù)較大的基因分別是β2微球蛋白基因(Beta-2-microglobulin,B2M)與細胞色素C氧化酶7C亞基基因(Cytochrome c oxidase subunit 7C,COX7C)。另外,活化T細胞核因子1基因(Nuclear factor of activated T cells 1,NFATC1)、內(nèi)皮素轉化酶1基因(Endothelin converting enzyme 1,ECE1)、鈣蛋白酶抑制蛋白基因(Calpastatin,CAST)在藏黃牛與三江黃牛中表達量也較高,差異倍數(shù)較大,且FDR值較低。

對篩選出的差異表達本進行聚類分析,每列代表一個樣品，每行代表一個轉錄本，顏色越紅表示表達量越高，顏色越藍表示表達量越低。結果如圖所示(圖5)。在縱向來看,藏黃牛(HAC)與三江黃

圖中的每個點代表一個基因,紅色和綠色的點分別表示上調(diào)與下調(diào)的顯著差異表達基因,黑色的點表示沒有差異。

Each point in the figure represents a gene, and the red and green dots represent significant up-regulated expressed genes and down-regulated genes, the black dots indicate no difference.

圖4 差異表達基因火山圖
Fig.4 Volcano plot of differentially expressed genes

圖5 差異表達基因聚類分析圖熱圖Fig.5 Differential expression gene clustering analysis heat map

牛(LAC)的3個生物學重復樣本聚類在一起,說明生物學重復樣品間具有很高的相關性。從縱向來看顯著性差異基因可以分為上下2個基因簇,結果顯示上調(diào)與下調(diào)基因的聚類是成功的,與篩選出的差異表達基因結果相符合,并且上調(diào)基因之間及下調(diào)基因間相關性均很高[27]。

2.6 差異表達基因的GO功能富集分析

608個顯著差異表達mRNA經(jīng)過功能富集分析在分子功能、細胞組分、生物學過程分別富集到327,253,2 191個條目。結果顯示,在生物學過程本體中,細胞發(fā)育過程條目富集程度最高；在分子功能本體中,結合與轉錄因子活性條目富集程度最高；在細胞組分本體中,細胞前緣條目富集程度最高。由圖6可知,以上下調(diào)基因進行分類可以富集在50個條目(Term),其中參與分子功能、細胞組分、生物過程的條目分別為11,18,21個。參與分子功能的條目中,涉及參與基因數(shù)目最多的3個條目分別為：細胞過程(Cellular process)條目,其上調(diào)基因274個，下調(diào)基因80個；單一生物過程(Single-organism process)條目其上調(diào)基因239個，下調(diào)基因72個；代謝過程(Metabolic process)條目，其上調(diào)基因190個，下調(diào)基因56個。參與細胞組分的條目中,涉及參與基因最多的條目分別為:核酸結合轉錄因子活性(Nucleic acid binding transcription factor activity條目,其上調(diào)基因27個，下調(diào)基因7個；結合(Binding)條目,上調(diào)基因256個，下調(diào)基因71個；催化活性(Catalytic activity)條目,其上調(diào)基因103個，下調(diào)基因有34個。參與生物過程的條目中細胞器(Organelle)、細胞(Cell)、細胞部分(Cell part)條目參與基因數(shù)目最多,其參與基因數(shù)目分別為上調(diào)基因212個，下調(diào)基因65個，上調(diào)基因258個，下調(diào)基因77個，上調(diào)基因258個，下調(diào)基因77個。

圖6 差異基因的GO富集分析柱狀圖Fig.6 GO enrichment analysis histogram of differential genes

2.7 差異mRNA的KEGG通路富集分析

經(jīng)過KEGG通路富集分析,差異表達轉錄本中顯著富集通路共有8條,分別是cGMP-PKG信號通路、T細胞受體信號通路、FcγR介導的吞噬作用、甲狀腺激素信號通路、軸突指導、磷脂酶D信號通路、白細胞跨內(nèi)皮遷移、肌動蛋白細胞骨架調(diào)節(jié)通路(圖7)。其中差異表達基因NFATC1顯著富集在cGMP-PKG信號通路與T細胞受體信號通路中(圖8)。

圖7 KEGG富集氣泡圖Fig.7 KEGG enriched bubble chart

圖8 cGMP-PKG信號通路Fig.8 cGMP-PKG signaling pathway

2.8 RT-qPCR驗證

隨機選取6個差異表達基因,分別是CAST、TANGO2、GUK1、TAF6、SQRDL、RPS14做RT-qPCR驗證,根據(jù)結果(圖9)可知,轉錄組測序數(shù)據(jù)與熒光定量結果基因的表達水平變化趨勢一致,判定轉錄組測序數(shù)據(jù)可靠。

圖9 轉錄組測序結果與熒光定量結果比較分析Fig.9 Comparison analysis of the RNA-Seq and RT-PCR result

3 討論與結論

近年來隨著高通量測序技術地迅速發(fā)展,轉錄組測序技術被廣泛地應用在基因的表達及調(diào)控研究中[27-32]。通過轉錄組技術對高原低氧環(huán)境中動物適應性的研究已經(jīng)在藏豬、藏雞[14]、田鼠和牦牛[33]等動物中有所研究。本試驗首次采用轉錄組測序技術對藏黃牛與三江黃牛心臟組織中差異表達轉錄本進行研究,為研究藏黃牛低氧適應性機制提供理論基礎。

3.1 轉錄組測序數(shù)據(jù)分析

本研究使用Illumina HiSeqTM 4000測序平臺構建了藏黃牛與三江黃牛心臟組織轉錄組本庫6個,經(jīng)過數(shù)據(jù)處理后獲得了大量的高質(zhì)量數(shù)據(jù)。對比到參考基因組后藏黃牛獲得唯一比對的reads比例大于86%,且對比率均大于87%；三江黃牛唯一比對的reads比例大于91.5%,對比率均大于92%,三江黃牛與參考基因組對比率高于藏黃牛。通過對mRNA表達水平的豐度研究,發(fā)現(xiàn)三江黃牛mRNA表達量要高于藏黃牛,各個樣品間相關性熱圖可以看出各個樣品間的相關性較好,同樣三江黃牛樣品的相關性高于藏黃牛,重復性散點圖得到的結果與其一致,各結果均表明測序數(shù)據(jù)穩(wěn)定且質(zhì)量較高。

3.2 差異表達mRNA分析

本研究對藏黃牛與三江黃牛心臟組織mRNA差異表達分析,發(fā)現(xiàn)在藏黃牛與三江黃牛中表達量最高,且表達差異倍數(shù)較大的B2M、COX7C、NFATC1、ECE1、CAST基因可能是與低氧適應性相關的基因。研究發(fā)現(xiàn)缺氧環(huán)境會導致肌纖維生長緩慢[34],CAST基因是與肉質(zhì)性狀相關的基因,推測藏黃?？赡芡ㄟ^增大心肌面積,增加心肌的密度,減小心肌纖維直徑等組織學改變來適應高原低氧環(huán)境。

細胞色素C氧化酶 (Cytochrome c oxidase,Cox)是線粒體中呼吸鏈電子傳遞的終末復合物,在能量的產(chǎn)生與代謝中起重要作用[35],COX7C基因是發(fā)揮能量代謝作用中一個起調(diào)節(jié)作用的重要亞基。研究發(fā)現(xiàn)其上有與線粒體基因表達轉錄因子作用的位點,許多魚類在低溫環(huán)境時通過增加細胞色素C氧化酶活性來抵御寒冷[36],該基因功能與線粒體活動密切相關。本研究中COX7C基因在藏黃牛與三江黃牛中顯著差異,而且在三江黃牛中下調(diào)并且表達量最高,這與藏黃牛通過增強能量代謝來抵御寒冷的特征相符[34],間接說明了藏黃牛對線粒體的調(diào)度與利用率高于三江黃牛。

活化T細胞核因子1基因(Nuclear fator of activated T cell 1,NFATC1)是Rel轉錄因子家族中的一員,其通過與Ca2+結合的途徑進入靶細胞核內(nèi)從而影響DNA的復制與轉錄,誘導心內(nèi)組織的分化,從而形成心內(nèi)膜和主肺動脈瓣等。研究證明NFATC1與心臟的房室分隔形成有關[37],進一步研究證明NFATC1可以通過降低血管內(nèi)皮生長因子表達而使心內(nèi)膜細胞分化,從而最終形成心臟瓣膜[38]。本研究中NFATC1基因在藏黃牛與三江黃牛中差異顯著且表達量較高,可以推測其是與低氧適應性相關的基因,并且它還參與cGMP-PKG信號通路與T細胞受體信號通路,進一步證明其可能是與低氧適應性相關的基因。

內(nèi)皮素轉化酶 1(Endothelin-converting enzyme-1,ECE)是生物體內(nèi)合成過程中重要的限速酶,其基因CEC1是與低氧適應性相關的重要候選基因之一,ECE主要位于心血管系統(tǒng)內(nèi),屬于中性肽鍵內(nèi)切酶(NEP)家族,它有3種異構體并且廣泛分布于組織內(nèi),研究表明,體內(nèi)在低氧誘導情況下可以通過上調(diào)ECE1的表達從而促進內(nèi)皮素-1的合成,這與本研究中ECE1基因在藏黃牛與三江黃牛中差異表達且在藏黃牛中上調(diào)的結果相符[39]。王延東[40]對ECE1基因在藏豬、萊蕪豬、大白豬的研究中發(fā)現(xiàn),ECE1基因在藏豬中的表達明顯高于其他2個平原豬,在藏豬10個組織中均有表達并且在心臟組織與肺臟組織的表達中最高。心肺組織是機體對抗低氧環(huán)境的重要器官,本研究中ECE1基因在藏黃牛心臟轉錄組中顯著表達且上調(diào),更進一步證明了ECE1基因可能是黃牛與低氧適應性相關的候選基因。

3.3 GO功能顯著富集分析

據(jù)報道,在生物體內(nèi)氧氣的運輸主要由紅細胞與血紅蛋白負責,同樣在細胞內(nèi)低氧環(huán)境時可以誘導線粒體發(fā)生自噬以適應低氧環(huán)境[41],牦牛、藏豬、藏雞等高海拔動物,在低氧環(huán)境時通過不斷增加血液中紅細胞和血紅蛋白含量[2]、增大血紅細胞體積[30]等生理生化指標的改變來適應低氧環(huán)境。本研究發(fā)現(xiàn),GO功能富集分析結果共涉及3大類別的50個功能條目,在生物學過程、分子功能、細胞組分3個本體中富集程度最高的條目分別是細胞發(fā)育過程、結合與轉錄因子活性和細胞前緣條目,說明富集程度高的條目都是與細胞過程有關的功能條目；差異表達mRNA功能主要富集在與細胞各項功能相關的條目上,因此,推測藏黃牛血紅蛋白在長期低氧環(huán)境中有獨特的攜帶氧氣的能力來適應低氧環(huán)境。

3.4 KEGG富集通路分析

經(jīng)過多重效驗校正后的KEGG分析顯著富集的通路有8條,其中最顯著富集的通路是cGMP-PKG信號通路(cGMP-PKG signaling pathway)是與心血管疾病相關的通路。cGMP:環(huán)磷酸鳥苷酸(cyclic guanosine 3′,5′ monophosphate)為一種第二信使,由鳥苷酸環(huán)化酶作用于GTP(三磷酸鳥苷)而形成,可以在多種組織細胞中發(fā)揮生物功能；PKG:環(huán)磷酸鳥苷酸依賴的蛋白激酶,是一種真核細胞中廣泛存在的絲/蘇氨酸蛋白激酶[42]。在高海拔地區(qū)由于氧氣含量低,容易引起心臟衰竭等臨床癥狀。研究表明,在心臟發(fā)生衰竭時,可以通過β3-AR(腎上腺素能受體)作用于cGMP-PKG-P38信號通路使心肌細胞產(chǎn)生巨噬細胞移動抑制因子(Macrophage migrationinhibitory factor)而避免心肌細胞凋亡[43]。心肌在缺血時可以通過再灌注治療使其可以較快的恢復心臟血液供應,但是不可避免的對心臟產(chǎn)生損傷。有研究表明,NO-cGMP-PKG通路可能與保護作用相關,其可以通過腦內(nèi)神經(jīng)產(chǎn)生特異性結合物質(zhì)激活cGMP使心肌內(nèi)濃度升高進而激活PKG,此活動是通過NO-cGMP-PKG通路實現(xiàn)的[44]。

KEGG富集通路中,肌動蛋白細胞骨架調(diào)節(jié)通路(Regulation of actin cytoskeleton)顯著富集。研究表明,它們是與心血管發(fā)育、肌肉組織結構、肌肉細胞分化等相關的通路,肌動蛋白骨架在維持細胞形態(tài)、調(diào)節(jié)細胞粘連等方面具有重要功能,而粘著斑結則參與細胞膜受體與肌動蛋白骨架之間的結構連接[21]。趙啟南等[45]在研究蒙古馬高負荷運動后肌肉組織轉錄組差異分析發(fā)現(xiàn),差異表達基因被富集于擴張型心肌病、心肌收縮、鈣信號通路、肌動白骨架調(diào)節(jié)等與運動相關的多條通路中。此結果與本研究相似,高強度運動時由于肺部氣體交換不及時,會在體內(nèi)形成缺氧癥狀,對機體各組織器官造成損傷,高海拔地區(qū)由于環(huán)境惡劣,氧分壓低造成氧氣含量少,高海拔地區(qū)動物逐漸形成了獨特的低氧適應機制,因此認為肌動蛋白細胞骨架調(diào)節(jié)通路可以作為與低氧適應性相關的通路。

本試驗通過高通量測序技術對藏黃牛與三江黃牛心臟組織轉錄組測序分析,獲得轉錄組測序數(shù)據(jù)穩(wěn)定且質(zhì)量良好；經(jīng)過篩選共得到顯著差異表達基因608個,其中上調(diào)基因467,下調(diào)基因141個,其中表達量較高且差異倍數(shù)較大的COX7C、ECE1、B2M、NFATG、CAST等基因可能是與低氧適應性相關的基因；經(jīng)過GO富集與KEGG Pathway信號通路分析,GO富集結果可知細胞發(fā)育過程(Cellular developmental process)、結合(Binding)、細胞前緣(Cell leading edge)等條目顯著富集；KEGG富集分析得到多條相關通路,多重效驗后顯著富集的通路共有8條,其中cGMP-PKG信號通路與肌動蛋白細胞骨架調(diào)節(jié)通路可能是與低氧適應性相關的通路。本研究為探究動物低氧適應性機制提供理論基礎。