基于16S rDNA扩增子测序技术分析四川省奶牛乳房炎乳汁中微生物菌落多样性

肖璐，谢晶，于吉锋，廖党金，曹冶，林毅，李兴玉，潘梦，叶健强，叶勇刚，康润敏

(四川省畜牧科学研究院动物遗传育种四川省重点实验室，四川 成都 610066)

奶牛乳房炎是影响全球奶业发展的主要疾病之一，奶牛乳房炎是由微生物感染、物理或化学刺激等因素引起的奶牛乳房乳腺组织的炎症反应。据报道与统计，全世界约有1/3的奶牛患有各种类型的乳房炎，有30%的奶牛在泌乳期间至少患1次临床型乳房炎[1-2]。奶牛乳房炎可分为临床型乳房炎和亚临床型乳房炎(又称隐性乳房炎)，该病的致病因素复杂、发病率高、治愈反复发病率高，由于病原微生物的入侵，影响奶牛乳腺健康，导致产奶量下降、乳品质降低，困扰着奶牛业的经济效益与健康发展，给全球的奶牛产业带来了巨大的经济损失[3]。引发奶牛乳房炎的病原微生物包括：球菌、杆菌、支原体、真菌、酵母和病毒等，根据来源和传播途径可以分为：传染型致病原和环境致病原，传染型致病原主要包括无乳链球菌、停乳链球菌、金黄色葡萄球菌和支原体等，环境型致病原主要包括大肠杆菌、乳房链球菌、克雷伯氏菌、绿脓杆菌、产气肠杆菌、变形杆菌、凝固酶阴性葡萄球菌、牛棒状杆菌、真菌及酵母菌等[4]。研究表明导致临床中奶牛乳房炎发生的致病菌以金黄色葡萄球菌、无乳链球菌、停乳链球菌和大肠杆菌为主，占到引起奶牛乳房炎发病率的80%～90%[5]，但不同地方导致奶牛乳房炎主要致病菌流行的趋势也存在显著性差异，且耐药性也不尽相同[6-7]。在牛奶乳房炎的防治中，养殖户为了最大程度地减少乳房炎带来的经济损失，滥用和乱用抗生素十分严重，导致牛奶中残留大量的抗生素，给食品安全和人体健康均造成了危害[8-9]，且已有报道指出奶牛乳房炎的防治已经成为奶牛养殖业中使用抗生素的主要原因[10-11]。乳房炎的发生与多种病原微生物相关，但乳房菌群之间存在某种共生、互生、竞争、寄生等关系，可通过宏基因组测序分析，直接获取乳汁中的菌群信息。本研究利用16S rDNA扩增子测序技术对四川省成都市、眉山市和凉山州3个地区的10家养殖场的奶牛乳房炎乳汁和正常乳汁微生物菌群多样性进行比较分析，分析两者直接的菌群结构差异，掌握四川省奶牛乳房炎流行主要致病菌，为针对性地使用抗生素、保障食品安全和人民健康提供科学依据。

1 材料与方法

1.1 样品采集

2018年6月分别从眉山市洪雅县、成都市青白江区、凉山州西昌市某规模化奶牛场随机选择健康奶牛(体细胞数<20万/mL)共计20头作为对照组；临床型乳房炎发病奶牛(体细胞数>100万/mL)共计15头作为试验组，样品信息详见表1。对奶牛乳房进行消毒后，以无菌EP管采集对照组和试验组的各奶牛的新鲜乳汁，随即将采集后的乳汁编号放入冰盒中低温保存送往实验室，其中正常奶牛乳汁编号为ZCY1～ZCY20，临床乳房炎奶牛乳汁样编号为RFY1～RFY15。

1.2 主要试剂

血液/细胞/组织基因组DNA提取试剂盒购自天根生化科技(北京)有限公司，Phusion®High-Fidelity PCR Master Mix with GC Buffer、NEB Nextr®Ultra TMDNA Library Prep Kit for Illumina 建库试剂盒、高效高保真酶均购自New England Biolabs，普通琼脂糖DNA回收试剂盒购自天根生化科技(北京)有限公司。

