用18S rRNA 高通量测序分析比较万峰湖丰水期和枯水期生物多样性

秦燕，范波，陈应美，卢彪

（兴义民族师范学院生物与化学学院，贵州 兴义 562400）

万峰湖位于黔、滇、桂三省交界处，因建设天生桥电站拦截南盘江而成的淡水湖泊，积水后水位深度约170 m，最大库容可达108 亿m3[1]。万峰湖是珠三角地区重要的水源供给地，因而水质的好坏直接影响周围地区的用水安全。目前有关万峰湖的研究大多集中在水质污染方面，特别是工业废水排放导致的汞污染。对于万峰湖浮游生物资源调查的研究较少，并且研究手段也仅限于传统镜检的分析方式[1-3]。目前国内对于万峰湖所有真核生物资源的研究几乎处于空白阶段。本研究对万峰湖真核生物多样性进行分析研究，了解其资源分布，将对今后万峰湖真核生物多样性研究提供参考。

本研究于2020 年10 月（丰水期）和2021 年5月（枯水期）分别对万峰湖上游、中游和下游湖中进行采样调查，利用高通量测序技术，分析并比较万峰湖真核生物的群落结构组成，通过多样性指数和均匀度分析、PCA 和物种丰度热图分析万峰湖真核生物群落结构多样性，为万峰湖水质的评价和后续富营养化的防治提供数据支持和理论依据。

1 材料与方法

1.1 采样点的选择

万峰湖是河流型湖泊，充分考虑到水流的方向和水位的特点，于2020 年10 月（丰水期）和2021年5 月（枯水期）分别对万峰湖上游、中游和下游湖中的A-L 共12 个采样点浅水区采样进行镜检鉴定，而后对DB（大坝）、MJW（梅家湾）、PM（坡母）、SL（石林）、ZYF（榨油坊）、BJ（巴结）、YYT（野鸭滩）、GB（革布）分别标记为1～8 的八个位置不同深度的水体进行采样，HC 的取样位置在SL 附近，采样位置分布见图1。表层水体的位置为0～2 m 水深处，中层水体的位置为5～10 m 水深处，底层水体的位置为15 m 水深左右，同一位置不同水深的水采样后混合。

图1 万峰湖采样点分布图Fig.1 The map of sampling sites distribution in Wanfeng Lake

1.2 样本采集与预处理

1.2.1 样本采集

用5 L 带盖的不锈钢采样器采集不同水深位置的水样，分别装于5 L 的干净矿泉水桶中。溶解氧（DO）和水温（WT）用溶氧仪（HQ30 型，美国哈希）现场检测；pH 用pH 试纸检测。

1.2.2 样本预处理

从万峰湖8 个采样点，分别采集水样2 000 mL。水样用真空泵抽滤，滤膜直径为0.45 μm，将载有样本的滤膜立即放在-80℃冰箱保存。

1.3 高通量测序

取适量滤膜上的样品，液氮冷冻后，组织破碎仪中进行破碎。取适量洗过之后的液体12 000 r/min离心2 min，留取沉淀进行提取。以OMEGA 试剂盒E.Z.N.ATMMag-Bind Soil DNA Kit 提取试剂盒（M5635-02）提取DNA 后质检。利用Qubit3.0 DNA检测试剂盒对基因组DNA 精确定量,以确定PCR反应DNA 模板的量。PCR 所用的引物如下：扩增区域为18SV4；测序平台Miseq 2x300 bp，18SV4 引物序 列 18SV4F：GGCAAGTCTGGTGCCAG；18SV4R：ACGGTATCTRATCRTCTTCG。

巢式PCR 进行两轮PCR 扩增，反应体系30 μL：2×HieffRobust PCR Master Mix 15 μL，引物F 和引物R 各1 μL，模板20～30 ng，加水9～12 μL 至总体积30 μL。PCR 程序如下，第一轮PCR 反应条件为：94℃，3 min；然后94℃30 sec、45℃20 sec、65℃30 sec 5 个循环，然后94℃20 sec、55℃20 sec、72℃30 sec 20 个循环，72℃延伸5 min，10℃保存。第二轮PCR 反应条件：95℃3 min，然后94℃20 sec、55℃20 sec、72℃30 sec 5 个循环，72℃延伸5 min，10℃保存。

Illumina MiseqTM测序：PCR 产物通过2%琼脂糖凝胶电泳检测并回收，用Qubit3.0 荧光定量仪进行文库浓度测定，而后送往上海生工进行高通量测序，测序平台为Illumina MiseqTM。

