基于高通量测序分析牛呼吸系统疾病鼻咽微生物多样性和功能预测研究

张静，张涛，卢昱希，任荣清，黄波，周思旋，王婧，赵春萍，史开志，张雄*

(1.贵州省农业科学院畜牧兽医研究所/畜禽重大疫病监测防治重点实验室，贵州 贵阳 550002；2.贵州省农业科学院畜牧兽医研究所/兽用中草药研究室，贵州 贵阳 550002)

牛鼻咽部是连接外部环境与机体内环境的重要部位，也是牛呼吸系统疾病(bovine respiratory disease, BRD)致病微生物入侵机体的切入点。据调查，我国养牛业中65%的疾病为BRD，其中肉牛场的感染率达90%以上，死亡率达35%[1]。而这些被诊断患有BRD的牛主要以抗生素治疗为主，大量抗生素的使用产生耐药菌株，最终导致无药可治。因此，系统性地了解患BRD牛鼻咽部微生物群落特征及致病性病原体有助于定向治疗呼吸系统疾病，减少抗菌药物滥用。然而，牛鼻咽部是一个半开放的环境，过境微生物种类繁多，且存在许多不可离体培养的条件性致病微生物，利用传统的培养技术和常规PCR检测来探索牛上呼吸道微生物多样性、分析其功能具有一定局限性；细菌16S rRNA、真菌18S rRNA和真菌核糖体基因转录间隔区(internal transcribed spacer，ITS)测序技术的发展使基于培养方法难以或不可能检测到的微生物物种的鉴定和定量成为可能[2-4]。

近年来，国内外利用微生物扩增子测序对牛消化系统微生物组成和性质，以及微生物群落结构的变化如何影响宿主的健康进行了广泛研究[2，5]。与消化系统微生物群在胃肠道中的作用相似，呼吸道微生物群也是影响呼吸系统健康的重要因素[6]。健康的牛呼吸道菌群处于动态平衡，主要由变形杆菌门(Proteobacteria)、厚壁菌门(Firmicutes)、软壁菌门(Tenericutes)、放线菌门(Actinobacteria)和拟杆菌门(Bacteroidetes)5个门水平下的500余种细菌种属组成，每种细菌种属的相对含量较低[7]。牛鼻咽微生物多样性及每类细菌的相对丰度随牛群所处的生产系统、年龄及健康状况的变化而变化[8]。Holman等[6]采用454焦磷酸测序技术评估了引种第1天和60 d后健康的牛上呼吸道微生物多样性，刚进入饲养场的牛群鼻咽微生物群主要由假单胞菌属(Pseudomonas)、希瓦氏菌属(Shewanella)、不动杆菌属(Acinetobacter)和肉杆菌属(Carnobacterium)等构成，易发生呼吸系统疾病；饲养60 d后牛鼻咽部微生物则逐渐趋于一致，主要以葡萄球菌属(Staphylococcus)、分支杆菌属(Mycobacterium)、曼氏杆菌属(Mannheimia)和莫拉菌属(Moraxella)为主，2个时间点牛群鼻咽微生物群之间存在显著差异。研究还发现，乳杆菌(Lactobacillus)/乳球菌(Lactococcus)丰度增加表明宿主处于健康状态，而溶血性曼氏杆菌(M.haemolytica)、支原体(Mycoplasmabovis)及多杀性巴氏杆菌(Pasteurellamultocida)等相对丰度增加则表示宿主处于疾病状态[9-10]。Zeineldin等[11]也发现，门水平下，BRD牛鼻咽部变形杆菌门(32.12%)、放线菌门(38.2%)和梭杆菌门(3.86%)比健康牛多；属水平下，BRD牛不动杆菌属(Acinetobacter，12.54%)、梭菌属(Fusobacterium，3.71%)和巴氏杆菌属(Pasteurella，2.38%)含量比健康牛高，支原体和莫拉菌属的含量显著高于健康牛。在BRD致死的牛肺灌洗液和肺组织中，致病性病原菌，如溶血性曼氏杆菌(M.haemolytica)、支原体及昏睡嗜血杆菌(H.somnus)的相对丰度都大于5%，其他种属在1%～5%之间，而健康的牛肺灌洗液中未发现支原体[3，12]。由此可知，鼻咽微生物的组成及丰度能直接影响牛下呼吸道中的微生态平衡，导致牛呼吸道疾病加剧。部分研究还表明真菌在慢性呼吸道疾病的发生和发展过程中也具有重要作用，且易产生复杂的耐药性[13]。张传师等[14]总结了国外的相关研究主要集中在鼻咽部细菌类微生物组与BRD的关系上，缺乏对牛呼吸道真菌微生物组学的研究，在国内则未见对患呼吸系统疾病牛鼻咽部微生物群的相关研究报道。因此，从鼻咽部细菌组和真菌组群落组成出发，全面了解并阐述微生物组与BRD之间的关系，有助于进一步探讨疾病发病机制，辅助临床防治工作。

