李秋玲，齐 颖，王 琛，张一名，王新妤，尚校兰，贾永红，李美茹，储明星

(1.廊坊师范学院 生命科学学院，河北省动物多样性重点实验室，廊坊市细胞工程与应用研究重点实验室，河北 廊坊 065000； 2.中国农业科学院 北京畜牧兽医研究所，农业农村部动物遗传育种与繁殖(家禽)重点实验室，北京 100193)

热应激是动物机体受到环境高温刺激所产生的一系列非特异性应答反应的总和[1]。随着温室效应加剧，夏秋季高温气候出现频率逐渐升高，奶牛热应激问题日益凸显，严重影响奶业生产，成为奶牛夏季生产的一个重大难题。中国荷斯坦牛原产于欧洲，耐寒性较强，耐热性较差，夏季容易发生热应激，轻者引起新陈代谢紊乱[2]，生产性能、乳品质[3]、繁殖力、免疫力下降[4]，生产年限缩短，重者机体功能衰竭导致死亡，给奶业生产带来巨大的经济损失。据统计，美国奶业每年因热应激造成的损失高达9亿美元[5]。

近年来，我国奶牛单产不断提高，但产奶量与耐热程度高度负相关，另外，选种时主要考虑生产性能使奶牛的耐热性不断下降。随着动物营养学的快速发展，阻碍生产性能提高的因素已不再是日粮营养供应问题，更多的是奶牛饲养环境的限制。热应激作为一个重要的环境因素已引起广大学者的高度关注，热应激发生时动物体内各种基因调控的关系有待深入研究，从而揭示热应激的作用机制和规律。本研究拟通过Illumina测序平台对奶牛热应激期和非热应激期乳腺组织转录组文库进行高通量测序和生物信息学分析，研究热应激奶牛乳腺组织差异表达基因及其调控通路，为进一步研究奶牛发生热应激的分子机制，耐热奶牛分子育种提供理论依据。

1 材料与方法

1.1 试验材料

根据系谱资料、奶牛生产性能(DHI)测定数据和牛场记录数据选取性别相同，饲养环境、年龄及胎次相同的8头中国荷斯坦牛作为实验牛。根据本实验室前期研究对奶牛热应激和非热应激期生理参数比较结果[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的碱基数占整条read的50%以上)。

1.5 比对序列的构建和定量

