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

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

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

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

1 材料與方法

1.1 試驗材料

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

1.2 總RNA提取和鑒定

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

1.3 文庫制備和測序

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

1.4 原始序列的處理

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

1.5 比對序列的構建和定量

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

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

1.6 差異基因篩選

通過Noiseq軟件包[8]計算出每個基因的差異倍數(shù)以及偏離度概率值?；虻牟町惐稊?shù)MA=log2((Treat-avg)/(Control-avg),絕對差值DA=|Control-avg-Treat-avg|，如果MA、DA均明顯偏離噪音背景數(shù)據集，則該基因屬于差異表達基因。按照差異倍數(shù)≥2、同時偏離概率值≥0.8的標準篩選差異表達基因。

1.7 GO功能和Pathway通路分析

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

對差異表達基因的KEGG信號通路進行分析[9]，以P≤0.05為閾值定義在差異表達基因中顯著富集的通路，確定差異表達基因參與的主要生化代謝途徑與信號轉導途徑。

2 結果與分析

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

2.1 Reads數(shù)據分析

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

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

2.2 測序序列與基因組比對結果

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

2.3 基因差異表達分析

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

2.4 差異表達基因的GO分類與富集分析

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

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

Table 1 Statistics of sequencing reads aligned to reference genes

分組Group樣品編號Sample No.總讀段數(shù)目Total reads總堿基對Total base-pairs(占比Proportion)總比對讀段Total mapped reads(占比Proportion)唯一比對讀長Unique match(占比Proportion)非熱應激組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%)熱應激組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%)

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

2.5 差異表達基因的KEGG富集分析

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

3 討論

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

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

表2 熱應激奶牛乳腺組織差異表達基因KEGG富集分析結果

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

通路編號Pathway ID通路描述DescriptionP值P-value差異表達基因Names of DEGsko00232咖啡因代謝Caffeine metabolism0.000222LOC540707ko04060細胞因子-細胞因子受體相互作用Cytokine-cytokine receptor interaction0.000414IL6ST、CCL2、IL1B、CCL20、CCL8、PRLR、IL6ko04978礦物質吸收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類風濕性關節(jié)炎Rheumatoid arthritis0.002156CCL2、IL1B、CCL20、IL6ko00040戊糖和葡萄糖醛酸酯相互轉換Pentose and glucuronateinter conversions0.005388UGP4(transcript variant X4)、UGP2(transcript variant X2)

續(xù)表2 Continued Table 2

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

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

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

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