本研究以2019—2020年期间采集的贵州省内3个牛场患BRD的牛鼻咽微生物为研究对象，利用16S rRNA V3-V4区和ITS1区Illumina NovaSeq测序，横向比较不同地区患病牛鼻咽部微生物群落组成并预测其功能，以期为系统了解养牛场发病病因，深入研究鼻咽微生物的组成与宿主呼吸系统健康之间的关系，开发新型微生态制剂及指导牛呼吸系统疾病防制工作提供基础数据。

1 材料与方法

1.1 样品采集

R1组共14份样品，2020年4月采自贵州省盘县某牛场，该牛场患病牛临床表现为中等以上呼吸系统疾病症状，具体表现为咳嗽、鼻腔分泌物多，严重时鼻腔有出血现象，持续时间10～20 d；R2组共6份样品，2020年5月采自贵州省盘县某牛场，该牛场患病牛临床表现为咳嗽，鼻腔分泌物增加，食欲减退，精神不佳的症状，持续时间为5～15 d；R3组共12份样品，2019年12月采自贵州省大方县某牛场，该牛场患病牛临床表现为流鼻涕、轻微咳嗽、食欲减退，持续时间为10～15 d。用无菌棉签采集患病牛鼻腔拭子，于无菌管中保存，放置液氮罐中送回实验室于-80 ℃保存备用。

1.2 样品中DNA的提取及PCR扩增

利用十六烷基三甲基溴化铵(CTAB)法对鼻拭子样品基因组DNA进行提取。之后利用琼脂糖凝胶电泳检测DNA的纯度和浓度，取适量的样本DNA于离心管中，使用无菌水稀释样本至1 ng/μL。以稀释后的基因组DNA为模板，扩增细菌16S V3-V4高变区(341F：5′-CCTAYGGGRBGCASCAG-3′， 806R：5′-GGACTACNNGGGTATCTAAT-3′)和真菌ITS1区(ITS1-1F-F：5′-CTTGGTCATTTAGAGGAAGTAA-3′，ITS1-1F-R：5′-GCTGCGTTCTTCATCGATGC-3′)，文库构建及Ilumina NovaSeq测序由北京诺禾致源科技股份有限公司协助完成。

1.3 测序数据处理

将每个样本测序获得的读段(reads)进行拼接，得到的拼接序列为原始标签(raw tags)数据。将原始标签序列经过严格的过滤处理得到高质量的标签数据(clean tags)。利用QIIME[15]软件对高质量的标签数据进行质量控制，具体操作如下：①标签截取：将原始标签从连续3个以上碱基质量小于或等于19的第1个低质量碱基位点截断；②标签长度过滤：标签数据经过截取后得到的标签数据集，进一步过滤掉其中连续高质量碱基长度小于标签长度75%的标签序列。质控后的标签序列去除其中的嵌合体序列，得到最终的有效数据(effective tags)进行数据分析。

1.4 生物信息学分析

