基于RNA测序分析热应激奶牛乳腺组织中miRNA在乳蛋白合成调控中的潜在作用

高胜涛 马 露 卜登攀

(中国农业科学院北京畜牧兽医研究所，动物营养学国家重点实验室，北京 100193)

热应激(heat stress，HS)影响奶业生产[1]和奶牛健康[2]。热应激可降低奶牛干物质采食量(DMI)、产奶量和乳蛋白产量[3-6]，但DMI下降只能解释30%～50%的产奶量下降[7]。本团队前期研究表明，热应激诱导的DMI下降只能解释65%的乳蛋白产量下降[5]，从而表明热应激对乳腺组织的乳蛋白合成具有直接的影响。

乳腺组织的乳蛋白合成前体物供应减少是热应激引起乳蛋白合成减少的重要原因之一[5]。对乳腺组织进行RNA测序(RNAseq)分析发现，热应激可降低乳腺组织中乳蛋白合成相关编码基因、转录翻译关键调控基因及氨基酸和葡萄糖转运体相关基因的表达；对所有差异表达基因进行生物信息学分析表明，乳腺组织中代谢相关通路(尤其是碳水化合物及脂代谢相关通路)均被热应激抑制；而上游调节因子分析显示炎症反应相关通路被激活，表明炎症反应诱导驱动的乳腺组织中整体代谢的抑制，可能也是热应激时乳蛋白合成减少的原因之一[8]。

乳腺组织中miRNA是存在数量最多的miRNA，可以在转录后水平通过介导mRNA的翻译或同源mRNA的降解调节靶基因的表达[9-11]。大量研究表明，miRNA在细胞生长与分化、增殖与凋亡的生物过程中发挥重要的作用[12-13]。此外，miRNA还参与了包括氧化应激、低氧应激等多种应激反应[14-15]。Li等[16]在热应激荷斯坦奶牛乳腺组织中鉴定出27个相对于非热应激牛差异表达的miRNA。然而关于miRNA表达谱的变化在热应激诱导的乳蛋白合成下降中的作用，目前仍未见报道。因此，本研究采用RNAseq分析了热应激对乳腺组织中miRNA表达谱变化的影响，并阐述了差异表达的miRNA在乳蛋白合成调控中的作用，为进一步验证热应激条件下miRNA对乳蛋白合成的调控作用奠定理论科学基础。

1 材料与方法

1.1 试验设计

本试验所用奶牛乳腺组织样品取自前期已发表试验[5,8]。选取4头健康经产荷斯坦泌乳奶牛[泌乳日龄：(101±10) d；体重：(574±36) kg；产奶量：(38±2) kg/d；2胎，怀孕1～2个月]，采用2×2交叉试验设计。本试验在中国农业科学院北京畜牧兽医研究所动物营养学国家重点实验室昌平基地4个大家畜人工环境控制舱(Kooland，北京库蓝科技有限公司；4.0 m×3.0 m×2.5 m；温度15～40 ℃；相对湿度25%～85%；光照0～800 lx，连续可调)内进行。试验共分2期，每期包含预试期和试验期各9 d，2期间隔30 d。预试期所有试验动物饲养于热中性环境[温度：20 ℃；相对湿度：55%；温热指数(THI)：65.5]，且自由采食。试验期配对限饲组(PFTN组，n=4)奶牛继续饲养于热中性环境，热应激组(HS组，n=4)奶牛饲养于热应激环境(06：00—18：00 h：36 ℃，18：00—06：00 h：32 ℃；相对湿度：40%；THI=84.0/79.2)。THI计算参考Buffington等[17]的方法。THI设置参考热应激阈值[18]。PFTN组奶牛的某天饲喂量(采食量)占该组牛预试期DMI平均值百分比与热应激组此前1天DMI占该组牛预试期DMI平均值百分比(26.8%～61.0%)保持一致，即PFTN组相对于HS组错后1 d开始和结束试验以及采集样品，保证试验期内HS组和PFTN组每天DMI保持一致[5,8]。

1.2 乳腺组织样品采集

