熱應(yīng)激對(duì)中國荷斯坦牛乳腺組織基因表達(dá)及信號(hào)通路的影響

李秋玲，齊 穎，王 琛，張一名，王新妤，尚校蘭，賈永紅，李美茹，儲(chǔ)明星

(1.廊坊師范學(xué)院 生命科學(xué)學(xué)院，河北省動(dòng)物多樣性重點(diǎn)實(shí)驗(yàn)室，廊坊市細(xì)胞工程與應(yīng)用研究重點(diǎn)實(shí)驗(yàn)室，河北 廊坊 065000； 2.中國農(nóng)業(yè)科學(xué)院 北京畜牧獸醫(yī)研究所，農(nóng)業(yè)農(nóng)村部動(dòng)物遺傳育種與繁殖(家禽)重點(diǎn)實(shí)驗(yàn)室，北京 100193)

熱應(yīng)激是動(dòng)物機(jī)體受到環(huán)境高溫刺激所產(chǎn)生的一系列非特異性應(yīng)答反應(yīng)的總和[1]。隨著溫室效應(yīng)加劇，夏秋季高溫氣候出現(xiàn)頻率逐漸升高，奶牛熱應(yīng)激問題日益凸顯，嚴(yán)重影響奶業(yè)生產(chǎn)，成為奶牛夏季生產(chǎn)的一個(gè)重大難題。中國荷斯坦牛原產(chǎn)于歐洲，耐寒性較強(qiáng)，耐熱性較差，夏季容易發(fā)生熱應(yīng)激，輕者引起新陳代謝紊亂[2]，生產(chǎn)性能、乳品質(zhì)[3]、繁殖力、免疫力下降[4]，生產(chǎn)年限縮短，重者機(jī)體功能衰竭導(dǎo)致死亡，給奶業(yè)生產(chǎn)帶來巨大的經(jīng)濟(jì)損失。據(jù)統(tǒng)計(jì)，美國奶業(yè)每年因熱應(yīng)激造成的損失高達(dá)9億美元[5]。

近年來，我國奶牛單產(chǎn)不斷提高，但產(chǎn)奶量與耐熱程度高度負(fù)相關(guān)，另外，選種時(shí)主要考慮生產(chǎn)性能使奶牛的耐熱性不斷下降。隨著動(dòng)物營養(yǎng)學(xué)的快速發(fā)展，阻礙生產(chǎn)性能提高的因素已不再是日糧營養(yǎng)供應(yīng)問題，更多的是奶牛飼養(yǎng)環(huán)境的限制。熱應(yīng)激作為一個(gè)重要的環(huán)境因素已引起廣大學(xué)者的高度關(guān)注，熱應(yīng)激發(fā)生時(shí)動(dòng)物體內(nèi)各種基因調(diào)控的關(guān)系有待深入研究，從而揭示熱應(yīng)激的作用機(jī)制和規(guī)律。本研究擬通過Illumina測(cè)序平臺(tái)對(duì)奶牛熱應(yīng)激期和非熱應(yīng)激期乳腺組織轉(zhuǎn)錄組文庫進(jìn)行高通量測(cè)序和生物信息學(xué)分析，研究熱應(yīng)激奶牛乳腺組織差異表達(dá)基因及其調(diào)控通路，為進(jìn)一步研究奶牛發(fā)生熱應(yīng)激的分子機(jī)制，耐熱奶牛分子育種提供理論依據(jù)。

1 材料與方法

1.1 試驗(yàn)材料

根據(jù)系譜資料、奶牛生產(chǎn)性能(DHI)測(cè)定數(shù)據(jù)和牛場(chǎng)記錄數(shù)據(jù)選取性別相同，飼養(yǎng)環(huán)境、年齡及胎次相同的8頭中國荷斯坦牛作為實(shí)驗(yàn)牛。根據(jù)本實(shí)驗(yàn)室前期研究對(duì)奶牛熱應(yīng)激和非熱應(yīng)激期生理參數(shù)比較結(jié)果[6]，采集樣本，即分別在3月(溫度15～20℃，非熱應(yīng)激期，non-heat stressed，NHS)和8月(溫度30～38℃，熱應(yīng)激期，heat stressed，HS)屠宰實(shí)驗(yàn)牛采集乳腺組織樣本。NHS對(duì)照組編號(hào)為A、B、C和D，HS組樣本編號(hào)為E、F、G和H，保存于-80 ℃，用于提取組織總RNA。

1.2 總RNA提取和鑒定