利用Uparse[16]对所有样本的有效序列以97%的一致性聚类成为分类操作单元(operational taxonomic units，OTUs)，为检验测序数据是否达到饱和且符合后续分析要求，从样本中随机抽取一定测序量的数据，统计它们所代表OTUs数目，以抽取的测序数据量与对应的物种数来构建稀释曲线(rarefraction curve)，并根据样本中OTUs的相对丰度绘制等级聚类曲线(rank abundance curve)。使用QIIME软件进行样本复杂度分析(alpha diversity)，包括Chao1，ACE，Shannon，Simpson和Observed species共5个指数；计算Unifrac距离，进行多组间比较分析(beta diversity)并绘制PCoA图。对OTUs序列进行物种注释，用Mothur[17]和QIIME分别对细菌和真菌进行物种注释，在门和属分类水平上，统计不同组中的物种丰度和多样性，为进一步挖掘组间微生物菌落的差异，选用t检验对3组样本两两之间的微生物菌落组成进行统计分析，筛选出各组存在显著差异的菌落。

1.5 16S rRNA和ITS功能预测分析

为探讨牛鼻咽部微生物群落在其体内发挥的功能，使用PICRUSt[18]和FunGuild[19]方法分别对细菌和真菌进行了功能预测，利用R软件绘制每组样本的丰度信息热图，并从功能差异层面进行聚类。

1.6 数据的统计与分析

Alpha多样性指数组间差异分析选用Tukey检验，结果表示为“平均数±标准差”，P<0.05表示差异显著，P<0.01表示差异极显著。

2 结果与分析

2.1 序列拼接质控与OTUs聚类分析

经16S rRNA和ITS基因高通量测序，平均测得牛上呼吸道内细菌原始序列为87 696条，经质控后平均得到64 128条有效序列；R1、R2和R3组对应的平均有效序列分别为64 880、63 545和63 959条。平均获得真菌原始序列105 973条，经质控后，平均得到有效序列65 472条；R1、R2和R3组对应的平均有效序列分别为64 285、65 560和66 572条(表1)。在97%相似度标准下聚类分析，分别得到细菌OTU为2 851、2 358和2 728个，其中3组样本共有的OTU数目为1 683，R1、R2和R3组特有的OTU数目分别为474、225和353；得到的真菌OTU分别为1 766、1 213和2 125个，其中3组样本共有的OTU数目为825个，R1、R2和R3组特有的OTU数目分别为263、139和607个(图1)。

表1 牛鼻咽部细菌16S rRNA基因和真菌ITS高通量测序基本信息表

图1 3组牛鼻咽部细菌16S rRNA(A)和真菌ITS(B)聚类OTUs韦恩图

2.2 微生物群落的Alpha多样性分析

16S rRNA和ITS测序数据的稀释曲线和等级聚类曲线结果表明，各组样本量足够，测序深度已经基本覆盖到样品中所有物种(图2、图3)。Alpha多样性指数反映样品微生物群落的多样性。其中，Chao1、ACE 指数反映菌群丰度，由表2可知，R1、R2和R3组之间上呼吸道细菌及真菌菌群丰度相似，R2组细菌和真菌的菌群丰度都最高，细菌菌群丰度最低的是R3组，真菌菌群丰度最低的是R1组。菌群的多样性由Simpson和Shannon指数反映，数据显示，细菌群落多样性最高的是R2组，R1组次之，R3组最低；R1和R2组菌群多样性显著(P<0.05)高于R3组，R1组与R2组之间无显著差异；真菌菌落多样性最高的是R3组，其次是R2组，最低的是R1组；R1组菌群多样性显著(P<0.05)低于R3组，R1组与R2组之间无显著差异。表明3组样品中菌群丰度虽然相似，但菌落多样性存在较大差异。

图2 不同组别中的16S rRNA稀释曲线(A)和等级聚类曲线(B)

图3 不同组别中的ITS稀释曲线(A)和等级聚类曲线(B)

表2 3组样本的Alpha多样性指数

2.3 微生物群落的Beta多样性分析