表1 样品信息

1.3 乳汁中微生物总DNA提取

从无菌的EP管中取乳汁样本200 μL按照血液/细胞/组织基因组DNA提取试剂盒的说明书提取乳汁中微生物总DNA，使用紫外分光光度仪(NanoDrop 2000)和琼脂糖凝胶电泳判断DNA样品的浓度与纯度。

1.4 乳汁中微生物16S rDNA的V3-V4高变区扩增

选择16S rDNA中的V3-V4区域作为目的片段进行扩增，引物序列为515F：5′-CCTAYGGGRBGCASCAG-3′，806R：5′-GGACTACNNGGGTATCTAAT-3′，引物由生工生物工程(上海)股份有限公司合成，并在上下游引物上分别添加不同的barcode。PCR反应体系为30 μL：Phusion Master Mix(2×)15 μL，上下游引物(2 μmol/L) 3 μL，gDNA(1 ng/μL) 10 μL，H2O 2 μL；反应程序：98 ℃ 1 min；98 ℃ 10 s，50 ℃ 30 s，72 ℃ 30 s，30个循环；72 ℃ 5 min。PCR产物经2%琼脂糖凝胶电泳检测其浓度和纯度，选择大小为400～450 bp的条带，并回收纯化目的片段送公司测序分析。

1.5 DNA文库构建和上机测序

使用NEBNext®UltraTMDNA Library Prep Kit for Illumina建库试剂盒进行DNA文库建立，将构建好的文库经过Qubit定量和文库检测合格后，使用HiSeq进行上机测序，测序委托北京诺禾致源生物信息科技有限公司完成。

1.6 测序信息分析

①对下机获得的原始数据(raw tags)进行拼接、过滤得到高质量数据(clean tags)，再经过Qiime的Tags质量控制流程和Tag去除嵌合体序列处理，最终获得有效数据(effective tags)。基于有效数据，利用Uprse v7.0.1001对有效数据进行OTUs(operational taxonomic units)聚类和物种分类分析，筛选出代表性OTUs序列，利用Mothur方法与SILVA的SSUrRNA数据库进行物种注释分析，获得不同分类水平上物种信息和基于物种的丰度分布情况；②使用Qiime软件(Version 1.9.1)计算Alpha多样性指数，获得35个样本菌群丰度和多样性；③利用R软件(Version 2.15.3)绘制稀释曲线，确定测序深度，选用wilcox检验对2个组Alpha多样性指数组间进行差异分析。使用R vegan包的anosim函数进行Anosim分析，验证2组组间菌落结构的差异是否显著大于组内差异，从而判断分组是否有意义；④基于R软件进行Metastats分析和Benjamini and Hochberg False Discovery Rate方法对各分类水平2组间物种丰度数据进行分析，筛选出2个组间存在显著差异的菌群。分析结果用均用“平均值±标准差”来表示。

2 结果

2.1 乳房炎和正常奶牛乳汁样本测序和OTU注释

通过Illumina Highseq测序后35个样本共获得2 839 426个原始数据(表2)，经过拼接、过滤、质控和去除嵌合体后有1 697 422个有效数据进入后期数据分析，其中乳房炎乳汁组和正常乳汁组每个样本分别获得(51 438±6 613)和(44 577±4 058)个有效数据。利用Uparse软件对35个样本中所有的有效数据进行聚类，以97%的一致性为原则将序列聚类成为OTUs。结果显示：试验组和对照组分别有(46 480±7 381)和(44 158±4 072)个有效数据被聚类到(1 509±488)和(249±243)个OTUs中。从表3可见，2个组有效数据注释水平随着物种分类等级的下降不断减少，其中试验组和对照组分别由61%和49%的有效数据能注释到属水平，而仅仅只有33%和42%的有效数据被注释到种水平上。

表2 35个样本测序数据和OTU聚类

表3 35个样本注释后在各分类的序列数目