于每期的试验期第10天采集乳腺组织样品，详细采集过程参照Gao等[8]和Bu等[19]所述方法。在样品采集前，使用2%盐酸利多卡因对乳腺组织的左后区或右后区进行皮下局部麻醉。样品采集完成后使用灭菌纱布进行压迫止血，直至切口不再有明显的血液流出，随后使用米歇尔夹(11 mm；Henry Schein公司，美国)进行皮肤快速封闭。样品采集后连续7 d监测奶牛直肠温度，产奶量及采食量。组织样品采集后立即放液氮中冷冻保存。

1.3 总RNA提取及miRNA测序

1.4 数据分析

首先移除原始数据中包含poly-N、5′接头污染、无3′接头及低质量的序列。同时计算原始数据的Q20、Q30和GC碱基含量。移除3′端接头后，使用miRdeep2(2.0.0.5)[21]对长度为18～30 nt的miRNA序列进行提取，并使用默认参数与miRbase(20.0)参考miRNA进行比对识别已知miRNA及新miRNA。使用Bowtie对测序所得小RNA序列与牛的基因组(Bos taurus Bos_taurus_UMD_3.1.1)比对[22]，识别鉴定并移除miRNA中的tRNAs、rRNAs、snoRNAs。另外，只能映射到基因组重复位点的miRNA也被移除。使用miRanda[23]、PITA[24]、和RNAhybrid[25]对已知miRNA和新miRNA进行靶基因预测。并选择2种或2种以上软件中都成功预测的靶基因作为本研究中miRNA的靶基因，并用于下游分析。