1.4 高通量测序数据分析和处理

1.4.1 数据预处理

下机测序得到的是双端序列数据含有barcode序列以及测序时加入的引物和接头序列，首先需要使用Cutadapt（v1.18）去除引物接头序列[4]，再根据PE reads 之间的overlap 关系，使用PEAR（v 0.9.8）将成对的reads 拼接（merge）成一条序列[5]，然后按照barcode 标签序列识别并区分样品得到各样本数据，最后使用PRINSEQ（v 0.20.4）对各样本数据的质量进行质控过滤[6]，得到各样本有效数据。

1.4.2 OTU 聚类

OTU（Operational Taxonomic Units）是在系统发生学或群体遗传学研究中，人为给某一个分类单元（品系、属、种、分组等）设置的统一标志。通过Usearch（v 11.0.667）软件聚类分析物种的种类和数目。将序列按照彼此的相似性分归为许多小组，一个小组就是一个OTU。对所有序列进行OTU 划分，通常对97%相似水平的OTU 进行生物信息统计分析[7,8]。

1.4.3 OTU 物种注释和统计

为了得到每个OTU 对应的物种分类信息，利用blastn 将序列与对应数据库进行比对，筛选出序列的最佳比对结果，并对比对结果进行过滤，默认满足相似度＞90%且覆盖率＞90%的序列被用来后续分类。18S 使用Blast 比对NCBI 18S 数据库（http://ncbi.nlm.nih.gov/）。最终分别在各个分类水平：domain（域）、phylum（门）、class（纲）、order（目）、family（科）、genus（属）、species（种）上统计各样本的群落组成。

1.4.4 OTU 的分析

基于OTU 丰度，每个取样点的OTU 数目的Venn图由R 的VennDiagram 软件R（V1.6.20）绘制[9]；为了展示不同取样点OTU 组成的不同，进行了主成分的分析（principal component analysis，PCA）。此步骤运用的软件是R 中的“ade4”包（v 1.7-13）[10]；OTU等级丰度曲线能可视化的得到物种的丰富程度和均匀度。物种的丰富度可以从X 轴上不同物种的数目得到。图中X 轴表示不同OTU 的等级，Y 轴表示OTU 的相对丰度，等级曲线绘制来自软件R 中的gplots 包（v3.0.1.1）；物种的注释中每个分类阶元（域、门、纲、目、科、属、种）上的OTU 的数目或者不同样本中OTU 的数目都是基于NCBI 数据库，RDP分类2.2[11]，贝叶斯算法得出的。通过比较每个OTU中的代表序列和数据库中的序列，进行相似性比对和物种注释，从而在不同的分类阶元的水平上计算每个样本中的群落组成；多样性指数分析可以估计环境群落的物种丰度和多样性[12]。本研究中计算群落分布丰度（Community richness）的指数有：Sobs—观察到的丰度（http://www.mothur.org/wiki/Sobs），Chao1估算（http://www.mothur.org/wiki/Chao），ACE估算（http://www.mothur.org/wiki/Ace）。计算群落分布多样性（Community diversity）的指数有：Shannon 指数（http://www.mothur.org/wiki/Shannon），Simpson 指数（http://www.mothur.org/wiki/Simpson）和覆盖率（http://www.mothur.org/wiki/Coverage）。计算群落分布均匀度（Community evenness）的指数有：Shannoneven—基于Shannon 指数的均匀度测量（Pielou 均匀度指数J）；稀释性曲线中数据计算使用软件Mothur（v1.43.0）完成[12]，稀释性曲线的绘制通过软件R 中的vegan（v 2.5-6）包完成；β 多样性分析主要用来评估样本间物种多样性的不同。β 多样性分析是通过软件phyloseq（v.1.30.0）完成的[13]，使用R 的gplots package 进行距离热图绘制；主坐标分析中所用软件为phyloseq；共线性关系图由软件R 的circlize 包绘制[14]。

1.5 传统镜检验证

选取万峰湖上游、中游和下游A-L 共12 个采样点浅水区采样进显微镜镜检，对各采样点的藻类组成和数量进行统计。

2 结果与分析

2.1 高通量测序分析结果

2.1.1 OTU 的分析

