一种基于环境dna宏条形码技术的底栖动物完整性评价方法
技术领域
本发明属于生态环境评价与治理领域，更具体地说，涉及一种基于环境DNA宏条形码技术的底栖动物完整性评价方法。
背景技术
环境DNA(eDNA)是指在环境样品中所有被发现的不同生物的基因组DNA的混合，包括水样DNA、沉积物DNA、生物源DNA；生物源DNA是指有生物组织提取出的DNA，例如皮肤、毛发等。DNA条形码(DNA barcode)是指生物体内能够代表该物种的、标准的、有足够变异的、易扩增且相对较短的DNA片段，已经成为生态学研究的重要工具。
DNA宏条形码(Metabarcoding)技术是一项传统DNA条形码与高通量测序技术相结合后衍生出的研究方法，能够无创地调查来自许多生态系统的物种丰富度。利用Metabarcoding技术对环境样本中的eDNA进行检测，可以获取环境样本中DNA所属物种的分类学信息和基因功能信息，进而研究不同生态环境中物种的多样性，追溯物种来源、发展以及变化，目前在环境生物学研究中得到广泛运用。
底栖动物种类繁多，数量巨大，占据着多种不同的生态位，发挥着不同的生态作用，直接或间接影响不同营养级生物类群的数量及分布。已有大量研究利用底栖动物生物完整性指数(Benthic Index of Biotic Integrity,B-IBI)来评价水生态健康状况。由生态环境部发布的《河流水生态环境质量监测与评价技术指南(征求意见稿)》和《湖库水生态环境质量监测与评价技术指南(征求意见稿)》中也将B-IBI指数纳入淡水生态系统健康评价体系中。但是目前的底栖动物完整性指数构建多基于显微镜下的形态学物种识别数据，监测效率低、耗时耗力，严重依赖于鉴定人的专业知识和经验，且识别物种数目有限，无法反映环境梯度变化下底栖动物遗传组成的变化，用于构建B-IBI指数的候选参数较少，且缺乏代表性，无法满足日益增长的生态环境监测和评价需求。
发明内容
1.要解决的问题
针对目前底栖动物完整性指数构建方法存在监测效率低、耗时耗力，严重依赖于鉴定人的专业知识和经验，且识别物种数目有限，无法反映环境梯度变化下底栖动物遗传组成的变化，用于构建B-IBI指数的候选参数较少，且缺乏代表性的问题，本发明提供了一种基于环境DNA(Environmental DNA，eDNA)宏条形码(Metabarcoding)监测方法构建的底栖动物生物完整性评价方法，用于更加高效，经济和高分辨率地评价水生态健康状况。
2.技术方案
为了解决上述问题，本发明所采用的技术方案如下：
本发明提供了一种底栖动物细胞色素C氧化酶I(COI)区域通用宏条形码扩增和测序引物对，包括一条上游引物和一条下游引物，其中：
所述的上游引物为mlCOIintF：GGWACWGGWTGAACWGTWTAYCCYCC；
所述的下游引物为dgHCO2198：TAAACYTCAGGRTGACCRAARAAYCA。
本发明还提供了一种基于eDNA宏条形码技术的底栖动物生物完整性评价方法，该方法包括如下步骤：
(1)参考点(Reference sites，R)和受损点(Impaired sites，I)的确定；
(2)采集参考点和受损点的样品并提取环境DNA；
(3)以提取的环境DNA为模板，扩增条形码片段、测序并分析所属物种类别；
(4)根据分析结果计算定义的多种生物完整性指数候选指标；
(5)根据步骤(4)的计算结果进行筛选最终的生物完整性指数指标；
(6)根据步骤(5)中确定的最终指标计算底栖动物生物完整性得分；
(7)根据底栖动物生物完整性得分划分底栖动物健康等级。
优选地，步骤(1)中，参考点的选择能够代表监测和评价的水域内，未受人为干扰或所受人为干扰较小的具有最优生物及生境状态的点位，依据水质情况、营养化程度、水生植物覆盖、人为干扰强度确定，其余监测位点为受损点；所述参考点需要满足水质相对较好、富营养化程度相对较轻、水生植被覆盖度较高、人为干扰强度相对较低等条件。
优选地，步骤(1)中，参考点需要满足的条件为：
湖库的具体筛选指标为：(i)总氮、总磷等营养盐指标达到地表水Ⅳ类水标准；(ii)综合营养状态指数小于监测样本的25％分位数；(iii)有沉水植被分布；(iv)样点水域无航道、养殖和娱乐等功能，受水利工程影响小；
河流的具体筛选指标为：(i)总磷≤0.3mg/L，氨氮≤2.0mg/L；(ii)监测断面上游来水方向有湖、荡等前置缓冲区或沿岸有水生植物分布；(iii)监测断面所在河段无航道、养殖和娱乐等功能，受水利工程影响小。
优选地，步骤(2)中，样品采集包括：采集水体底部沉积物，清洗，冷冻保存；环境DNA提取包括解冻，乙醇浸提，离心，提取DNA。
优选地，样品采集包括：用彼得森采泥器、三角拖网等工具采集水体底部沉积物，筛网、过滤、冲洗，基本沥干混合物水分，-40℃到-80℃冷冻保存，冷冻时间为1-3天。
优选地，环境DNA提取采用本领域常规手段进行。更进一步地，环境DNA提取包括：解冻后称重，1g混合物换算为1mL，加入3～5倍体积的无水乙醇，乙醇浸提3～30日后，摇匀，静置后抽取2mL上层液体，经过真空离心后，获得干燥的组织残渣，用试剂盒提取环境DNA，更进一步地，离心条件为8000转/分，离心时长为5分钟。
优选地，步骤(3)中，扩增和测序包括：采用上述COI引物对扩增步骤(2)中的环境DNA，并对PCR产物进行高通量测序，通过qiime、usearch等软件平台处理测序数据，获得每个样本对应的操作分类单元(Operational Taxonomic Units，OTU)信息，并采用NCBI数据库对OTU进行物种注释。
优选地，步骤(4)中，定义的生物完整性指数候选指标需要包括群落丰富度、物种优势度、物种组成、物种耐污值、遗传多样性。
优选地，步骤(4)中，定义的生物完整性指数候选指标不少于27种。
优选地，步骤(4)中，定义的生物完整性指数候选指标及计算公式包括如下：
优选地，步骤(5)中，筛选最终的生物完整性指数指标包括：
a)分布范围筛选：计算每个指数分布的最小值，四分位数(Q25)，中位数，四分位数(Q75)，最大值及值>0的点位占比，保留在80％以上点位的值>0且最大值>Q75>中位数>Q25>最小值的指数，舍弃值域区间相对狭窄的指标；
b)箱体判别分析：采用箱线图法分析进入判别能力分析的各参数在参照状态和受损状态之间的分布情况，即比较参照状态和受损状态25th～75th分位数范围，即箱线图箱体重叠情况(IQ)，保留IQ≥2的指标；
c)相关性分析：对步骤b)中保留的指标进行Pearson相关分析，检验各指标反映信息的独立性，选常用且代表性强的指标作为优先保留指数，筛选掉和优先保留的指数高度相关的指标，所述高度相关是指两个指标间的相关系数|r|大于0.750；优选地，所述优先保留指数包括M01_总分类单元数，M12_颤蚓科_仙女虫科占比，M27_平均遗传多样性，M02_软体动物分类单元或M18_软体动物占比；M24_BMWP指数或M25_ASPT指数或M26_BI指数，因生态学上一般认为物种丰富度越高，生物完整性更高，生态系统越稳定，因此优先保留M01_总分类单元数，其包含的信息量最大，也是目前较为常用的生物指数；颤蚓科和仙女虫科在eDNA生物监测数据中的相对丰度更高，代表性更强，且颤蚓科和仙女虫科是耐污种，对环境具有很好的指示意义，因此M12_颤蚓科_仙女虫科占比可进一步保留；因M27_平均遗传多样性通过序列差异反映每个位点的遗传多样性，遗传多样性更高，分子水平的生物完整性更高，群落对外界压力的抵抗能力更强，可予以保留；软体动物对外界胁迫非常敏感，因此可选择性保留M02_软体动物分类单元数和M18_软体动物占比；M24_BMWP指数，M25_ASPT指数和M26_BI指数是目前较为常用的考虑物种类群耐污值的综合指数，可选择性保留一个。
d)最终候选指标：根据上述步骤a)、b)、c)确定最终候选指标。
优选地，步骤(6)中，底栖动物生物完整性的得分计算包括：
a)指数分值Mi计算：对干扰越强，值越低的指数，以监测样本值的95％分位数为期望值(即生态目标值)，指数分值等于监测值除以期望值；对于干扰越强，值越高的指数，则以5％分位数为期望值，计算方法为：(最大值-监测值)/(最大值-期望值)。若计算结果大于1，按1计，小于0，按0计。
b)B-IBI计算：
B-IBI＝∑Mi
c)B-IBI归一化
为方便以同一尺度进入总体评价，对底栖动物完整性指数进行归一化，其中B-IBIq95为所有监测位点B-IBI的95％分位数：
优选地，步骤(7)中，底栖动物完整性指数健康等级的评价，以监测样本B-IBI归一化值的95％分位数为“优”的健康分级标准，再采用4分法进行分级，标准如下：
3.有益效果
本发明与现有技术相比，其有益效果在于：
(1)本发明提供的一种基于eDNA宏条形码技术的底栖动物生物完整性评价方法，与传统方法相比，eDNA技术通常所需的人力、物力和时间更少；该方法操作易学，无需具有很高的分类鉴定能力；底栖动物中许多个体偏小，且分类特征不明显，eDNA宏条形码监测可识别大量隐蔽物种，灵敏度更高；通过序列差异识别形态学容易忽略混淆的物种，提高物种鉴定的准确率。
(2)本发明提供的一种基于eDNA宏条形码技术的底栖动物生物完整性评价方法，高分辨率地评估不同人类活动干扰强度下水生态健康状况。基于eDNA宏条形码监测可以直接从物种以下水平的分子分类单元变化(如OTU)，并可通过比较不同点位的遗传多样性，来评估底栖动物的完整性，评价结果的分辨率更高。
(3)本发明提供的一种基于eDNA宏条形码技术的底栖动物生物完整性评价方法，该方法可以标准化操作，监测和评价结果可比性强，有助于大尺度时间和空间尺度上的比对。该方法采取标准化的分子生物学和生物信息学流程，不同样品、不同人员的实验结果具有可比性，有助于我国重点流域的生态恢复状况的评估和各大流域间水生态健康状况的比较。
附图说明
图1是实施例1筛选后候选指数箱体分布状况；
图2是实施例1候选指标间的Pearson相关性；
图3是实施例1太湖流域湖库底栖动物完整性指数检验，其中：(a)参考点与受损点的区分度；(b)本方法与基于形态学监测的底栖动物完整性评价的一致性；
图4是实施例2筛选后候选指数箱体分布状况；
图5是实施例2候选指标间的Pearson相关性；
图6是实施例2太湖流域河流底栖动物完整性指数检验，其中：(a)参考点与受损点的区分度(b)本方法与基于形态学监测的底栖动物完整性评价的一致性。
具体实施方式
下面结合具体实施例对本发明进一步进行描述。
实施例1
本实施例样本采集自江苏太湖流域的30个湖泊和水库位点，以下分步骤对实施方法进行具体介绍。
(1)参考点和受损点确定
结合区域水生态特点及环境管理目标，根据水质相对较好、富营养化程度相对较轻、水生植被覆盖度较高、人为干扰强度相对较低的确定条件，根据参照状态的遴选条件，最终选择塘马水库、沙渚、浦庄、椒山、大溪水库南、漫山、东西氿作为太湖流域湖泊和水库(湖库)底栖动物完整性指数构建的参照点，其余点位均视为受损点。
(2)样品采集和环境DNA提取
用彼得森采泥器、三角拖网等工具采集水体底部沉积物，经筛网过滤、冲洗后，基本沥干混合物水分，保存在塑封袋中，-20℃保存。解冻后称重，1g混合物换算为1mL，加入3倍体积的无水乙醇，乙醇浸提7日后，摇匀，稍微静置后抽取2mL上层液体，经过真空离心后，获得干燥的组织残渣，用OMEGA Water试剂盒提取环境DNA。
(3)DNA条形码片段扩增和高通量测序与分析
使用COI引物对(上游引物为mlCOIintF：GGWACWGGWTGAACWGTWTAYCCYCC；下游引物为dgHCO2198：TAAACYTCAGGRTGACCRAARAAYCA)扩增步骤(2)中的环境DNA，并对PCR产物进行高通量测序，通过qiime、usearch等软件平台处理测序数据，获得每个样本对应的OTU信息，并采用NCBI数据库对OTU进行物种注释。
(4)底栖动物生物完整性指数的候选指标及计算
基于eDNA宏条形码监测数据筛选的底栖动物生物完整性指数的候选指标如下表1，包括群落丰富度参数6个，物种优势度参数2个，物种组成参数15个，物种耐污值参数3个，遗传多样性参数1个。
表1底栖动物生物完整性指数的候选指标及计算
(5)底栖动物生物完整性指数的候选指标筛选
a)分布范围筛选：计算每个指数分布的最小值，Q25，中位数，Q75，最大值及值>0的点位占比，保留在80％以上点位的值>0且最大值>Q75>中位数>Q25>最小值的指数，初步分析显示M03_EPT分类单元数、M22_EPT占比在20％以上的点位均为0，被舍弃，其余参数共计25项进入下一步分析；
b)箱体判别分析：采用箱线图法分析进入判别能力分析的各参数在参照状态和受损状态之间的分布情况。即比较参照状态和受损状态25th～75th分位数范围，即箱线图箱体重叠情况(IQ)，只有IQ≥2的参数才作进一步分析。结果显示，共15个候选指数IQ≥2，这些参数的箱体分布状况见图1。
c)相关性分析：对15个生物指标进行Pearson相关分析，检验各指标反映信息的独立性，两个指数间的相关系数|r|大于0.750表示两个指数间高度相关，对高度相关的指标，取其中一个即可，结果如图2所示。优先保留M01_总分类单元数，因此剔除与之高度相关的M24_BMWP指数(r＝0.87)；优先保留M27_平均遗传多样性，因此剔除与之高度相关的M08_前3位优势度(r＝-0.81)、M09_香农多样性指数(r＝0.81)、M10_Pielou均匀度指数(r＝0.78)、M11_Simpson指数(r＝0.76)；优先保留M12_颤蚓科_仙女虫科占比；优先保留M25_ASPT指数，因此剔除与之高度相关的M02_软体动物分类单元数(r＝0.95)；优先保留M18_软体动物占比，因此剔除M13_寡毛类占比(r＝-0.9)、M16_寡毛纲_摇蚊科占比(r＝-0.99)，M19_腹足纲占比与软体动物重复，故剔除。
d)最终候选指标：保留M01_总分类单元数、M07_第一优势度、M12_颤蚓科_仙女虫科占比、M18_软体动物占比、M25_ASPT指数、M27_平均遗传多样性这6个参数共同构成B-IBI。
(6)底栖动物生物完整性计算，划分健康等级
a)B-IBI计算：对干扰越强，值越低的指数，以监测样本值的95％分位数为期望值(即生态目标值)，指数分值等于监测值除以期望值；对于干扰越强，值越高的指数，则以5％分位数为期望值，计算方法为：(最大值-监测值)/(最大值-期望值)。若计算结果大于1，按1计，小于0，按0计，其具体参数计算公式见表2，计算值大于1的取1。
表2太湖流域湖库底栖动物完整性指数参数计算公式
b)归一化
为方便以同一尺度进入总体评价，对底栖动物完整性指数进行归一化，其中B-IBIq95为所有监测位点B-IBI的95％分位数：
c)健康等级
以监测样本B-IBI归一化值的95％分位数为“优”的健康分级标准，再采用4分法进行分级，太湖流域湖库底栖动物完整性指数体系分级标准见表3。
表3太湖流域底栖动物生物完整性指数体系分级标准
(7)底栖动物生物完整性评价结果的可靠性
a)参考点和受损点的区分度
根据构建的太湖流域湖库底栖动物完整性指数体系，计算各点位B-IBI值并绘制箱线图，如图3a，可见参照点与受损点的箱体无重叠，即上述流程构建的B-IBI值可以区分参照点与受损点的底栖动物完整性状态。
b)和形态学监测获得的生物完整性评价结果的一致性
本方法基于环境DNA宏条形码监测获得的底栖动物完整性评价结果和基于形态学监测获得的B-IBI具有显著正相关性(图3b)。
实施例2
样本采集自太湖流域的37个河流和溪流位点。经与实施例1相同的步骤得到37个河流和溪流位点的底栖动物完整性分析结果如下：
(1)根据参照状态的遴选条件，最终选择钓渚大桥、入江水道_戴楼衡阳、旧县、航管站、太浦闸、落蓬湾、赵屯作为太湖流域河流底栖动物完整性指数构建的参照点，其余非参照点位则均为受损点。
(2)经分布范围筛选，所有指数可变范围较宽，对胁迫响应的变化空间较大，可以参与构建完整性指标体系，故所有指数共27项都进入箱体判别分析。
(3)经箱体判别分析，结果显示，共13个候选指数IQ≥2，部分参数的箱体分布状况见图4。
(4)13个生物指标经相关性分析，如图5所示，M09_香农多样性指数作为多样性指数之一，应当纳入评价体系，予以保留，因此剔除与之高度相关的M08_前3位优势度(r＝-0.96)；指示物种方面，保留M12_颤蚓科_仙女虫科占比和M18_软体动物占比(r＝-0.78)，因此剔除M13_寡毛类占比(r＝0.83)、M16_寡毛纲_摇蚊科占比(r＝0.81)；由于M20_瓣鳃纲占比和M21_蚬属占比与M18指数略有重复，因此剔除。
(5)最终保留M01_总分类单元数、M09_香农多样性指数、M12_颤蚓科_仙女虫科占比、M18_软体动物占比(r＝-0.78)这4个参数共同构成B-IBI。B-IBI各参数计算公式如下表4，计算值大于1取1。
表4太湖流域河流底栖动物完整性指数参数计算公式
(6)为方便以同一尺度进入总体评价，对底栖动物完整性指数进行归一化，其中B-IBIq95为所有监测位点B-IBI的95％分位数：
(7)经归一化，太湖流域河流底栖动物完整性指数体系分级标准同表3。
(8)根据构建的太湖流域河流底栖动物完整指数体系，计算个点位B-IBI并绘制箱型图如图6a，与基于形态学监测获得的B-IBI相关性结果如图6b，结果均表明，基于上述流程构建的B-IBI值适用于太湖流域河流底栖动物完整性评价。