2.2 乳房炎和正常奶牛乳汁样本的多样性分析

通过样品的多样性分析，可以反映样品内的细菌群落的丰富度和多样性，由稀释曲线可知当测序序列条数达到30 000时，稀释曲线趋于平滑衍生，OUT的数目接近饱和，表明现在所获得的数据已经足够反应样本菌群的多样性情况，且试验组和对照组35个样本的覆盖率(goods coverage)值均在0.99以上，两者共同说明本次所测序分析的35个样品微生物群落的检测比率接近饱和，目前的测序量能够覆盖样本中绝大部分物种。使用Qiime软件计算获得35个样本的Alpha多样性指数，通过2组样品的观测物种数(Observed species)、香农指数(Shannon index)、辛普森多样性指数(Simpson index)、Chao1指数(Chao1)，文库覆盖度(Good’ coverage)见表4，经过wilcox检验结果显示：试验组的观测物种数值显著高于对照组(P<0.05)，两者相差6～7倍，说明患乳房炎组奶牛乳汁中含有的菌落种类的数目远比健康奶牛乳汁组更加丰富；试验组香农和辛普森指数均显著高于对照组(P<0.05)，说明2组菌落分布多样性存在显著差异；而从Chao1和ACE指数来看，试验组也显著高于对照组(P<0.05)，说明2组中菌群的丰度也存在显著差异。

表4 35个样本Alpha多样性指数分析

2.3 乳房炎和正常奶牛乳汁样本分组合理性分析

2组菌落结构差异与组间结构差异的相关性分析结果表明：2组样品的相关系数(R值)为0.939(图1)，该结果证明组间菌落结构的差异显著大于组内菌落结构的差异(P<0.001)，由此可知分组合理。

图1 Anosim分析2组间菌落结构差异分析

2.4 乳房炎和正常奶牛乳汁组间微生物菌群多样性差异分析

利用Uprse v7.0.1001对有效数据进行OTUs聚类和物种分类分析表明：在门水平上，试验组和对照组样品中最为丰富的菌群均为：变形菌门(Proteobacteria)、厚壁菌门(Firmicutes)、拟杆菌门(Bacteroidetes)和放线菌门(Actinobacteria)，试验组分别占样品序列总数比例依次是：33%、31%、7.6%、6.8%，对照组占样品序列总数的比例依次为： 61%、35%、1.9%、0.2%(图2)。在属水平上，链球菌属(Streptococcus)、叶杆菌属(Phyllobacterium)、假单胞菌属 (Pseudomonas)为试验组最为丰富的微生物菌群，分别占样品序列的7.6%、6.97%、4.25%，而乳球菌属 (Lactococcus) 、不动菌属 (Acinetobacter)、假单胞菌属 (Pseudomonas)则是正常乳汁组(对照组)最为丰富的微生物菌群，分别占样品序列的32.3%、5.0%、1.62%(图3)。

图2 门水平上对照组和试验组乳汁的菌落结构

图3 属水平上乳汁的菌落结构

Metastats分析和Benjamini and Hochberg False Discovery Rate方法对各分类水平2组间物种丰度数据的分析结果显示：2组间在门水平上有27个门的相对丰度存在显著差异，其中除放线菌门 (Actinobacteria) 对照组的相对丰度显著高于试验组外(P<0.05)，其余26个菌门，包括变形菌门(Proteob-acteria)、厚壁菌门(Bacteroidetes)、绿弯菌门(Chlo-roflexi)、酸杆菌门(Acidobacteria)、栖热链球菌门(Deinococcus-Thermus)、广古菌门(Euryarchaeota)、芽单胞菌门(Gemmatimonadetes)、螺旋体门(Spirochaetes)、软壁菌门(Tenericutes)、产氧光菌门(Oxyphotobacteria)、疣微菌门(Verrucomicrobia)、黑水仙菌门(Melainabacteria)、梭杆菌门(Fusobacteria)、纤维杆菌门(Fibrobacteres)、Candidatus Yanofskybacteria、Kiritimatiellaeota、互养菌门(Syne-rgistetes)、迷踪菌门(Elusimicrobia)、Calditrichaeota、Acetothermia、奇古菌门(Thaumarchaeota)、黏胶球形菌门(Lentisphaerae)、衣原体门(Chlamydiae)、泉古菌门(Crenarchaeota)、Diapherotrites、Margulisbacteria、装甲菌门(Armatimonadetes)，其在试验组的相对丰度均显著高于对照组(P<0.05)。