利用TOPHAT v2.1.0软件[7]将clean reads与参考序列进行比对，参考基因组来源于Ensembl，序列为Bostaurus的UMD3.1版本(http://www.ensembl.org.Bos_taurus/ Info/Index/)，并统计reads与每一个参考序列的比对结果。用RSEM(RNASeq by Expectation Maximization)工具进行基因以及转录本的表达定量。表达定量的结果以FPKM(fragments per kilobase per million fragments)为单位，具体计算公式如下：

设FPKM(A)为基因A的表达量，则C为唯一比对到基因A的fragments数，N为唯一比对到参考基因的总fragments数，L为基因A编码区的碱基数。

1.6 差异基因筛选

通过Noiseq软件包[8]计算出每个基因的差异倍数以及偏离度概率值。基因的差异倍数MA=log2((Treat-avg)/(Control-avg),绝对差值DA=|Control-avg-Treat-avg|，如果MA、DA均明显偏离噪音背景数据集，则该基因属于差异表达基因。按照差异倍数≥2、同时偏离概率值≥0.8的标准筛选差异表达基因。

1.7 GO功能和Pathway通路分析

利用Gene Ontology数据库(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，达到了文库构建的要求，可以进行后续试验。

2.1 Reads数据分析

由于原始测序数据可能包含低质量序列、接头序列等，为了保证信息分析结果的可靠性，将原始reads进行分类，分别统计低质量、含接头、含N reads数目。经perl脚本过滤处理后对照组个体及热应激个体分别得到图1所示结果。

A-D，NHS样本；E-H，HS样本。A-D， Samples of NHS; E-H， Samples of HS.图1 原始数据过滤统计Fig.1 Statistics of raw data filtering

2.2 测序序列与基因组比对结果

比对结果显示，NHS组和HS组个体与参考基因组的比对率、唯一比对率较高(表1)，由此推断各样本文库覆盖率较高，表明转录组测序在建库时样品代表性较高，可用于后续试验分析。

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 测序序列与参考基因的比对统计

Table 1 Statistics of sequencing reads aligned to reference genes

分组Group样品编号Sample No.总读段数目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]。有机体的生长、发育和存活都依赖于细胞增殖和凋亡的平衡。热应激持续作用于泌乳奶牛，会对机体产生不良的累积反应。高温抑制奶牛乳腺上皮细胞正常生长，促使细胞凋亡，影响乳蛋白合成[11-12]。干奶期的奶牛发生热应激可使分娩前乳腺发育减缓，导致下一个泌乳期的产奶量降低[13]。与非热应激期相比，泌乳前、中、后期的奶牛在整个热应激期的平均产奶量分别下降15.41%、12.56%、11.87%[14]。本研究以春夏季中国荷斯坦牛的乳腺组织为研究对象，探讨其在奶牛发生热应激过程中发挥的作用及分子机制。选取的奶牛个体之间没有亲缘关系，可避免亲缘关系引起的基因表达差异；热应激组和非热应激组奶牛年龄相同，可避免因年龄差异引起的基因表达差异。

图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类风湿性关节炎Rheumatoid arthritis0.002156CCL2、IL1B、CCL20、IL6ko00040戊糖和葡萄糖醛酸酯相互转换Pentose and glucuronateinter conversions0.005388UGP4(transcript variant X4)、UGP2(transcript variant X2)

续表2 Continued Table 2

PRLR和IGF1能参与奶牛泌乳调控。PRL与其受体PRLR的结合可以激活酪氨酸激酶2，催化信号转导和转录活化蛋白5(STAT5)磷酸化，进而调控乳蛋白基因的表达[15]。PRLR基因位点的遗传多样性与牛乳中脂肪酸组成相关[16]。血液IGF-1及其结合蛋白浓度升高，可有效促进奶牛的泌乳活动[17]。IGFl可介导生长激素启动JAK-STAT信号通路调节泌乳[18]。热应激对血浆IGF1浓度没有显著影响，但热应激降低了肝脏IGF1量[19]，表明在热应激状态，肝脏生长激素-类胰岛素生长因子轴可能发生了解偶联作用。本研究发现，96个基因在热应激组和非热应激组奶牛中发生差异表达，其中PRLR和IGF1表达降低，说明热应激导致奶牛产奶量下降可能与PRLR和IGF1表达下降有关。

热应激除了降低奶牛产奶性能，还会抑制体液和细胞免疫，易发生乳房炎、子宫内膜炎等疾病。FGFR1为酪氨酸激酶受体家族中的一员，研究表明，FGFR1在乳腺癌、乳腺浸润性导管癌、肺鳞状细胞癌及大肠癌等肿瘤中异常表达[20]，使细胞增殖调节发生紊乱，激活细胞内多种信号通路，导致疾病发生[21]。本研究中，热应激组奶牛FGFR1基因表达显著升高，暗示热应激组奶牛易发生疾病。LPO作为抗菌肽具有免疫调节功能，为乳腺非特异性抵抗细菌系统之一，对革兰氏阳性菌和阴性菌均有抑制或杀灭作用。对副结核分枝杆菌感染牛进行唾腺转录组分析，结果发现LPO表达发生变化[22]。Hsp蛋白为广泛存在于原核和真核生物中受应激诱导表达的高度保守性蛋白，又称为伴侣蛋白。按相对分子质量不同，可分为Hsp40、Hsp60、Hsp70和Hsp90等。Hsp40蛋白除了作为分子伴侣还具有多种功能，近些年发现Hsp40蛋白在先天性免疫中具有重要作用[23]。本研究中，LPO和Hsp40表达显著升高，这些结果表明，热应激影响了奶牛的免疫机能。

一般来说，生物体内的重要生理过程需要多个基因共同参与，发挥相关功能，进行一系列调控。本研究差异表达基因GO富集分析中，细胞组分、生物学过程和分子功能类别中富集了大量基因。KEGG通路分析可进一步了解这些基因的生物学功能。因此，本研究在GO分析的基础上又进行了信号通路分析，发现6条疾病相关通路明显富集，包括细胞因子-细胞因子受体相互作用信号通路。细胞因子-细胞因子受体的相互作用可调节细胞的生长分化，参与免疫、炎症反应。奶牛热应激过程中，细胞因子-细胞因子受体的相互作用发挥了重要的调控作用。该通路中显著富集的PRLR与其配体催乳素(PRL)的相互作用，被认为可作为奶牛对热应激生理反应的介导因子。另外，NOD样受体信号通路被显著富集。NOD样受体为一类位于细胞质的模式识别受体，在先天性免疫应答中起着十分重要的作用。NOD样受体被激活后，通过一系列的信号通路，可诱导各种炎症因子的释放。水牛发生热应激后，该通路中IL1B和IL6表达发生明显变化[24]。犊牛热应激后，血液中NOD样受体信号通路被显著富集[25]，说明热应激对奶牛免疫应答的影响可能与细胞因子-细胞因子受体相互作用信号通路和NOD样受体信号通路有关，还需进一步实验验证。

本研究利用RNA-Seq测序技术研究热应激对中国荷斯坦牛乳腺组织基因表达的影响，共发现96个差异表达基因，初步揭示了热应激奶牛基因表达差异。对差异表达基因进行通路分析发现，6条通路与疾病相关，另外，与免疫应答相关的细胞因子-细胞因子受体相互作用和NOD样受体信号通路明显富集，暗示其在奶牛热应激免疫应答中可能发挥重要作用，为进一步了解奶牛热应激的免疫反应分子机制提供了参考。