基于Illumina MiSeq的麋鹿粪样菌群多样性分析

李俊芳，白 冰，孟庆辉，程志斌，单云芳，刘 田，钟震宇

（1.北京麋鹿生态实验中心，北京，100076；2.北京市科学技术研究院资源环境研究所，北京，100089）

麋鹿（Elaphurus davidianus）是国家一级重点保护野生动物，曾在中国本土灭绝［1］。自1985 年从英国重引入以来，已建立80 余处迁地保护地，总数量近12 000 只。经种群调查发现，迁地保护麋鹿存在繁殖障碍和死亡率高等问题，个别种群发展难以维持［2-3］，消化系统疾病是导致麋鹿大量死亡的最主要原因，主要表现为由产气荚膜梭菌（Clostridium perfringens）、大肠埃希氏菌（Escherichia coli）和多杀性巴氏杆菌（Pasteurella multocida）等引起的出血性肠炎。1995—2009 年，每年都有出血性肠炎等消化系统疾病导致麋鹿死亡的报道［2］，因此，有必要对麋鹿消化道疾病进行深入研究，而了解麋鹿胃肠道菌群的结构和组成是研究的前提。

动物肠道菌群形成一个庞大而复杂的微生态平衡系统，在维持胃肠道健康方面起着重要作用［4］。近年来，对野生动物肠道微生态的研究得到了广泛关注［5］。肠道细菌大多不可培养，常规方法无法全面检测菌群组成，需用测序技术研究麋鹿消化道正常菌群分布。本研究应用Illumina MiSeq 对石首麋鹿国家级自然保护区、辽宁辽阳千山鹿场、北京麋鹿生态实验中心和河北滦河上游国家级自然保护区的32 只麋鹿粪便菌群多样性进行全面检测，旨在为今后麋鹿消化道疾病的研究提供有价值的基础数据。

1 材料与方法

1.1 样品采集

32 只麋鹿（表1）均为正常健康个体，各采集新鲜粪便20 g，快速放入自封袋内，转至液氮速冻，置-80 ℃保存。北京麋鹿生态实验中心（MLY）12 份样品，河北滦河上游国家级自然保护区（MLWC）6份样品，辽宁辽阳千山鹿场（LY）12 份样品，石首麋鹿国家级自然保护区（SS）2 份样品，其中MLY-12 为1 月龄麋鹿粪便样品，其余为成年麋鹿粪便样品。

表1 麋鹿粪便样本信息Tab.1 Information of feces samples of Père David’s deer

1.2 试验方法

1.2.1 细菌总DNA提取

每份样本取200 mg，使用粪便基因提取试剂盒（天根生化科技（北京）有限公司）提取细菌总DNA，用1%琼脂糖凝胶电泳检测DNA 的提取质量，使用NanoDrop 2000测定DNA浓度和纯度。

1.2.2 Illumina MiSeq测序流程

细菌基因组文库构建及测序均由上海美吉生物医药科技有限公司完成。具体操作步骤如下：合格的DNA 片段经电泳回收后，在Illumina MiSeq 测序平台上，使用338F（5'-ACTCCTACGGGAGGCAGCAG-3'）和806R（5'-GGACTACHVGGGTWTCTAAT-3'）对16S rRNA 基因V3—V4 可变区进行PCR 扩增（PCR仪：ABI GeneAmp®9700 型）。扩增程序：95 ℃预变性3 min，27 个循环（95 ℃变性30 s，55 ℃退火30 s，72 ℃延伸30 s），然后72 ℃稳定延伸10 min，最后4 ℃保存。将同一样本的PCR 产物混合后，使用2%琼脂糖凝胶回收PCR 产物，利用AxyPrep DNA Gel Extraction Kit（Axygen Biosciences，Union City，CA，USA）回收产物纯化，并用Quantus™ Fluorometer（Promega，USA）对回收产物检测定量。使用NEXTflexTMRapid DNA-Seq Kit（Bioo Scientific，美国）建库。构建好的文库利用Illumina 公司的MiSeq PE300 平台测序。

1.3 数据分析

