基于宏基因组测序技术对绵羊呼吸道微生物群落的研究

2022-11-17刘海霞韩大勇朱爱文周明夏王步忠耿昕之

中国饲料 2022年18期

刘海霞，韩大勇，闫伟，朱爱文，周明夏，王步忠，吴昊，耿昕之

（1.江苏农牧科技职业学院，江苏泰州 225300； 2.江苏西来原生态农业有限公司，江苏泰州 225328； 3.江苏省阜宁县农业农村局，江苏盐城 224400）

宏基因组学是以高通量测序技术为基础的组学研究方法，宏基因组测序技术可以全方位测定分析样品中所含微生物的种群结构、进化关系、功能活性和多样性等。郭玉梅等（2022）利用宏基因组测序技术研究发热病人咽部微生物菌群，鉴定出4079个物种，归属于97个门，951个属，发现试验组较对照组存在消化链球菌、拟杆菌、支原体等物种差异，β-内酰胺类酶类基因为耐药基因。宏基因组技术可以通过挖掘微生物基因组信息，分析其功能及代谢通路，进而挖掘菌群的种群结构、功能特性及种群间的相互关系（彭玺等，2022；孙波等，2021；许悦等，2020；马海霞等，2015；贺纪正等，2008）。本研究采用宏基因组测序技术，对绵羊呼吸道样本进行微生物菌群多样性的测定分析，指导防治绵羊呼吸道感染性疾病，为绵羊健康养殖提供数据支撑。

1 材料与方法

1.1 材料

1.1.1 样品采集采集8只试验用绵羊鼻涕样本（健康羊、呼吸道症状病羊各4只）。鼻拭子轻轻蘸取鼻腔中间段涕样，采集完成后，将鼻拭子插入预置2 mL PBS溶液的进口冻存管中，液氮速冻，-20℃保存。

1.1.2 主要仪器、试剂 Illumina测序仪（6000，NovaSeq），高效拭子基因组DNA提取试剂盒（WH0215，离心柱型，百奥莱博），TIANSeq DNA Library Prep Kit文库构建试剂盒（NG103，天根生物），TruSeq PE Cluster Kit v3-cBot-HS（PE-401- 3001，Illumina）。

1.2 方法

1.2.1 样品DNA提取采用基因组DNA提取试剂盒提取样本中的基因组DNA，然后分析提取的基因组DNA的纯度和完整性（琼脂糖凝胶电泳），再DNA浓度和质量进行检测（Nanodrop 2000）。

1.2.2 文库构建与检测用超声波破碎仪将检测合格的DNA样品随机打断，片段长度约为350 bp，修复补平打断的DNA片段末端，在3’端添加碱基“A”，将片段末端转换为黏性末端，利用碱基互补在黏性末端两侧添加包含Index序列的DNA接头。磁珠筛选收集一定长度范围内的目的片段，PCR扩增，在目的片段末端添加Index，构建完成测序文库。使用Aglient 2100生物芯片分析系统检测文库片段大小，并使用Step One PlusTM Real-Time PCR System，以P5、P7接头作为引物进行QPCR定量。

1.2.3 上机测序将检测合格的不同测序文库按照目标下机数据量的需求和有效浓度进行pooling，在Illumina测序平台上运行双端测序程序，得到150 bp的双端测序reads。

1.2.4 数据分析将测序得到的原始数据进行质控后得到有效数据。使用CD-HIT、MetaGeneMark、Bowtie2等软件进行所有样品的reads的Metagenome组装、物种注释、功能数据库注释、抗性基因注释等。DIAMOND软件将Unigenes与从NCBI NR数据库中抽提出的微生物序列进行比对，genes与各功能数据库进行比对进行物种注释，DFAMOND软件将Unigenes与各功能数据库进行比对，进行功能注释，CARD数据库提供的Resistance GeneIdentifier（RGI）软件将Unigenes与抗性基因数据库（CARD）进行比对，完成抗性基因注释。

2 结果与分析

2.1 测序结果概述本研究测序在Illumina NovaSeq测序平台完成，共获得625347.73 Mbp的原始数据，经过质控得到536874.34 Mbp的有效数据，混合组装单样品、多样品，共得到989926098 bp的Scaftigs。使用MetaGeneMark软件进行基因预测，共得到561403个开放阅读框（ORFs），去冗余后共获得418423个ORFs，完整基因的个数为96776，所占比例为20.13%。非冗余基因集注释到属的比例为72.21%，非冗余基因集注释到门的比例83.36%。使用DIAMOND软件对非冗余基因集进行KEGG、eggNOG、CAZy、CARD等功能数据库进行注释，有572318个ORFs比对到KEGG数据库，560255个ORFs比对到eggNOG数据库，10112个ORFs比对到CAZy数据库。通过CARD注释，比对到CARD数据库有299个基因。

2.2 物种注释通过物种丰度分析，试验组共得到103个门，159个纲，413个目，840个科，1789个属，4790个种；对照组共得到99个门，159个纲，409个目，836个科，1780个属，4783个种。试验组与对照组比较发现，变形菌纲丰度对照组显著高于试验组；乳杆菌目丰度对照组显著高于试验组，而杆菌目丰度试验组显著高于对照组；葡萄球菌科丰度试验组显著高于对照组，链球菌科丰度对照组显著高于试验组；乳球菌属丰度试验组显著高于对照组，葡萄球菌属丰度对照组显著高于试验组；金黄色葡萄球菌丰度试验组显著高于对照组，格氏乳球菌丰度对照组显著高于试验组。综合分析，试验组致病菌丰度相对较高，对照组益生菌丰度相对较高。

