基于轉(zhuǎn)錄組測(cè)序篩選黃牛低氧適應(yīng)性相關(guān)差異基因

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

(1.青藏高原動(dòng)物遺傳資源保護(hù)與利用教育部重點(diǎn)實(shí)驗(yàn)室,西南民族大學(xué),四川 成都 610041； 2.青藏高原動(dòng)物遺傳資源保護(hù)與利用四川省重點(diǎn)實(shí)驗(yàn)室,西南民族大學(xué),四川 成都 610041)

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

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

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

1 材料和方法

1.1 樣品采集

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

1.2 總RNA提取、cDNA文庫(kù)構(gòu)建及Illumina測(cè)序

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

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

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

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

使用 Cufflink(http://cufflinks.cbcb.umd.edu/ howitworks.html)軟件根據(jù)比對(duì)結(jié)果來重構(gòu)轉(zhuǎn)錄本,對(duì)比重構(gòu)轉(zhuǎn)錄本與參考基因組(ARS-UCD1.2)序列(篩選轉(zhuǎn)錄本長(zhǎng)度≥200 bp且外顯子數(shù)目≥2),如果新轉(zhuǎn)錄本與參考基因組對(duì)比上,則得到樣本已知的轉(zhuǎn)錄本和新的轉(zhuǎn)錄本。重構(gòu)轉(zhuǎn)錄本可能由于延長(zhǎng)了轉(zhuǎn)錄本注釋的5′或者3′端,因此優(yōu)化了轉(zhuǎn)錄本結(jié)構(gòu)[18]。 利用 rMATS 軟件(http://rnaseq-mats.sourceforge.net/index.html)進(jìn)行可變剪切事件的分析。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 差異表達(dá)mRNA的 GO富集與KEGG信號(hào)通路分析

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

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

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

1.5 RT-qPCR檢測(cè)

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

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

2 結(jié)果與分析

2.1 測(cè)序數(shù)據(jù)質(zhì)量評(píng)估及分析

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

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

注:HAC1、HAC2、HAC3分別是藏黃牛的3個(gè)生物學(xué)重復(fù)；LAC1、LAC2、LAC3是三江黃牛的3個(gè)生物學(xué)重復(fù)。表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 轉(zhuǎn)錄組測(cè)序數(shù)據(jù)對(duì)比分析

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

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

2.3 新轉(zhuǎn)錄本預(yù)測(cè)

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

表4 轉(zhuǎn)錄本統(tǒng)計(jì)匯總Tab.4 Transcript statistics summary

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

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

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

對(duì)所有樣本轉(zhuǎn)錄本的表達(dá)量豐度進(jìn)行比較,其中橫坐標(biāo)為log10(FPKM)，數(shù)值越高表示轉(zhuǎn)錄本表達(dá)量越高；縱坐標(biāo)為轉(zhuǎn)錄本的豐度，即為對(duì)應(yīng)橫軸表達(dá)量的轉(zhuǎn)錄本數(shù) / 檢測(cè)已表達(dá)轉(zhuǎn)錄本的總數(shù)。結(jié)果如圖2顯示,三江黃牛3個(gè)樣本的轉(zhuǎn)錄本表達(dá)量豐度明顯高于藏黃牛樣本轉(zhuǎn)錄本的表達(dá)量。

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

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

圖3 樣品間表達(dá)量相關(guān)性熱圖Fig.3 Correlation heat map of gene expression level in samples

2.5 mRNA差異表達(dá)分析

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

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

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

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 差異表達(dá)基因火山圖
Fig.4 Volcano plot of differentially expressed genes

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

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

2.6 差異表達(dá)基因的GO功能富集分析

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

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

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

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

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

圖8 cGMP-PKG信號(hào)通路Fig.8 cGMP-PKG signaling pathway

2.8 RT-qPCR驗(yàn)證

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

圖9 轉(zhuǎn)錄組測(cè)序結(jié)果與熒光定量結(jié)果比較分析Fig.9 Comparison analysis of the RNA-Seq and RT-PCR result

3 討論與結(jié)論

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

3.1 轉(zhuǎn)錄組測(cè)序數(shù)據(jù)分析

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

3.2 差異表達(dá)mRNA分析

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

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

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

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

3.3 GO功能顯著富集分析

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

3.4 KEGG富集通路分析

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

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

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