通过Barcode 将数据分拣后拼接成Tag。利用QIIME软件对原始测序数据初筛，保留长度在200 bp左右的序列，将Barcode 序列和引物序列同时去除，用Mothur 软件将序列合并，去掉干扰数据，用uclust和furthest 对样品的有效序列进行聚类，以97%的一致性聚类为特征序列（OTUs），挑选每个OTU 中丰度最高的序列与Greengenes比对，再用FastTree法进行多序列比对，构建OTU表。在QIIME中，将嵌合体和Singletons 去除，计算样品序列及OTU 数，在界、门、纲、目、科、属和种水平上进行物种注释，并绘制各样品在门水平上物种组成的柱状图，运用去除嵌合体和Singletons 后的OTU 表，计算菌群多样性及丰度指数，并绘制稀释曲线，结合样本的OTU 种类及丰度，计算样品间的加权UniFrac 矩阵，利用非加权组聚类法UPGMA（unweighted pair group method with arithmetic mean）在MEGA 5.0 软件中，对样品的共享菌群进行聚类分析，用R 软件生成相应的聚类热图（heatmap）和主成分分析（PCA）图，同时计算Alpha多样性指数，包括香农-维纳指数（Shannon-Wiener index）及辛普森指数（Simpson index）。

2 结果

2.1 各样品序列及OTU

共检测麋鹿粪便样本32个，得到1 438 677条有效序列，平均每个样品为（44 959±12 348）条序列，其中MLY-4 含有效序列最多，为68 077 条，LY-2 含有效序列最少，为29 573条。序列拼接优化后，按相似度大于97%的标准对有效序列进行OTU 聚类，共获得31 459个物种分类的OTUs，平均每个样品OTU个数为（983±240）（表2）。样品质控后的序列长度为420～460 bp（图1）。

图1 麋鹿粪便样品质控后的序列长度Fig.1 Sequence length distribution of fecal samples of Père David’s deer after quality control

表2 测序序列和OTU数量统计结果Tab.2 Statistical results of sequencing and OTU quantity

2.2 样品OTU稀释曲线

从样本中随机抽取一定数量的个体，以个体数与物种数建模，比较不同样本中菌群的丰富度。图2显示，抽取样品的稀释曲线趋于平缓，表示该试验结果较好地覆盖了大部分菌群。

图2 肠道菌群丰富度稀释曲线Fig.2 Richness rarefaction curve of get microbiota

2.3 Shannon-Wiener曲线

根据样本测序量在不同测序深度构建Shannon-Wiener 指数曲线。图3 中Shannon-Wiener 曲线已经进入平台期，说明测序数据量足够大，可以反映本试验需要的绝大多数的微生物信息。

图3 肠道菌群Shannon-Wiener稀释曲线Fig.3 Shannon-Wiener index rarefaction curve of get microbiota

2.4 核心菌群结构及组成分析

用QIIME 软件对样品不同分类水平进行物种注释，对菌群结构进行详细的分类整理，共获得12 个门，74个属。

在门水平上，物种注释共获得11 个门和1 个未分类的门，检测到的99.7%微生物群落为细菌，未标记的微生物仅占0.3%。优势菌群（相对丰度大于1%）从高到低分别为厚壁菌门（Firmicutes）、拟杆菌门（Bacteroidetes）、变形菌门（Proteobacteria）、螺旋体门（Spirochaetes）、放线菌门（Actinobacteria）、软壁菌门（Tenericutes）、疣微菌门（Verrucomicrobia）、梭杆菌门（Fusobacteria）和其他（others），其中拟杆菌门和厚壁菌门为优势菌群。图4 显示，不同样品的菌群门类基本相同，但是组成比例不同，同一菌门的相对丰度差异较大。

图4 麋鹿肠道菌群在门水平上的物种相对丰度Fig.4 Species relative abundance of gut microbiota in Père David’s deer at the phylum level

在属水平上，共获得74 个属，其中27 个（36.5%）尚未被注释出属名。优势菌属（相对丰度大于1%）一共有15 种，占89%，从高到低分别为瘤胃球菌属（Ruminococcus）、不动杆菌属（Acinetobacter）、拟杆菌属（Bacteroides）、普雷沃氏菌属（Prevotella）、鞘氨醇杆菌属（Sphingobacterium）、副拟杆菌属（Parabacteroides）、假单胞菌属（Pseudomonas）；而双歧杆菌属（Bifidobacterium）和乳杆菌属（Lactobacillus）检测的相对丰度低，在MLY-12 粪便样品中却是优势菌属。

2.5 样品间共享菌群分析

运用Heatmap 将样本间丰度相似性进行聚类，得到共享菌群分析图（图5），其中MLY-12样本的细菌群落组成与其他样本明显不同，LY-12 的菌群组成与其他样本的差异性也较大。

图5 样品间共享菌群分析Fig.5 Analysis of shared microbiota between samples

2.6 主成分分析