基于加权的Unifrac距离进行PCoA分析，分析3组样品间微生物的相似性。由图4A可知，细菌类主坐标PC1和PC2的贡献率分别为总变量的38.43%和17.81%，3组样品并没有各自单独聚类，说明9组样品细菌微生物群落组成相似。而真菌类微生物PCoA分析结果显示(图4B)，主坐标PC1和PC2的贡献率分别为总变量的45.69%和23.87%，R1和R2组聚为一类，R3组则单独聚类，表明R3组中的多数患病牛上呼吸道真菌群落组成与R1、R2组差异较大。

图4 不同地区患病牛上呼吸道中细菌类微生物群落PCoA图(A)和真菌类微生物群落PCoA图(B)

2.4 不同地区患病牛上呼吸道中的微生物品种鉴定及菌群结构分析

3个不同牛场患病牛上呼吸道中细菌类微生物门水平分布如表3所示，主要以变形杆菌门(Proteobacteria)、拟杆菌门(Bacteroidetes)、梭杆菌门(Fusobacteria)、厚壁菌门(Firmicutes)、软壁菌门(Tenericutes)及放线菌门(Actinobacteria)为主；与R1组相比，R2组变形杆菌门显著(P<0.05)减少，R3组拟杆菌门、厚壁菌门及放线菌门显著减少；与R2组相比，R3组变形杆菌门显著增加，拟杆菌门则显著减少。属水平结果显示，3组样品中排前10的主要菌群属组成相似，但相对丰度相差较大。R1组中主要以拟杆菌属(Bacteroides)、不动杆菌属、卟啉菌属(Porphyromonas)及梭菌属(Fusobacterium)，其中不动杆菌属和拟杆菌属显著高于R3组；R2组中主要以卟啉菌属、拟杆菌属、梭菌属及不动杆菌属，其中不动杆菌属显著高于R3组，R3组中主要以莫拉菌属(Moraxella)、支原体属(Mycoplasma)、梭菌属及卟啉菌属，其中莫拉菌属显著高于R1和R2组。

表3 细菌分类门水平和排前10的属水平的相对分布

真菌类微生物门水平分布如表4所示，3组样品中，主要以子囊菌门(Ascomycota)和担子菌门(Basidiomycota)为主，R1组和R2组内的担子菌门极显著高于(P<0.01)R3组。属水平上，3组样品排前10的主要菌群属丰度相差较大。R1组主要以克柔氏菌属(Issatchenkia)、隐毛孢耳属(Apiotrichum)、假丝酵母属(Diutina)、毛孢子菌属(Trichosporon)和曲霉菌属(Aspergillus)为主；R2组主要以克柔氏菌属、隐毛孢耳属和毛孢子菌属为主；R3组主要以曲霉菌属、(Millerozyma)和克柔菌属为主。与R1组相比，R2组假丝酵母属相对丰度极显著减少；与R3组相比，R1组和R2组克柔氏菌属丰度增加，曲霉菌属丰度减少，差异极显著。

表4 真菌分类门水平和排前10的属水平的相对分布

2.5 上呼吸道微生物功能预测

用PICRUSt对不同地区患病牛上呼吸道细菌类群进行功能预测。由图5可知，KEGG通路二级功能层面上丰度排名前35的预测功能基因及它们在每个组中的丰度信息。R1组样品有18个功能基因通路丰度减少，17个功能基因通路丰度增加；R2组样品有12个功能基因通路丰度减少，23个功能基因通路丰度增加；R3组样品有21个功能基因丰度减少，14个功能基因丰度增加。R1组和R2组功能基因丰度富集相似；然而，R3组中功能基因丰度增加的通路在R1和R2组中的丰度降低。其中，R3组中丰度增加的功能基因通路有：疾病感染、遗传信息呈递、转录、细胞过程和信号、信号转导、癌症和神经退行性疾病等相关通路；丰度减少的功能基因通路有：核苷酸代谢、细胞生长和死亡、消化系统、代谢疾病、氨基酸代谢和碳水化合物代谢等相关通路。结果表明，R1组和R2组患病牛感染的病原能够导致代谢类的相关通路丰度增加，其中R2组丰度增加的通路比R1组多；疾病感染、信号转导等相关通路丰度减少。而R3组患病牛感染的病原功能基因丰度富集的结果与R2组相反。