OTU 的聚类分析主要是选出与代表序列相似性在97%以上的序列。在丰水期的样本中共得到363 个OTUs，OTU 的韦恩图显示七个取样点的样本共同有的OTUs 有238 个，另外YYT 还有2 个特有的OTUs，PM 和SL 各有1 个OTU（图2-A1）。在枯水期的样本中252 个OTUs，韦恩图中显示三个样本中共有的有131 个OTUs，另外BJ 和HC 共有的有36 个，BJ 和ZYF 共有的有28 个，ZYF 和HC 共有的有15 个，而BJ、HC 和ZYF 特有的分别有9个、6 个和27 个（图2-A2）。

稀释性曲线主要用来评估测序产生的数据量是否能足以覆盖种群中所有物种。当曲线趋于平缓时，说明测序的数据量比较合理，Alpha 多样性指数的稀释曲线（图2-B1、B2）均表明稀释曲线渐渐趋于平缓，说明样本的测序数据量足够充分，测序深度基本全覆盖，也可以说明在我们测序的样本中基本没有物种检测不到，这保证了我们研究数据的可靠性。

基于物种相对丰度的主成分分析表明，在我们取样的不同样本中种群的组成有一定的差别，如丰水期中主成分PC1 和PC2 的值为50.31%和22.66%，而在枯水期中主成分PC1 和PC2 的值分别为66.29%和33.71%（图2-C1、C2）。

图2 丰水期（A1-D1）和枯水期（A2-D2）样本中OTU 分析Fig.2 OTU analysis of samples in wet season（A1-D1）and dry season（A2-D2）

OTU 等级图表主要从两个方面展示物种的多样性：物种的丰度和物种的均匀度。在水平方向上，曲线的宽度反映物种的丰度，曲线跨越的水平轴越宽，物种的丰度越高。曲线的走向越平缓，说明物种的分布差别越小，曲线倾斜度越大，物种分布的多样性越大。从丰水期和枯水期的曲线比较，可以得出：丰水期曲线的跨度宽，则其物种的丰富度高，枯水期曲线的倾斜度大，则其物种的多样性大，均匀度低。

2.1.2 Alpha 多样性分析

Alpha 多样性主要是对样本中物种多样性进行分析，利用OTU 物种和丰度计算出Ace、Chao、Shannon、Simpson 指数和Shannoneven。Chao 和Ace是丰度指数，也即OTU 的数目；Shannon、Simpson 和Coverage 主要用来描述群落分布多样性特征，因此Shannon 指数越大，Simpson 指数越小，则样本中物种多样性越高；、Shannoneven 是群落分布均匀度指数。

比较丰水期和枯水期的数据，丰水期的Chao，Ace 指数明显高于枯水期，说明丰水期物种的丰度要比枯水期丰富的多。丰水期的Shannon 指数低于枯水期的Shannon 指数，而丰水期的Simpson 指数高于丰水期的Simpson 指数，且差异显著，说明丰水期群落分布多样性低于枯水期（见表1）。而丰水期物种的均匀度指数低于枯水期，则群体的多样性低于枯水期。

表1 丰水期和枯水期样品Alpha 多样性指数比较表Tab.1 Comparison of alpha diversity index of samples in wet season and dry season

2.1.3 Beta 多样性分析

为了进一步展示样本中物种多样性的不同，用主坐标分析来展示不同样本见的不同（图3）。如果两个样本距离很近，则说明两个样本的物种组成很相似。图3 中不同形状的点代表不同的样本，PCoA1指引起的样本的最大的不同的主要的坐标成分，解释度值分别为丰水期40.83%和枯水期75%。PCoA2是其次，解释度值分别为24.17%和25%。由图3 可知，丰水期样本中，样本DB、MJW 与BJ 距离较近，SL 和PM 距离较近，而样本GB 与YYT 距离较近，说明DB、MJW 与BJ 样本中的物种组成相似，SL 和PM 样本、样本GB 和YYT 中的物种组成相似，枯水期的三个样本之间距离均较远，说明它们之间的群落组成均差别较大。

图3 丰水期和枯水期主坐标分析比较Fig.3 Comparison of principal co-ordinates analysis in wet season and dry season

2.1.4 在门水平上比较丰水期和枯水期真核浮游生物的组成

在门的水平上比较，丰水期样本中含有22 种真核浮游生物，前五种分别为节肢动物门，绿藻门，壶菌门，硅藻门和子囊菌门。枯水期样本中含有15种真核浮游生物，前五种分别为绿藻门，节肢动物门，轮虫动物门，壶菌门，硅藻门和纤毛亚门。两组样本中真核浮游生物的种类和相对丰度见图4，优势物种丰度热图见图5。