在属水平上，2组有350个属的相对丰度存在显著性差异，其中有113个属属于变形菌门(Proteobacteria)，105个属属于厚壁菌门(Firmicutes)，48个属属于放线菌门(Actinobacteria)，45个属于拟杆菌门(Bacteroidetes)，为2组相对分析显著差异属主要聚集的门。从2组中相对丰度差异最显著的前20个属中可见，对照组乳酸链球菌属(Lactococcus)、乳杆菌属(Lactobacillus)、果胶杆菌属(Pectobacterium)、柠檬酸杆菌属(Citrobacter)和亚特兰大杆菌属(Atlantibacter)的相对丰度显著性的高于试验组(P<0.05)，而试验组的链球菌属(Streptococcus)、隐秘杆菌属(Trueperella)、水栖菌属(Enhydrobact-er)、葡萄球菌属(Staphylococcus)、肉杆菌属(Atopostipes)、短波单胞菌属(Brevundimonas)、气球菌属 (Aerococcus)、Guggenheimella属、嗜盐单胞菌属(Halomonas)、魏斯氏菌属(Weissella)、泰氏菌属(Tissierella)、拟杆菌属(Bacteroides)、不活动粒菌属(Ignavigranum)、费克蓝姆氏菌属(Facklamia)和Iamia属的相对丰度显著高于对照组(P<0.05)。

在种水平上，2组有213个种的相对丰度存在显著性差异(P<0.05)。从2组相对丰度差异最显著的前20个种中可见，对照组乳酸乳球菌(Lactococcuslactis)、乳酸杆菌(Lactobacillushelveticus)、约氏不动杆菌(Acinetobacterjohnsonii)、绿色气球菌(Aeroeoccusviridans)和蜂房哈夫尼菌(Hafniaalvei)的相对丰度显著性高于试验组，而试验组的无乳链球菌(Streptococcusagalactiae)、nosocomialis不动杆菌(Acinetobacternosocomialis)、化脓隐秘菌(Trueperellapyogenes)、奥斯陆莫拉氏菌(Moraxellaosloensis)、泡囊短波单胞菌(Brevundimonasvesicularis)、食窦魏斯氏菌(Weissella.cibaria)、丙杆菌(Caulobacter)、凝固酶阴性葡萄球菌(Staphylococcuspettenkoferi)、海迪茨氏菌(Dietziamaris)、乳微杆菌(Microbacteriumlacticum)、双孢梭菌 (Clostridiumdisporicum)、结膜干燥棒状杆菌(Corynebacteriumxerosis)、Piceasmithiana菌、塔氏不动杆菌(Acinetobactertowneri)、绿色气球菌(Aerococcusviridans)和戴尔福特菌(Delftiatsuruhatensis)的相对丰度显著性高于对照组(P<0.05)。

3 讨论

奶牛乳房炎是全球公认的奶牛群体中最常见、治疗费用最昂贵的疾病，给全球奶牛养殖业带了巨大的经济损失，而我国的奶牛乳房炎发病情况更加严峻，严重制约了奶牛业的健康发展。奶牛乳房炎的病因非常复杂，病原繁多，有超过150种致病微生物可导致奶牛乳房炎的发生，较为常见的有20多种[11-13]。宏基因组学测序技术的出现，突破了传统微生物技术的局限性，使得很多难培养及含量低的细菌可以被测出，通过对正常与患乳房炎奶牛菌群结构的研究，以期发现奶牛乳房炎与乳汁菌群之间可能存在的因果关系。本研究依托16S rDNA扩增子测序技术，对采自四川省奶牛主要养殖区成都、眉山、凉山州3地的10家奶牛场中的20份健康奶牛乳汁和15份乳房炎病牛乳汁进行了测序，横向比较了患乳房炎病牛与正常奶牛乳汁之间的菌群差异，并分析患乳房炎病牛乳汁中的主要致病菌，探索奶牛乳房炎发生的微生物群落特征差异。