图5 不同地区患病牛上呼吸道中细菌类微生物功能注释聚类热图

基于真菌的物种分类，FunGuild可以获得相应真菌的生态功能。用FunGuild对不同地区患病牛上呼吸道真菌类群进行功能预测。由图6可知，3组样品的上呼吸道中均存在动物病原菌富集及丰度增加；其中R1和R2组中的动物病原菌分别有2种，R1组种有1种，3组样品中的病原菌种类各不相同。

图6 不同地区患病牛上呼吸道中真菌类微生物功能注释聚类热图

3 讨论

3.1 呼吸系统疾病牛鼻咽部菌群多样性

本研究通过高通量测序技术对患BRD的牛鼻咽部菌群结构及组成进行分析，并基于细菌和真菌注释数据库进行菌群功能预测。不同场中微生物Alpha多样性分析和组间微生物群落Beta多样性分析结果表明，不同发病牛场细菌和真菌丰度相似；同一地区不同牛场菌群多样性及菌落结构组成相似，但不同地区牛场之间菌群多样性存在显著差异，菌落结构组成也有所不同。由此可知，牛慢性呼吸系统疾病受地区环境因素影响较大，不同地区导致发病的诱因存在差异。临床上，需进一步采用分子生物学技术或微生物菌属注释确认病原后进行针对性的治疗，否则会导致牛场抗生素滥用，产生耐药性。除此之外，随着生物安全控制体系在畜禽疫病防控上获得的显著成效，牛场生物安全控制体系的建立也应该被重视，使其成为控制和预防地区性牛呼吸系统疾病蔓延的关键。

本研究结果显示，患病牛呼吸道细菌类微生物主要由6个门类组成，与Zeineldin等[11]研究的BRD犊牛鼻咽优势菌群(变形杆菌门32.12%、放线菌门38.2%和梭杆菌门3.86%)组成相比，本研究的样品中还存在拟杆菌门(19.19%～47.15%)、厚壁菌门(9.59%～16.65%)和软壁菌门(1.06%～7.93%)；而放线菌门相对丰度较低，梭杆菌门相对丰度较高，可能是由不同地区及不同环境所致。与正常肉牛鼻咽部微生物组成[20]相比，患BRD牛能检测到拟杆菌门和梭杆菌门的存在，放线杆菌门和软壁菌门的相对丰度增加，厚壁菌门和变形杆菌门相对丰度类似，说明BRD可导致部分鼻咽部菌群组成及相对丰度变化。属水平上，患BRD牛鼻咽部微生物相对丰度均小于20%，其中包括莫拉菌属(1.04%～14.99%)、卟啉菌属(5.16%～19.35%)、拟杆菌属(2.12%～16%)、不动杆菌属(0.74%～10.92%)、支原体属(0.68%～7.82%)和曼氏杆菌属(0.19%～3.06%)。与大多数研究结果[8]一致，支原体属、莫拉菌属、曼氏杆菌属、巴氏杆菌属在不同地区的患BRD牛鼻咽部中均存在，推测这3种菌属可能在牛呼吸系统疾病中起到关键作用，多重感染的具体作用机制仍需进一步研究。

研究表明，真菌在人类呼吸系统疾病的发生、发展过程中扮演着重要角色，如真菌与喘息病、耐药性等相关，也参与了宿主黏膜免疫反应，并证实了呼吸系统疾病的发生会导致真菌群的变化[21]。本研究对患BRD牛鼻咽部真菌类微生物进行了探究，发现与慢性呼吸系统疾病相关的优势菌属，如克柔氏菌属(2.28%～45.1%)、曲霉菌属(1.19%～16.26%)等，表明患BRD的牛呼吸系统会受到致病性真菌的入侵，最终导致真菌性感染。

3.2 呼吸系统疾病牛鼻咽部菌群中的致病菌