提取熱應(yīng)激和非熱應(yīng)激期奶牛乳腺組織總RNA，利用核酸測(cè)量儀檢測(cè)濃度及純度，變性聚丙烯酰胺凝膠電泳鑒定所提總RNA質(zhì)量，用Agilent 2100 BioAnalyzer生物分析儀進(jìn)行質(zhì)檢，分析RNA完整性(RNA integrity number，RIN)。

1.3 文庫制備和測(cè)序

用Poly (T)寡聚核苷酸提取Poly (A)+RNA，進(jìn)行片段化處理。去核糖體RNA后，純化剩余RNA，片段化和用于cDNA合成。反轉(zhuǎn)錄合成cDNA第一鏈，反轉(zhuǎn)錄合成cDNA第二鏈。對(duì)cDNA進(jìn)行末端修復(fù)并添加測(cè)序接頭 (adapter)。對(duì)連接產(chǎn)物切膠純化再進(jìn)行PCR擴(kuò)增(大小200～300 bp)，獲得cDNA文庫，利用HiSeq2000進(jìn)行高通量測(cè)序。

1.4 原始序列的處理

通過去接頭、去低質(zhì)量、去污染等過程對(duì)原始reads進(jìn)行過濾，過濾的步驟如下：(1)去除含adapter的reads；(2)去除含N比例大于10%的reads；(3)去除低質(zhì)量reads(質(zhì)量值Q≤10的堿基數(shù)占整條read的50%以上)。

1.5 比對(duì)序列的構(gòu)建和定量

利用TOPHAT v2.1.0軟件[7]將clean reads與參考序列進(jìn)行比對(duì)，參考基因組來源于Ensembl，序列為Bostaurus的UMD3.1版本(http://www.ensembl.org.Bos_taurus/ Info/Index/)，并統(tǒng)計(jì)reads與每一個(gè)參考序列的比對(duì)結(jié)果。用RSEM(RNASeq by Expectation Maximization)工具進(jìn)行基因以及轉(zhuǎn)錄本的表達(dá)定量。表達(dá)定量的結(jié)果以FPKM(fragments per kilobase per million fragments)為單位，具體計(jì)算公式如下：

設(shè)FPKM(A)為基因A的表達(dá)量，則C為唯一比對(duì)到基因A的fragments數(shù)，N為唯一比對(duì)到參考基因的總fragments數(shù)，L為基因A編碼區(qū)的堿基數(shù)。

1.6 差異基因篩選

通過Noiseq軟件包[8]計(jì)算出每個(gè)基因的差異倍數(shù)以及偏離度概率值?；虻牟町惐稊?shù)MA=log2((Treat-avg)/(Control-avg),絕對(duì)差值DA=|Control-avg-Treat-avg|，如果MA、DA均明顯偏離噪音背景數(shù)據(jù)集，則該基因?qū)儆诓町惐磉_(dá)基因。按照差異倍數(shù)≥2、同時(shí)偏離概率值≥0.8的標(biāo)準(zhǔn)篩選差異表達(dá)基因。

1.7 GO功能和Pathway通路分析