图4 丰水期（A）和枯水期（B）真核浮游生物的种类和相对丰度比较（在门的水平上）Fig.4 Comparison of species and relative abundance of eukaryotic plankton in wet season（A）and dry season (B）(at the phylum level）

图5 丰水期（A）和枯水期（B）真核浮游生物优势物种丰度热图（在门的水平上）Fig.5 The heatmap of dominant species abundance of eukaryotic plankton in wet season (A）and dry season (B）(at the phylum level）

在门的水平上，比较同一取样点样本在枯水期和丰水期的物种种类和比例，显示出明显的差异。以BJ 为例，节肢动物门在丰水期的比例为45.46%而枯水期的比例为39.54%；绿藻门从丰水期的2.05%上升到枯水期的11.12%；硅藻门在丰水期没有检测到，而在枯水期比例为2.44%（图6）。

图6 以BJ 为例，比较丰水期（A）和枯水期（B）真核浮游生物的种类及百分比（在门的水平上）Fig.6 Comparison of the species and percentage of eukaryotic plankton in wet season (A）and dry season(B）(at the phylum level）,taking BJ as example

2.1.5 在种的水平上比较丰水期和枯水期真核浮游生物的组成

在种的水平上比较，丰水期的七个样本中共得到147 种真核浮游生物，其中145 种是所有样本共有的核心物种，PM 和SL 各有1 个样本特有的物种；在枯水期的三个样本中，共得到187 种真核浮游生物，其中108 种是三个样本共有的核心物种，有21 种是BJ 和HC 所共有，19 种是BJ 和ZYF 共有，10 种是HC 和ZYF 共有，另外BJ、HC 和ZYF 特有的物种分别有6 种、5 种和17 种（图7）。

图7 丰水期（A）和枯水期（B）真核浮游生物物种分布Venn图（在种的水平上）Fig.7 Venn diagram of eukaryotic plankton species distribution in wet season (A）and dry season (B）(at the species level）

比较丰水期和枯水期物种，两组数据都含有的物种有马索隐藻、拟蛋白核隐藻和光滑鱼鳞藻三种。丰水期的优势物种为右突新镖水蚤、甲藻、马索隐藻、海洋尾蚴和中剑水蚤，所占的比例分别为41.8%、10.3%、8.3%、6.92%和3.28%。枯水期的优势物种为真镖水蚤、马索隐藻、拟二叉角甲藻、辐球藻和蛋白核隐藻，所占的比例分别为35.67%、23.68%、4.97%、4.47%和2.83%（图8）。

图8 以BJ 为例，丰水期（A）和枯水期（B）真核浮游生物物种的种类和百分比（在种的水平上）Fig.8 Species and percentage of eukaryotic plankton in wet season (A）and dry season (B）(at the species level）,taking BJ as example

2.1.6 共线性关系

在门的水平上分析比较丰水期和枯水期，不同样本中优势物种的共线性关系。图9 仅显示丰度最高的前10 个物种分类，共线性分析的结果表明，丰水期7 个取样点已知的真核浮游生物中，节肢动物门的丰度最高，7 个取样点的比例相当，其次是绿藻门和壶菌门；而在枯水期的样本中，绿藻门的丰度最高，取样点HC 中的数量最多，其次是节肢动物门，取样点BJ 中的数量最多，再次是轮虫动物门，取样点ZZF 中的数量最多。

图9 丰水期（A）和枯水期（B）真核浮游生物的共线性关系比较图Fig.9 Comparison of collinearity between eukaryotic plankton in wet season（A）and dry season（B）

2.2 镜检验证丰水期和枯水期绿藻的种类和丰度

比较枯水期和丰水期优势物种，发现绿藻门是共有优势物种，并且真核浮游生物中绿藻指数还是检测水质营养型的重要指标[2]。通过显微镜镜检，从图1 中A-L 共12 个采样点浅水区采样进行定性定量分析，共检测到9 种绿藻，并对其丰度进行统计，网状水网藻、水绵属、集团刚毛藻、衣藻为万峰湖优势种（表2）。而后分别在丰水期和枯水期取样比较其分布情况，网状水网藻、水绵和衣藻在较多取样点有分布。拟新月藻原变种和独球藻在枯水期分布多，丰水期只有一个取样点检测到；而隆起葡萄藻原变种只在枯水期的4 个位点有分布，丰水期未检测到（表3）。

表2 万峰湖12 个取样点绿藻的组成Table 2 The species of green algae in 12 sampling sites of Wanfeng Lake

表3 万峰湖枯水期和丰水期9 种绿藻分布比较Fig.3 The distribution of nine species of green algae in Wanfeng Lake in dry season and wet season

3 讨论

本研究的创新点在于首次利用高通量测序技术对万峰湖丰水期和枯水期真核浮游生物多样性调查，对OTUs 的分析如Alpha 多样性指数稀释曲线充分说明了测序深度足够，数据合理。通过比较万峰湖枯水期和丰水期群落分布多样性和群落分布均匀度发现，不同时间的样本存在显著的差异，说明湖体季节的变化会影响浮游生物的数量和种类。沈韫芬[15]在《微型生物检测新技术》书中曾指出生物多样性指数越大，水质越好。Shannon 指数2.0～3.0 代表水质轻度污染，1.0～2.0 代表水质中度污染，我们的数据显示万峰湖枯水期和丰水期的Shannon 值分别为3.03 和2.62，属于轻度污染或无污染状态。

镖水蚤隶属于节肢动物门桡足类，在万峰湖丰水期和枯水期均为优势物种，是鱼类的饵料，多分布于营养型水体中。甲藻和隐藻是鞭毛藻类，含量也比较丰富，大多分布在富有有机质的水体中，枯水期取样时间在5 月，随着光照增多，湖水中浮游植物（如隐藻）快速大量繁殖，也推动了浮游动物如镖水蚤的增长。

万峰湖的藻类主要是绿藻，绿藻门从丰水期的2.05%上升到枯水期的11.12%，说明枯水期的湖水条件更利于藻类的繁殖。绿藻指数还是检测水质营养型的重要指标，我们还通过镜检比较丰水期和枯水期不同地点的绿藻组成，结果发现枯水期的绿藻种类和数量均明显高于丰水期。传统镜检分析时，为了保证样本的全面性，在不同季节的多个取样点取样，然而镜检仅鉴定出9 种绿藻，显著低于高通量测序鉴定的数量，这说明通过镜检的形态鉴定的方法鉴定出的生物的种类较少[1-3]，随机性大，存在人为误差，容易漏检一些数量少的物种，不能宏观和全面的鉴定生物的种类和数量。而高通量测序检测生物多样性的研究较为广泛的利用[16-20]，该技术灵敏度高，能检测出丰度低于万分之一的痕量物种。测序技术已经被比喻为“新时代显微镜”，通过测序和序列比对，可快速准确处理海量数据，减少人为误差。

很多生物多样性的研究都会结合多种影响因子对群落多样性进行研究。如陈美军[19]对太湖中处于不同营养水平湖区的真核浮游生物多样性进行研究；张莉[20]用ITS 高通量测序技术研究黄海中微型真核浮游生物多样性，发现温度和经纬度都会真核浮游生物多样性；李磊[21]研究了环境因子对万峰湖浮游生物时空分布的影响。本研究未能结合这些影响因子进行分析，后续的工作将会结合这些因子进行深入分析。

万峰湖的浮游生物的调查已有很多研究报道，之前的研究大多以浮游植物为指标评价万峰湖的水质和污染情况。本研究利用高通量测序技术，用18S rRNA 第一次全面调查了万峰湖真核浮游生物的组成，并比较了丰水期和枯水期真核浮游生物的组成，并对水质有指示作用的绿藻镜检鉴定，全面的调查了万峰湖真核浮游生物的组成，虽然未能结合一些影响因子分析，但为保护万峰湖提供科学依据，并对近几年对万峰湖的治污效果做初步的评价。

结论

（1）通过两次调查分别得到丰水期真核浮游生物16 个门，147 种，枯水期23 门，187 种。丰水期的优势物种为浮游节肢动物，主要包括右突新镖水蚤、甲藻和马索隐藻；枯水期优势物种为浮游节肢动物和绿藻，主要包括真镖水蚤、马索隐藻和拟二叉角甲藻。而镖水蚤在丰水期和枯水期均属于优势物种。

（2）比较丰水期和枯水期的Ace、Chao、Shannon、Simpson 指数和Shannoneven 数据，发现丰水期物种的丰度要比枯水期丰富，而丰水期群落分布多样性低于枯水期，且丰水期物种的均匀度指数低于枯水期。

（3）万峰湖枯水期和丰水期的Shannon 值分别为3.03 和2.62，属于无污染和轻度污染状态，说明近几年万峰湖的治污效果明显。