将不同样本OTU（97%相似性）主成分组成的差异反映在二维坐标图上（图6）。所有样品被分为3部分，其中SS 的2 个样品单独聚在一起，LY 的大部分样品聚在一起，而MLY、MLWC 和部分LY 系列样品聚在一起。

图6 麋鹿粪便样品微生物主成分分析Fig.6 Principal component analysis（PCA）in the fecal sample of Père David’s deer

2.7 菌群多样性

样品的多样性结果见表3。群体内香农-维纳指数为2.42～5.83，其中MLY-12 最低，LY-3 最高。辛普森指数为0.007 3～0.148 6，最高为MLY-12（0.148 6），最低为MLY-4（0.007 3）。

表3 样品间的Alpha多样性Tab.3 Alpha diversity index of between samples

3 讨论

本研究利用Illumina MiSeq测序检测了4个麋鹿栖息地的32 只麋鹿粪样微生物菌群多样性，获得了较全面的菌群信息。Shannon-Wiener 曲线已经进入平台期，说明测序数据量足够大，可覆盖32 只麋鹿的肠道微生物，能真实全面地反映正常肠道菌群组成。

OTU 聚类后共检测出11个门和1个未分类的门类，麋鹿肠道微生物主要来自拟杆菌门和厚壁菌门，这与前人对瘤胃微生物的研究结果［6-10］一致，Sundset 等［9-10］对挪威地区驯鹿（Rangifer tarandus）瘤胃微生物区系的研究表明，驯鹿瘤胃内细菌主要由拟杆菌门和厚壁菌门细菌组成。麋鹿和其他有蹄类动物粪便中的优势菌群为厚壁菌门和拟杆菌门，可能是因为它们以粗饲料为食［11］，胃肠道中的拟杆菌和厚壁菌在粗纤维的消化过程中发挥着重要作用［12］。在属水平上还有大量尚未命名的菌属，这与目前麋鹿微生物的研究水平有限有关，随着测序水平的不断进步，未来将会检测认知更多的菌属。

样品OTUs 的丰度显示MLY-12 的菌群丰度最低，样品间共享菌群的分析也显示MLY-12 与其他样品的共享菌群最少，同时其香农-维纳指数也最低。该样品采集于麋鹿中心1 只人工喂养的1 月龄麋鹿，研究结果表明1 月龄麋鹿与成年麋鹿的肠道菌群存在很大的差异，这个结果也有待进一步分析。

主成分分析中，所有粪便样品被分为3 部分，其中，石首麋鹿国家级自然保护区采集的粪便样品、大部分辽宁辽阳千山鹿场采集的粪便样品分别聚在一起，而北京麋鹿生态实验中心、河北滦河上游国家级自然保护区和部分辽宁辽阳千山鹿场粪便样品聚在一起。这可能与迁地保护时间长短有关（表1），迁地后，随着采食条件的不同，迁地的麋鹿肠道菌群有了进一步的分化。北京麋鹿生态实验中心自1985 年成立以来，致力于扩大麋鹿繁育种群及建立迁地种群，迄今麋鹿中心建立的麋鹿迁地保护场所占我国总麋鹿保护场所的71.74%，成为我国麋鹿种群恢复和迁地保护的主要技术支撑单位，发挥着举足轻重的作用。本研究采集了石首麋鹿国家级自然保护区、辽宁辽阳千山鹿场、河北滦河上游国家级自然保护区的麋鹿粪样，聚类分析结果表明，随着迁地时间缩短，聚类结果越集中。河北滦河上游国家级自然保护区于2016年重引入麋鹿10只，麋鹿群体形成时间较短，肠道菌群改变尚不明显，因此，河北滦河上游国家级自然保护区与麋鹿中心麋鹿粪便菌群聚类在一起，而湖北石首麋鹿种群形成时间最长，肠道菌群分化最明显，与其他样品的分离距离最远，区分也最明显。

4 小结

本研究应用Illumina MiSeq 对4 个麋鹿栖息地32 只麋鹿粪便样本的菌群多样性进行分析，共检测出1 438 677 条有效序列，平均每个样品有（44 959±12 348）条有效序列；在97%相似度条件下共获得31 459 个物种分类的OTUs，平均每个样品为（983±240）个OTUs；平均读长为426 bp。

在门的分级水平上，粪便细菌主要来自11 个门，其中厚壁菌门和拟杆菌门为优势菌门。不同迁地间麋鹿粪样细菌多样性存在差异。聚类结果可将粪便样品聚为3 个群落，随麋鹿迁出时间的延长，麋鹿肠道菌群聚类分析的远近越明显，这可能与迁出地的麋鹿取食改变有关。