利用Gene Ontology數(shù)據(jù)庫(http://www.geneontology.org/)對(duì)差異表達(dá)基因進(jìn)行GO功能顯著性富集分析。將整個(gè)基因組作為背景基因，通過超幾何檢驗(yàn)計(jì)算P值，以P值≤0.05為閾值，滿足此條件的GO term定義為差異表達(dá)基因中顯著富集的GO term。

對(duì)差異表達(dá)基因的KEGG信號(hào)通路進(jìn)行分析[9]，以P≤0.05為閾值定義在差異表達(dá)基因中顯著富集的通路，確定差異表達(dá)基因參與的主要生化代謝途徑與信號(hào)轉(zhuǎn)導(dǎo)途徑。

2 結(jié)果與分析

每個(gè)樣品總RNA濃度均大于200 ng·μL-1，RIN值范圍在7.4～8.4，達(dá)到了文庫構(gòu)建的要求，可以進(jìn)行后續(xù)試驗(yàn)。

2.1 Reads數(shù)據(jù)分析

由于原始測(cè)序數(shù)據(jù)可能包含低質(zhì)量序列、接頭序列等，為了保證信息分析結(jié)果的可靠性，將原始reads進(jìn)行分類，分別統(tǒng)計(jì)低質(zhì)量、含接頭、含N reads數(shù)目。經(jīng)perl腳本過濾處理后對(duì)照組個(gè)體及熱應(yīng)激個(gè)體分別得到圖1所示結(jié)果。

A-D，NHS樣本；E-H，HS樣本。A-D， Samples of NHS; E-H， Samples of HS.圖1 原始數(shù)據(jù)過濾統(tǒng)計(jì)Fig.1 Statistics of raw data filtering

2.2 測(cè)序序列與基因組比對(duì)結(jié)果

比對(duì)結(jié)果顯示，NHS組和HS組個(gè)體與參考基因組的比對(duì)率、唯一比對(duì)率較高(表1)，由此推斷各樣本文庫覆蓋率較高，表明轉(zhuǎn)錄組測(cè)序在建庫時(shí)樣品代表性較高，可用于后續(xù)試驗(yàn)分析。

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

熱應(yīng)激組和對(duì)照組之間的差異表達(dá)基因共有96個(gè)，與對(duì)照組相比，熱應(yīng)激組中上調(diào)的基因有46個(gè)，下調(diào)的基因有50個(gè)。熱應(yīng)激組和對(duì)照組基因表達(dá)散點(diǎn)圖見圖2。差異表達(dá)基因中顯著升高的前3個(gè)基因分別為成纖維細(xì)胞生長因子受體1(FGFR1)、顆粒體蛋白(GRN)和線粒體核糖體蛋白L49(MRPL49)，但未見這幾個(gè)基因在熱應(yīng)激方面的文獻(xiàn)報(bào)道，其在奶牛熱應(yīng)激中的作用機(jī)制還需進(jìn)一步研究。此外，與免疫相關(guān)的乳過氧化物酶基因(LPO)，與熱應(yīng)激相關(guān)的熱休克蛋白40基因(Hsp40)表達(dá)顯著升高，而與泌乳相關(guān)的催乳素受體基因(PRLR)和胰島素樣生長因子1基因(IGF1)表達(dá)顯著下降。

2.4 差異表達(dá)基因的GO分類與富集分析

通過對(duì)96個(gè)顯著差異表達(dá)的基因進(jìn)行GO富集分析，包括生物過程(biological process)、細(xì)胞組分(cellular component)以及分子功能(molecular function)。GO分類分析結(jié)果顯示，“細(xì)胞組分”共注釋到41個(gè)基因，“分子功能”共注釋到38個(gè)基因，“生物學(xué)過程”共注釋到36個(gè)基因(圖3)。其中富集最多的生物過程是細(xì)胞過程，細(xì)胞組分富集主要集中在細(xì)胞及細(xì)胞組分，分子功能富集最為顯著的為結(jié)合功能。

表1 測(cè)序序列與參考基因的比對(duì)統(tǒng)計(jì)

Table 1 Statistics of sequencing reads aligned to reference genes

分組Group樣品編號(hào)Sample No.總讀段數(shù)目Total reads總堿基對(duì)Total base-pairs(占比Proportion)總比對(duì)讀段Total mapped reads(占比Proportion)唯一比對(duì)讀長Unique match(占比Proportion)非熱應(yīng)激組A619400546194005400(100%)48106566(77.67%)44514720(71.87%)NHS groupB663271186632711800(100%)49767554(75.03%)46061857(69.45%)C559933365599333600(100%)43522351(77.73%)40297755(71.79%)D5822090458220904(100%)43982632(75.54%)40510113(69.58%)熱應(yīng)激組E667136886671368800(100%)50245709(75.32%)46694358(69.99%)HS groupF615957866159578600(100%)43909443(71.29%)40986696(66.54%)G573263445732634400(100%)44326144(77.32%)40918820(71.38%)H487264804872648000(100%)35720407(73.31%)33210434(68.16%)

黃色表示顯著上調(diào)基因，藍(lán)色表示顯著下調(diào)基因，灰色表示非顯著性改變基因。Yellow indicated significantly up-regulated genes， blue indicated significantly down-regulated genes and gray indicated non-significantly changed genes.圖2 基因表達(dá)散點(diǎn)圖Fig.2 Scattered plot of gene expression

2.5 差異表達(dá)基因的KEGG富集分析

通過計(jì)算KEGG各功能類別的P值來對(duì)差異表達(dá)基因進(jìn)行KEGG富集分析，富集條件為P<0.05。P值越小說明差異表達(dá)基因在該類別中富集越顯著。經(jīng)KEGG分析后共富集到18條信號(hào)通路，其中6條與疾病有關(guān)，9條與代謝有關(guān)。另外，與免疫應(yīng)答相關(guān)的細(xì)胞因子-細(xì)胞因子受體相互作用信號(hào)通路和NOD樣受體信號(hào)通路明顯富集(表2、圖4)。

3 討論

乳腺為奶牛泌乳的重要器官，熱應(yīng)激對(duì)乳腺的直接影響是高溫引起產(chǎn)奶量下降的重要原因[10]。有機(jī)體的生長、發(fā)育和存活都依賴于細(xì)胞增殖和凋亡的平衡。熱應(yīng)激持續(xù)作用于泌乳奶牛，會(huì)對(duì)機(jī)體產(chǎn)生不良的累積反應(yīng)。高溫抑制奶牛乳腺上皮細(xì)胞正常生長，促使細(xì)胞凋亡，影響乳蛋白合成[11-12]。干奶期的奶牛發(fā)生熱應(yīng)激可使分娩前乳腺發(fā)育減緩，導(dǎo)致下一個(gè)泌乳期的產(chǎn)奶量降低[13]。與非熱應(yīng)激期相比，泌乳前、中、后期的奶牛在整個(gè)熱應(yīng)激期的平均產(chǎn)奶量分別下降15.41%、12.56%、11.87%[14]。本研究以春夏季中國荷斯坦牛的乳腺組織為研究對(duì)象，探討其在奶牛發(fā)生熱應(yīng)激過程中發(fā)揮的作用及分子機(jī)制。選取的奶牛個(gè)體之間沒有親緣關(guān)系，可避免親緣關(guān)系引起的基因表達(dá)差異；熱應(yīng)激組和非熱應(yīng)激組奶牛年齡相同，可避免因年齡差異引起的基因表達(dá)差異。

圖3 差異表達(dá)基因的GO分類Fig.3 GO classification of differentially expressed genes

表2 熱應(yīng)激奶牛乳腺組織差異表達(dá)基因KEGG富集分析結(jié)果

Table 2 KEGG enrichment analysis of DEG between the mammary glands of heat stressed and non-heat stressed dairy cattle

通路編號(hào)Pathway ID通路描述DescriptionP值P-value差異表達(dá)基因Names of DEGsko00232咖啡因代謝Caffeine metabolism0.000222LOC540707ko04060細(xì)胞因子-細(xì)胞因子受體相互作用Cytokine-cytokine receptor interaction0.000414IL6ST、CCL2、IL1B、CCL20、CCL8、PRLR、IL6ko04978礦物質(zhì)吸收Mineral absorption0.000532S100A8、S100A9(transcript variant X2)、S100A9(transcript variant X3)、S100A12ko04970唾液分泌Salivary secretion0.001677LOC617219、GP2(transcript variant X1)、GP2、LPO(tran-script variant X1)、LPO(transcript variant X2)ko05323類風(fēng)濕性關(guān)節(jié)炎R(shí)heumatoid arthritis0.002156CCL2、IL1B、CCL20、IL6ko00040戊糖和葡萄糖醛酸酯相互轉(zhuǎn)換Pentose and glucuronateinter conversions0.005388UGP4(transcript variant X4)、UGP2(transcript variant X2)

續(xù)表2 Continued Table 2

PRLR和IGF1能參與奶牛泌乳調(diào)控。PRL與其受體PRLR的結(jié)合可以激活酪氨酸激酶2，催化信號(hào)轉(zhuǎn)導(dǎo)和轉(zhuǎn)錄活化蛋白5(STAT5)磷酸化，進(jìn)而調(diào)控乳蛋白基因的表達(dá)[15]。PRLR基因位點(diǎn)的遺傳多樣性與牛乳中脂肪酸組成相關(guān)[16]。血液IGF-1及其結(jié)合蛋白濃度升高，可有效促進(jìn)奶牛的泌乳活動(dòng)[17]。IGFl可介導(dǎo)生長激素啟動(dòng)JAK-STAT信號(hào)通路調(diào)節(jié)泌乳[18]。熱應(yīng)激對(duì)血漿IGF1濃度沒有顯著影響，但熱應(yīng)激降低了肝臟IGF1量[19]，表明在熱應(yīng)激狀態(tài)，肝臟生長激素-類胰島素生長因子軸可能發(fā)生了解偶聯(lián)作用。本研究發(fā)現(xiàn)，96個(gè)基因在熱應(yīng)激組和非熱應(yīng)激組奶牛中發(fā)生差異表達(dá)，其中PRLR和IGF1表達(dá)降低，說明熱應(yīng)激導(dǎo)致奶牛產(chǎn)奶量下降可能與PRLR和IGF1表達(dá)下降有關(guān)。

熱應(yīng)激除了降低奶牛產(chǎn)奶性能，還會(huì)抑制體液和細(xì)胞免疫，易發(fā)生乳房炎、子宮內(nèi)膜炎等疾病。FGFR1為酪氨酸激酶受體家族中的一員，研究表明，F(xiàn)GFR1在乳腺癌、乳腺浸潤性導(dǎo)管癌、肺鱗狀細(xì)胞癌及大腸癌等腫瘤中異常表達(dá)[20]，使細(xì)胞增殖調(diào)節(jié)發(fā)生紊亂，激活細(xì)胞內(nèi)多種信號(hào)通路，導(dǎo)致疾病發(fā)生[21]。本研究中，熱應(yīng)激組奶牛FGFR1基因表達(dá)顯著升高，暗示熱應(yīng)激組奶牛易發(fā)生疾病。LPO作為抗菌肽具有免疫調(diào)節(jié)功能，為乳腺非特異性抵抗細(xì)菌系統(tǒng)之一，對(duì)革蘭氏陽性菌和陰性菌均有抑制或殺滅作用。對(duì)副結(jié)核分枝桿菌感染牛進(jìn)行唾腺轉(zhuǎn)錄組分析，結(jié)果發(fā)現(xiàn)LPO表達(dá)發(fā)生變化[22]。Hsp蛋白為廣泛存在于原核和真核生物中受應(yīng)激誘導(dǎo)表達(dá)的高度保守性蛋白，又稱為伴侶蛋白。按相對(duì)分子質(zhì)量不同，可分為Hsp40、Hsp60、Hsp70和Hsp90等。Hsp40蛋白除了作為分子伴侶還具有多種功能，近些年發(fā)現(xiàn)Hsp40蛋白在先天性免疫中具有重要作用[23]。本研究中，LPO和Hsp40表達(dá)顯著升高，這些結(jié)果表明，熱應(yīng)激影響了奶牛的免疫機(jī)能。

一般來說，生物體內(nèi)的重要生理過程需要多個(gè)基因共同參與，發(fā)揮相關(guān)功能，進(jìn)行一系列調(diào)控。本研究差異表達(dá)基因GO富集分析中，細(xì)胞組分、生物學(xué)過程和分子功能類別中富集了大量基因。KEGG通路分析可進(jìn)一步了解這些基因的生物學(xué)功能。因此，本研究在GO分析的基礎(chǔ)上又進(jìn)行了信號(hào)通路分析，發(fā)現(xiàn)6條疾病相關(guān)通路明顯富集，包括細(xì)胞因子-細(xì)胞因子受體相互作用信號(hào)通路。細(xì)胞因子-細(xì)胞因子受體的相互作用可調(diào)節(jié)細(xì)胞的生長分化，參與免疫、炎癥反應(yīng)。奶牛熱應(yīng)激過程中，細(xì)胞因子-細(xì)胞因子受體的相互作用發(fā)揮了重要的調(diào)控作用。該通路中顯著富集的PRLR與其配體催乳素(PRL)的相互作用，被認(rèn)為可作為奶牛對(duì)熱應(yīng)激生理反應(yīng)的介導(dǎo)因子。另外，NOD樣受體信號(hào)通路被顯著富集。NOD樣受體為一類位于細(xì)胞質(zhì)的模式識(shí)別受體，在先天性免疫應(yīng)答中起著十分重要的作用。NOD樣受體被激活后，通過一系列的信號(hào)通路，可誘導(dǎo)各種炎癥因子的釋放。水牛發(fā)生熱應(yīng)激后，該通路中IL1B和IL6表達(dá)發(fā)生明顯變化[24]。犢牛熱應(yīng)激后，血液中NOD樣受體信號(hào)通路被顯著富集[25]，說明熱應(yīng)激對(duì)奶牛免疫應(yīng)答的影響可能與細(xì)胞因子-細(xì)胞因子受體相互作用信號(hào)通路和NOD樣受體信號(hào)通路有關(guān)，還需進(jìn)一步實(shí)驗(yàn)驗(yàn)證。

本研究利用RNA-Seq測(cè)序技術(shù)研究熱應(yīng)激對(duì)中國荷斯坦牛乳腺組織基因表達(dá)的影響，共發(fā)現(xiàn)96個(gè)差異表達(dá)基因，初步揭示了熱應(yīng)激奶牛基因表達(dá)差異。對(duì)差異表達(dá)基因進(jìn)行通路分析發(fā)現(xiàn)，6條通路與疾病相關(guān)，另外，與免疫應(yīng)答相關(guān)的細(xì)胞因子-細(xì)胞因子受體相互作用和NOD樣受體信號(hào)通路明顯富集，暗示其在奶牛熱應(yīng)激免疫應(yīng)答中可能發(fā)揮重要作用，為進(jìn)一步了解奶牛熱應(yīng)激的免疫反應(yīng)分子機(jī)制提供了參考。