通过对四川省部分奶牛场生产中正常乳汁和乳房炎乳汁的微生物菌群的差异分析，发现2组菌群在门水平上有27个门的相对丰度存在显著差异，在属水平上有350个属的相对丰度有显著差异，在种水平上有213个种的相对丰度存在显著性差异。从门的水平上看，曾学琴等[14]研究表明，乳房炎奶牛的牛奶中丰度排名前3的菌门依次是变形菌门、厚壁菌门、放线菌门；Catozzi等[15]发现，患乳房炎奶牛的牛奶中排名前3的菌门依次为厚壁菌门、变性菌门、拟杆菌门；陶娜拉[16]分析发现，强阳及弱阳性乳房炎奶牛乳汁中的优势菌门为变形菌门、拟杆菌门和厚壁菌门，其结果与本次调查结果不完全一致。四川地区患乳房炎奶牛的牛奶中排名前3的菌门依次为变形菌门、放线菌门及拟杆菌门。从丰度分布来看，在传染型致病原中，四川地区患乳房炎奶牛与正常奶牛乳汁中存在显著差异的致病菌为无乳链球菌、凝固酶阴性葡萄状球、化脓隐秘菌和绿色气球菌，而常见的金黄色葡萄球菌2组并无显著差异；在环境型致病原中2组乳汁中凝固酶阴性葡萄状球菌、化脓隐秘菌、绿色气球菌丰度也出现了显著升高。

葡萄球菌是引起奶牛乳房炎常见、重要病原，其中，金黄色葡萄球菌属于凝固酶阳性菌，溶血性葡萄球菌、琥珀葡萄球菌则属于凝固酶阴性菌。金黄色葡萄球菌和无乳链球菌作为奶牛乳房炎的主要致病菌，常常能从生产的乳汁中分离出，且金黄色葡萄球菌的阳性感染率通常会显著高于无乳链球菌，因此金黄色葡萄是导致奶牛乳房炎的首要致病菌[17-18]。然而，在本研究中仅观察到2组之间无乳链球菌存在丰度差异，而金黄色葡萄球菌的丰度不存在显著差异，由此推测该地区乳房炎流行的主要致病菌为无乳链球菌。在环境性病菌的分析中发现，2组乳汁中凝固酶阴性葡萄球菌、化脓隐秘菌、绿色气球菌丰度也出现了显著性差异，由于凝固酶阴性葡萄球菌导致的乳房炎临床症状不明显，因此该菌常常不受重视，但近年来大量研究中乳房炎奶样中分离获得的凝固酶阴性葡萄球菌甚至已经超过了金黄色葡萄球菌[19-20]，已成为许多国家奶牛乳房炎病原中的重要致病菌[21-22]。化脓隐秘杆菌是一种条件致病菌，常寄生于健康牛的乳房中，当奶牛免疫机能低下时便可导致奶牛乳房化脓性感染，继而引发乳房炎[23]。绿色气球菌作为一种机会性致病菌，常分布于奶牛场环境和皮肤中，部分菌株被证实对小鼠乳腺具有很强的致病性，说明该菌可能对奶牛乳房炎的发生具有一定的影响[24]。

综上，本试验分析发现，四川地区10家养殖场中临床型乳房炎奶牛和健康奶牛乳汁中无乳链球菌、凝固酶阴性葡萄球菌、化脓隐秘菌、绿色气球菌呈现显著差异，在奶牛乳房炎的防控中应重点关注这4种菌的感染。本研究为四川地区奶牛乳房炎的防控提供参考，为健康牛的标准制定提供微生物指标参考。