Lima等[9]利用高通量测序技术对健康和患BRD牛的鼻咽微生物进行了分析，与健康牛相比，患病牛鼻咽部中的微生物多样性降低，致病性病原菌如：牛支原体、溶血性曼氏杆菌、牛巴氏杆菌等相对丰度显著升高。本试验在不同地区患病牛鼻咽部中也检测到以上3种致病性细菌，除此之外，还检测到不动杆菌属和金黄杆菌属。研究发现，金黄杆菌感染常发生在慢性呼吸道疾病和免疫力低下的宿主中[22]；弓清梅等[23]发现脑膜脓毒金黄杆菌会导致下呼吸道疾病的发生，且对绝大多数抗生素耐药。产吲哚金黄杆菌的耐药性也极其复杂[24]。不动杆菌是一种机会性致病菌，常引起肺部炎症，其易感群体与金黄杆菌相似，该菌对市面上常用的抗菌素具有耐药性，导致临床治疗困难，病死率较高[25]。以上结果表明，除常见的致病性细菌外，金黄杆菌属和不动杆菌属可能导致患慢性呼吸道疾病牛和长期使用抗菌素的牛场的继发性感染。

本试验对患BRD牛鼻咽真菌菌群进行了测序分析，发现与呼吸系统疾病相关的真菌属，如：曲霉菌属、假丝酵母属和克柔氏菌属等，其中，曲霉菌属在慢性呼吸系统疾病中最为常见[13]。本研究3组样品中均检测曲霉菌属，R3组中相对丰度最高，表明曲霉菌属是R3组的主要致病性真菌。假丝酵母属也是引起下呼吸道感染的真菌之一，在R1组中的相对丰度大于5%。陈太方等[26]在研究抗真菌药物对下呼吸道感染的抑菌作用中发现假丝酵母属检出率是多种真菌中最高的，且具有一定的耐药性。克柔氏菌属在R1和R2组中的相对丰度均大于30%，该菌属条件性致病菌，高丰度的克柔氏菌可导致支气管肺炎的发生[27]。上述结果表明了在患呼吸系统疾病牛鼻咽微生物可检测出致病性真菌，可能影响呼吸道疾病的发展，具体致病机制仍需进一步探究。

3.3 呼吸系统疾病牛鼻咽部菌群功能预测

鼻咽部微生物作为牛机体微生态系统的重要组成部分，微生物组与机体的互作对维持宿主生理稳态和机体健康至关重要。在相似的的环境中，尽管群落结果发生了巨大变化，但其功能仍然相对稳定[28]。研究者对牛瘤胃微生物及肠道微生物功能进行了分析，发现功能集中在氨基酸运输和代谢、碳水化合物转运及代谢的相关基因上[2]。大多数的牛呼吸道疾病微生物研究只分析了结构组成、相对丰度及其与疾病的相关性，很少针对微生物测序结果进行功能预测[8]。基于PICRUSt和FunGuild数据库，本研究揭示了在患呼吸系统疾病的牛鼻咽部细菌类KEGG二级水平的通路中，R1组和R2组细菌类微生物导致代谢类的相关通路(氨基酸代谢、碳水化合物代谢、消化系统等)丰度增加，可能是由于细胞和结缔组织破坏，与肺功能受损和全身炎症相关[29]。其中消化系统相对丰度增加，可能与肺-肠轴通路有关[30]。R3组与其他2组的富集结果存在差异，功能基因丰度增加的通路主要包括疾病感染、遗传信息呈递和信号转导等，这可能与呼吸系统疾病的发生、微生物在呼吸道中增殖及宿主-微生物互作有关。真菌类群功能预测结果显示，3组样品的上呼吸道中均存在动物病原菌富集及丰度增加，表明牛呼吸系统疾病存在致病性真菌感染的可能，为临床治疗提供了方向。

4 结论

本研究从患BRD牛鼻咽微生物出发，分析了鼻咽微生物中细菌、真菌的群落组成、相对丰度及功能预测。鉴定出牛鼻咽部致病性细菌，如：牛支原体、巴氏杆菌、曼氏杆菌、金黄杆菌及不动杆菌；致病性真菌，如：曲霉菌属、假丝酵母属和克柔菌属，说明患呼吸系统疾病牛鼻咽微生物菌落组成丰富，且包含多种致病性微生物。