2.3 功能注释

2.3.1 样品KEGG注释结果本研究发现，样品中636981个不同基因被注释，从属于KEGG数据库6大类46个KEGG通路中，分别是新陈代谢（187129，29.37%）、人类疾病（173291，27.20%）、生命系统（169859，26.66%）、细胞进程（48558，7.62%）、环境信息处理（29505，4.63%）、遗传信息处理（28639，4.52%）。富集差异表达基因最多的6个KEGG通路，分别是信号传导（52127）、氨基酸代谢（42336）、免疫系统（27816）、转运与分解代谢（20982）、传染性疾病（20029）、感知系统（17920）。环境信息处理类通路中的信号传导通路、新陈代谢类通路中氨基酸代谢通路所含基因数最多。（如图1）

图1 基因KEGG代谢通路统计

2.3.2 样品eggNOG数据库注释结果 eggNOG数据库注释发现，本试验相对含量在前5位的功能大类是细胞内转运、分泌和囊泡转运，翻译后修饰、蛋白质翻转，复制、重组与修复，信号转导机制，转录，基因数分别为50128、24833、12680、12635、11206个。eggNOG功能分类统计结果（如图2）

图2 基因eggNOG功能分类

2.3.3 样品CAZy数据库注释结果经过与CAZy数据库比对，本试验样品中共鉴定出6905个碳水化合物活性酶，包含3238个糖基转移酶，2217个糖苷水解酶类，数量最多，占据样品检测数据总数的77%。1002个碳水化合物结合模块，258个碳水化合物酯酶，112个辅助氧化还原酶，78个多糖裂合酶最少，占据样品检测数据总数的0.63%。

2.3.4 样品抗性基因注释结果利用CARD数据库进行抗生素抗性基因的预测，绘制ARO分布热图，可反映各抗性基因类型ARO在各样品中的分布情况。如图3可知，葡萄球菌和埃希氏杆菌属耐药基因分布最广，存在于测定的8个样品中，耐药基因APH（3）-Ib只存在于样品D2中，耐药基因novA、antibiotic、msbA、Mycobacterium、macB广泛存在于除样品D1外的其他样品中。耐药基因smeR、aadA只存于B2和D1样品中。

图3 ARO分布热图

3 讨论

3.1 试验样品的物种注释随着生物信息学和现代分子生物学的发展，宏基因组测序技术的应用领域更广。为探究试验样本微生物组成特征，使用Diamond软件将unigene序列与NR数据库比对进行物种注释，可获得样本在分类水平各物种组成及丰度信息，通过物种丰度分析反映试验样本微生物结构组成多样性变化情况。李娟（2020）等预测到麦洼牦牛瘤胃微生物基因2065638个开放阅读框（ORFs）完整基因数为333371个，舍饲组厚壁菌门、拟杆菌门、软壁菌门、浮霉菌门相对放牧组丰度均增加，纤维杆菌门相对丰度舍饲组显著低于放牧组。本研究发现，绵羊呼吸道微生物菌群多样性特征明显，在门水平上，试验组优势菌群主要为厚壁菌门、拟杆菌门、变形菌门、放线菌门、螺旋体门，对照组优势菌群主要为厚壁菌门、拟杆菌门、放线菌门、梭杆菌门、变形菌门。

3.2 试验样品的功能注释通过对基因在KEGG生物通路数据库上进行注释分析，可以查看基因参与的生物通路、体现的生物功能。对比分析eggNOG数据库与基因翻译获得的蛋白质序列，获得相应的蛋白质功能注释信息。将基因序列与CAZy数据库进行比较，得到碳水化合物活性酶类的物种来源、酶功能EC分类、基因、蛋白序列及蛋白结构信息。Pitta等（2016）利用宏基因测序研究泌乳期奶牛瘤胃微生物群，发现基因序列对碳水化合物蛋白质代谢等近50%的代谢有重要作用。本研究从KEGG数据库分别注释到636981个基因，环境信息处理类通路中的信号传导通路所含基因数最多，其次是新陈代谢类通路中氨基酸代谢通路所含基因量。CAZy数据库比对发现，样品中共鉴定出碳水化合物活性酶6905个，其中，糖基转移酶3238个和糖苷水解酶类2217个，数量最多，占据样品检测数据总数的77%。

3.3 试验样品的抗性基因注释滥用抗生素会导致动物体内产生抗生素耐药性细菌，Ma等（2015）在猪、鸡和人类粪便中检测出了含量较高的四环素和红霉素等抗生素。Hitch等（2018）在绵羊瘤胃中检测到达托霉素和粘菌素。利用CARD数据库进行抗生素抗性基因预测，能反映各抗性基因类型ARO物种中的分布情况。本研究发现，葡萄球菌和埃希氏杆菌属耐药基因分布最广，存在于检测的8个样品中，耐药基因APH3-1b只存在于样品D2中，耐药基因novA、antibiotic、msbA、Mycobacterium、macB广泛存在于除了样品D1以外的其他样品中，耐药基因smeR、aadA只存于B2和D1样品中。