使用DESeq2[26]进行miRNA的差异分析。对P值进行Benjamin & Hochberg[27]矫正，并使用矫正后P值(Padj)<0.05作为差异miRNA的阈值。使用DAVID(用于注释、可视化和集成发现的数据库，Database for Annotation, Visualization, and Integrated Discovery，6.7；http://david.abcc.ncifcrf.gov/)[28]，对差异表达miRNA的靶基因进行基因本体(gene ontology，GO)生物过程和京都基因与基因组百科全书(Kyoto encyclopedia of genes and genomes，KEGG)的功能注释。

2 结 果

2.1 测序数据产出情况

本研究中每个组共3个文库构建成功，并平均产出24 847 562个错误率低于0.01%的序列。所有文库的干净数据中，配对碱基质量分数高于20和30(Q20和Q30)的平均比率分别为98.52%和96.83%，平均GC含量为48.72%(表1)。

表1 HS组和PFTN组奶牛乳腺组织miRNA测序产出数据质量Table 1 miRNA sequencing output data quality of dairy cows’ mammary tissue in HS group and PFTN group

2.2 已知miRNA和新miRNA识别及其靶基因预测

如表2所示，本研究共鉴定识别出已知miRNA 578个，新miRNA 139个。所有已知miRNA和新miRNA的潜在靶基因如补充材料1(https://doi.org/10.6084/m9.figshare.14315609.v1)中所示。

表2 已知miRNA和新miRNA与参考miRNA比对情况汇总Table 2 Statistics of known miRNA and novel miRNA predicting mapping to reference miRNA

续表2项目Items样品Samples映射的成熟体Mapped mature映射的发夹结构Mapped hairpin映射的单一小RNAMapped unique sRNA映射的所有小RNAMapped total sRNA新miRNANovel miRNA合计 Total1391441 0614 664HS16674138505HS28292166734HS393101217654PFTN17883141629PFTN28192174872PFTN3951042251 270

2.3 差异表达miRNA

如图1及补充材料2(https://doi.org/10.6084/m9.figshare.14315609.v1)所示，在本研究识别出的578个已知miRNA和139个新miRNA中，共16个miRNA在HS组奶牛乳腺中相对于PFTN组差异表达(Padj<0.05)。其中，12个miRNA在HS组奶牛乳腺组织中相对于PFTN组表达上调，4个miRNA在HS组奶牛PFTN乳腺组织中相对于奶牛表达下调(表3)。16个差异表达的miRNA在每个样品中的表达情况如图2所示。

红点代表在HS组奶牛乳腺中相对于PFTN组表达上调的miRNA，蓝点代表HS组奶牛乳腺中相对于PFTN组表达下调的miRNA，灰点代表HS组和PFTN组之间奶牛乳腺中没有显著上调或下调的miRNA。差异表达miRNA选择的阈值是Padj<0.05。Red points represent the up regulated miRNA in dairy cows mammary tissue of HS group compared with PFTN group, blue points represent the down regulated miRNA in dairy cows mammary tissue of HS group compared with PFTN group, grey points represent the not significantly up or down regulated miRNA in dairy cows mammary tissue between HS group and PFTN group. The selected threshold of differently expressed miRNA was Padj<0.05.图1 HS组和PFTN组奶牛乳腺组织中差异表达miRNA火山图Fig.1 Volcano plot of differently expressed miRNA in dairy cows’ mammary tissue between HS group and PFTN group

图2 差异表达miRNA在HS组和PFTN组奶牛乳腺组织中的热图Fig.2 Heatmap of differently expressed miRNA in dairy cows’ mammary tissue between HS group and PFTN group

表3 16个HS组奶牛乳腺组织中相对PFTN组差异表达miRNATable 3 Sixteen differently expressed miRNA in dairy cows’ mammary tissue of HS group compared with PFTN group

2.4 差异表达miRNA的靶基因的功能富集

本研究使用DAVID对下调和上调miRNA的靶基因分别进行功能分析(图3，P<0.05)。其中表达下调的miRNA的靶基因富集的GO生物过程包括：铁离子转运、典型Wnt信号通路、味觉、细胞内铁离子螯合、细胞增殖、核糖体大亚基组装、白细胞介素2生成负调控、蛋白质聚合、血小板活化、刺激响应、鸟苷三磷酸酶(GTPase)激活和氧化还原过程；KEGG通路主要包括：吞噬体、Ras信号通路、AMP激活的蛋白质激酶(AMPK)信号通路、赖氨酸降解、金黄色葡萄球菌感染、核糖体、趋化因子信号通路、醚类脂质代谢和氨基糖和核苷酸糖代谢。下调的miRNA的靶基因富集的GO生物过程包括：细胞迁移的正向调节、细胞对肿瘤坏死因子的响应、氧自由基响应、白细胞介素-1响应、呼吸气体交换、转化生长因子β受体信号通路、mRNA异化正向调控、干扰素-γ生成正相关、上皮细胞调控正向调控、细胞因子调控、血管内皮生长因子生成的正向调控、甲状腺激素刺激的细胞反应和三羧酸循环；KEGG通路主要包括：AMPK信号通路、代谢通路、溶酶体、柠檬酸循环、胰岛素信号通路、鞘脂类信号通路、哺乳动物雷帕霉素靶蛋白(mTOR)信号通路、叉头转录因子(FoxO)信号通路、肿瘤坏死因子信号通路、ErbB信号通路、胰岛素拮抗以及脂肪细胞因子信号通路。

GO_BP：基因本体生物过程；KEGG：京都基因与基因组百科全书；miRNA：差异表达miRNA的变化方向，其中“↑”代表表达上调，“↓”代表表达下调。绿色区域为被下调miRNA的靶基因富集的GO生物过程和KEGG通路，红色区域为被上调miRNA的靶基因富集的GO生物过程和KEGG通路。GO_BP: gene ontology biology process; KEGG: Kyoto encyclopedia of genes and genomes; miRNA: the change direction of the differently expressed miRNA, “↑” represents up regulated, “↓” represents down regulated. The green areas were GO biology process KEGG pathway of enriched target genes of down regulated miRNAs the red areas were GO biology process KEGG pathway of enriched target genes of up regulated miRNA.图3 HS组和PFTN组奶牛乳腺组织中差异表达miRNA靶基因的GO和KEGG功能富集分析Fig.3 GO and KEGG functional enrichment analysis of target genes of differently expressed miRNA in dairy cows’ mammary tissue between HS group and PFTN group

此外，本研究还对HS组奶牛乳腺组织中每个相对于PFTN组差异表达的miRNA的靶基因分别进行了功能富集分析(图4)。在16个差异表达的miRNA中，除bta-miR-184、bta-miR-151-3p、bta-miR-3578、bta-miR-451外，其余12个差异表达的miRNA的靶基因均通过DAVID成功富集(P<0.05)。bta-miR-132的靶基因显著富集的GO生物过程为吞噬的正向调控和炎症反应的负向调控。bta-miR-136的靶基因显著富集的KEGG通路为内分泌及其他因子调节的钙离子的重吸收。bta-miR-147的靶基因显著富集的GO生物过程为饥饿刺激的细胞响应，KEGG通路为癌症中的蛋白聚糖。bta-miR-218的靶基因显著富集的GO生物过程为细胞质钙离子浓度的正向调控。bta-miR-885的靶基因显著富集的GO生物过程为转录和DNA复制的调控，KEGG通路为糖酵解和糖异生。

GO_BP：基因本体生物过程；KEGG：京都基因与基因组百科全书。蓝条为下调的miRNA的靶基因所富集的GO生物过程和KEGG通路，红条为上调的miRNA的靶基因富集的GO生物过程和KEGG通路。黑色圆圈代表富集到相应生物过程或通路的基因的数量，圆圈越大，代表富集的基因数量越多。GO_BP: gene ontology biology process; KEGG: Kyoto encyclopedia of genes and genomes. The blue bars were the GO biology process and KEGG pathway enriched by the target genes of the down regulated miRNA, while the red bars were the GO biology process and KEGG pathway enriched by the target genes of the up regulated miRNA. The black circles represented the number of target genes enriched in the corresponding biological processes or pathways, and the larger the circles, the more genes were enriched.图4 HS组和PFTN组奶牛乳腺组织中每个差异表达miRNA的靶基因富集的GO生物过程和KEGG通路Fig.4 GO biology process and KEGG pathway of target genes of each differently expressed miRNA in dairy cows’ mammary tissue between HS group and PFTN group

bta-miR-365-5p的靶基因显著富集的GO生物过程为肽酶活性的负调节，KEGG通路为弓形虫病、胰腺癌、乙型肝炎、FoxO信号通路和慢性髓细胞白血病。bta-miR-382的靶基因显著富集的GO生物过程为变形生长因子β受体信号通路和糖异生，KEGG通路为胰岛素拮抗、胰高血糖素信号通路、糖酵解/糖异生、AMPK信号通路、脂肪细胞因子信号通路。bta-miR-502a的靶基因显著富集的KEGG通路为破骨细胞分化。bta-miR-493的靶基因显著富集的GO生物过程为骨骼肌纤维发育、蛋白水解、组蛋白H3-K9三甲基化负调控、单核细胞趋化、淋巴细胞趋化、白细胞介素-1的细胞反应、干扰素-γ的细胞反应、细胞氧化剂解毒，KEGG通路为磷脂酰肌醇信号系统、代谢通路、抗生素的生物合成、AMPK信号通路。bta-miR-338的靶基因显著富集的GO生物过程为蛋白质四聚化和肽或多糖抗原通过主要组织相容性复合体(MHC)Ⅱ类抗原的处理和呈递，KEGG通路为肺结核、弓形体病、甲型流感、人类T细胞白血病病毒1型(HTLV-1)感染、移植物抗宿主病和哮喘。bta-miR-424-3p的靶基因显著富集的GO生物过程为跨膜运输、线粒体转录启动子、核糖体rRNA导入线粒体、辐射响应、激素响应、复制叉保护、细胞内调节pH、血管内皮生长因子生成的正调控、聚腺苷酸化依赖的mRNA异化过程、DNA结合的负调控和胚胎发育，KEGG通路为小细胞肺癌和嘧啶代谢。bta-miR-6516的靶基因显著富集的GO生物过程为感官味觉感知、核糖体大亚基组装、刺激响应、蛋白质聚合、血小板激活、氧化还原过程、神经嵴细胞的迁移、白细胞介素-2生成的负调控、细胞生长的负调控、细胞死亡负调控、血管生成的负调控、铁离子转运、细胞增殖、典型Wnt信号通路和血管生成；KEGG通路为Ⅰ型糖尿病、Ras信号通路、癌症通路、黑色素生成、赖氨酸降解、HTLV-1感染、乙醚脂质代谢、补体和凝血级联反应和细胞周期。

3 讨 论

热应激可以降低乳蛋白的含量和产量在前期研究中已多次报道[3,5,29]。通过设置PFTN组，发现DMI下降只能部分解释热应激诱导的乳蛋白合成的减少[4-5,30]。截至目前，热应激诱导的乳蛋白合成的下降的原因仍不甚清楚，但可能涉及几个生物系统的变化，包括因营养重分配引起的乳蛋白前体物的供应不足[5,7]，由于乳腺组织能量供应不足导致的mTOR信号通路激活的抑制[5,31]，以及乳腺细胞内代谢的整体抑制和炎症反应激活[8]。

成熟的miRNA可以与靶mRNA结合，通过mRNA降解或翻译抑制在转录后水平调控靶基因的表达[32]。单个miRNA甚至可以靶向调节多达数百个在调节细胞稳态和对各种刺激的反应的过程中起主要作用的基因的表达[9,33-34]。Li等[16]在热应激和非热应激奶牛的乳腺中鉴定出483个牛已知miRNA和139个新miRNA。而本研究在HS组和PFTN组奶牛乳腺中鉴定出578个牛已知miRNA和139个新miRNA。其中bta-miR-451在Li等[16]的研究结果及本研究中都显著上调。Guelfi等[35]和Pasquariello等[36]研究表明，bta-miR-451在卵泡发育和卵母细胞成熟过程中具有重要的功能作用。但bta-miR-451在乳腺组织中的功能目前还没有相关报道。本研究对bta-miR-451的16个靶基因使用DAVID进行功能预测，但是所有富集的通路均不显著。然而，值得注意的是，异柠檬酸脱氢酶3催化亚基α(IDH3A)是bta-miR-451的靶基因之一。已知异柠檬酸脱氢酶可以催化异柠檬酸氧化脱羧生成2-氧戊二酸，是柠檬酸循环中的关键酶。May等[37]研究证实，IDH3A基因的功能缺失可抑制柠檬酸循环的进行，并抑制氧化磷酸化。有趣的是，本团队的前期研究发现，热应激可以很大程度上抑制乳腺组织中能量代谢相关通路的基因表达[8]，其中包括：氧化磷酸化、丙酸代谢和柠檬酸循环。Bionaz等[38]认为mTOR-胰岛素信号通路的激活需要能量，而mTOR信号通路的激活在调控乳蛋白合成中起着关键作用。因此，在本研究中，bta-miR-451的上调可能通过参与乳腺能量代谢的调控，从而参与乳蛋白合成下降的调控。

此外，bta-miR-885和bta-miR-184在Li等[16]和本研究中均出现差异表达。其中bta-miR-885表达量在本研究中显著下调，而在Li等[16]的研究结果中却显著上调。人类临床研究表明miR-885在肝癌中起抑癌作用，可调控肝癌细胞的增殖和从上皮细胞向间质细胞的转化[39]。此外，循环系统中miR-885在包括肝癌、肝硬化、乙肝、急性肝损伤或药物毒性等肝脏病理情况中也出现升高[40-41]。Ylioja等[42]的研究发现，奶牛血液中游离脂肪酸含量升高时，初乳中miR-885含量降低，表明肝脏功能与miR-885之间存在潜在的联系。本研究预测出28个bta-miR-885的靶基因。功能富集分析表明，bta-miR-885的靶基因显著富集于转录和DNA复制的调控和糖酵解/糖异生。然而，本团队前期研究表明HS组奶牛乳腺组织中糖酵解/糖异生相对于PFTN组显著下调。因此，bta-miR-885在热应激诱导的乳蛋白下降的过程中的作用，仍需要进一步的研究分析。

除bta-miR-885外，本研究中bta-miR-184在HS组奶牛乳腺组织中相对于PFTN组的表达量变化也与Li等[16]的研究结果不一致，但与Salama等[43]的结果一致，即热应激增加了bta-miR-184在奶牛乳腺组织中的表达量。本研究共预测得到7个bta-miR-184的靶基因，但靶基因的功能富集分析均不显著。Bu等[10]研究分析了奶牛乳腺上皮细胞中的miRNA表达谱，结果发现bta-miR-184是体外培养乳腺细胞中表达丰度最高的miRNA。Phua等[44]认为miR-184的激活可以抑制三阴性乳腺癌细胞系的体外增殖和自我更新，以及延缓原发肿瘤的形成和降低体内转移负担，提示过表达miR-184将导致细胞生长停滞。此外值得注意的是，本研究中另外2个具有癌细胞抑制功能的miRNA(bta-miR-218和bta-miR-493)也在HS组奶牛乳腺组织中相对于PFTN组表达量下调。其中miR-218被认为可以通过抑制细胞增殖和诱导细胞凋亡来抑制癌细胞的增殖[45-48]，miR-493可以抑制细胞生长和细胞迁移能力来抑制癌细胞[45,49-50]。有趣的是，Salama等[43]研究发现，热应激时促凋亡基因B细胞淋巴瘤/白血病-2相关X蛋白(Bax)的表达量出现上调。Tao等[46]也证实了热应激可降低干奶期奶牛乳腺细胞增殖率。综上所述，本研究中HS组奶牛乳腺组织中bta-miR-184、bta-miR-218和bta-miR-493表达量与PFTN组相比上调，可能通过影响乳腺细胞的增殖或更新，干扰乳腺组织中乳蛋白的合成。

本研究中bta-miR-6516的表达量在HS组奶牛乳腺组织中下调，并发现548个bta-miR-6516的靶基因，高于其他15个差异表达的miRNA。这或许表明bta-miR-6516在HS组奶牛乳腺组织中有更广泛的影响。目前，还没有关于bta-miR-6516或miR-6516的相关研究报道。本研究中对bta-miR-6516的靶基因的功能预测显示，其靶基因主要富集于细胞增殖或细胞生长相关生物过程或通路，这似乎与bta-miR-184、bta-miR-218和bta-miR-493的功能相反，但考虑到bta-miR-6516在乳腺组织中的表达量远低于bta-miR-184、bta-miR-218和bta-miR-493，因而热应激对乳腺细胞的增殖的抑制效应或许会强于促进效应。

在16个差异表达的miRNA中，bta-miR-151-3p在HS组和PFTN组奶牛的乳腺组织中均是表达量最高的miRNA。bta-miR-151-3p共预测发现5个潜在的靶基因。Lu等[47]在所有6种肾病亚型儿童肾组织中均发现miR-151-3p的表达量出现下调，提示miR-151-3p的异常表达可能参与了儿童肾脏疾病的发病机制，并与肾小球损伤相关。Mc Nally等[48]证实miR-151-3p的高表达与胆管癌切除术后总生存期的改善相关。然而，对于miR-151-3p在乳腺中的作用，目前尚缺乏研究。

除bta-miR-151-3p外，bta-miR-136和bta-miR-382也是在HS组和PFTN组奶牛乳腺组织中也高表达的miRNA，且在HS组奶牛乳腺组织中与PFTN组相比表达量上调。bta-miR-136的靶基因主要富集于内分泌及其他钙离子重吸收调节因子相关通路。Jin等[51]研究发现，添加低氧剂可显著促进体外培养细胞的凋亡，降低miR-136的表达，而过表达miR-136的细胞则可以缓解缺氧引起的细胞副作用，提示过表达miR-136可保护神经细胞不受缺氧应激的损伤。然而，本研究中HS组奶牛乳腺组织中bta-miR-136表达量出现上调的原因仍不清楚。bta-miR-382的靶基因主要富集于一些信号通路，及与糖异生和胰岛素拮抗相关的通路，这与我们之前的结果是一致的，即热应激抑制了乳腺组织中糖酵解/糖异生、胰高糖素信号通路、胰岛素拮抗和脂肪细胞因子信号通路相关基因的表达[8]。因此，bta-miR-382的上调可能与mRNA水平的调控具有协同作用。

本研究发现，HS组奶牛乳腺组织中bta-miR-3578表达量相对于PFTN组下调，共预测出bta-miR-3578的靶基因17个。然而，靶基因的富集分析的通路或生物过程均不显著，且目前没有关于bta-miR-3578功能的相关报道。

HS组奶牛乳腺组织中bta-miR-365-5p表达量相对于PFTN组上调。功能注释分析表明，bta-miR-365-5p的靶基因主要富集于癌症、弓形体病和肝炎相关的通路或生物过程。Mori等[52]证实，miR-365的表达可能在人胎盘发育和半异体胚胎的免疫保护中发挥重要作用。Chen等[53]的结果表明，miR-365在循环拉伸应变(cyclic tensile strain，CTS)处理的细胞中显著上调，并在CTS介导的软骨形成中发挥重要作用。但bta-miR-365-5p在HS奶牛乳腺中的潜在功能目前仍不清楚。

bta-miR-338的靶基因富集于与哮喘、肺结核、弓形体病等疾病相关的通路或生物过程。Nie等[54]和Zheng等[55]认为miR-338-3p可以靶向肝细胞和红细胞中糖酵解的关键酶——丙酮酸激酶，然后抑制肝癌细胞的Warburg效应。因此，本研究中bta-miR-338表达量的下调似乎提示HS组奶牛乳腺组织中糖酵解被激活。然而，在本团队前期研究发现，与PFTN组相比，HS组奶牛乳腺组织中糖酵解/糖异生通路被抑制，丙酮酸激酶的表达降低[8]。因此，bta-miR-338在热应激诱导的乳蛋白合成下降中的作用仍需进一步研究。此外，靶基因显著富集于破骨细胞分化相关通路的bta-miR-502a，目前有关研究也很有限，其在热应激诱导的乳蛋白合成下降中的作用仍需要进一步研究分析。

本研究中，在HS组奶牛乳腺组织中显著上调的bta-miR-132的靶基因显著富集于吞噬正调控和炎症反应负调控相关通路中。因此，HS组奶牛乳腺组织中bta-miR-132表达量的上调似乎表明HS组奶牛乳腺组织中炎症反应被激活。与此一致的是，本团队前期对HS组奶牛乳腺组织中差异表达mRNA的生物信息学分析表明，热应激诱导了奶牛乳腺组织中炎症反应的发生[8]。Aten等[56]认为miR-132可能在应激反应和焦虑的调节中发挥关键作用。Diao等[57]证实miR-132可以通过靶向细胞因子信号抑制因子5在慢性阻塞性肺疾病中促进人单核细胞样细胞和人支气管上皮细胞的炎症反应。综上所述，本研究结果似乎提示bta-miR-132可能在热应激诱导的乳腺组织的炎症反应中发挥了一定作用。

功能注释分析表明，bta-miR-147的靶基因显著富集于细胞饥饿反应和癌症中蛋白质多糖相关通路。Liu等[58]研究发现，miR-147b在脂多糖(LPS)刺激后的肺泡巨噬细胞中过度表达。Xu等[59]研究发现，miR-147在大鼠牙周炎诱导的炎症反应中表达显著上调，并证实miR-147对巨噬细胞活化的正向调控机制。然而，bta-miR-147在HS组奶牛乳腺组织中的作用无法根据以往的研究进一步推测。

bta-miR-424-3p的靶基因显著富集于与转录调控、细胞反应、嘧啶代谢等相关通路或生物过程。Dastmalchi等[60]证实miR-424-5p作为肿瘤抑制因子可抑制癌细胞增殖。然而，在食管鳞癌细胞中，miR-424被发现可通过影响细胞周期来促进细胞增殖[61]。因此，miR-424对细胞增殖的调控作用可能是多种多样的，且依赖于细胞类型。因而，bta-miR-424-3p对乳腺细胞增殖的影响还不能确定。

综上所述，在16个差异表达的miRNA中，似乎bta-miR-451、bta-miR-184、bta-miR-218、bta-miR-493、bta-miR-382和bta-miR-132更可能通过上述机制参与热应激下乳腺组织内代谢或生理功能的变化，进而影响乳蛋白合成。此外，对下调和上调miRNA的全部靶基因进行整体的功能注释分析，发现被热应激抑制的能量代谢、mTOR和脂肪细胞信号相关的通路，也被上调的miRNA的靶基因所富集。因而表明，miRNA和mRNA在调节热应激诱导的乳腺组织中代谢的变化过程中具有协同效应。

4 结 论

① HS组奶牛乳腺组织中有16个miRNA相对于PFTN组差异表达。

② bta-miR-451可能参与能量代谢变化的调控，bta-miR-184、bta-miR-218和bta-miR-493可能参与细胞增殖的调控，bta-miR-132可能参与炎症反应的调控，bta-miR-382可能与热应激下乳腺mRNA表达的调节有协同作用。

③ miRNA可能参与了热应激条件下乳腺的代谢和生理变化，为进一步验证miRNA在热应激诱导的乳蛋白合成减少过程中的作用提供了科学依据。