徐 念, 熊美华, 邵 科, 阙延福, 李键庸
水利部中国科学院水工程生态研究所, 水利部水工程生态效应与生态修复重点实验室&湖北省水生态保护与修复工程技术研究中心, 湖北 武汉 430079
淡水生态系统为人类文明发展提供了重要资源,然而面临着比陆地和海洋生态系统更严重的威胁[1-2],是全球生物多样性丧失最严重的区域之一,却缺乏足够的调查研究[3-5]. 物种是生物多样性的核心组成部分,对物种资源的调查和监测是生物多样性保护的重要基础[6-7],在全球生物多样性丧失日益严重的背景下,对生物多样性现状全面了解的需求尤为紧迫[8].
长江是我国淡水资源的宝库,但多年来受人类活动影响,长江生物多样性持续下降,水生态保护形势严峻. 2018年4月,生态环境部、农业农村部和水利部联合发布了《重点流域水生生物多样性保护方案》,明确提出开展长江等重点流域水生生物多样性本底调查的主要目标,以及在流域干流等水域开展鱼类及水生哺乳动物等物种的组成、分布和种群数量进行调查观测的重点任务. 同年9月,国务院办公厅发布了《关于加强长江水生生物保护工作的意见》,要求到2020年长江流域重点水域实现常年禁捕,强调提升监测能力,全面开展水生生物资源与环境本底调查,提高监测系统自动化智能化,加强大数据集成分析. 生物多样性研究是一项工作量巨大的长期任务,过去零散的调查数据已无法满足当前的研究和管理需求[7],需要在流域范围内获取更全面更系统的资料,因此对调查方法的科学性和有效性、调查结果的全面性和准确性,以及对大数据的综合分析能力提出了更高的要求.
环境DNA宏条形码(environmental DNA metabarcoding)指从环境样本中提取DNA,利用高通量测序分析获得大量DNA序列,通过序列检索比对检测环境中的多个物种[8],因其简单高效、灵敏度高等优势获得广泛关注[9]. 环境DNA宏条形码对物种及其丰度的检测在微生物中已形成标准化的流程,如16S和18S扩增子多样性分析,有力地补充了传统培养方法的缺失[10]. 进入21世纪初以来,通过环境DNA检测大型生物的技术逐渐兴起,水样环境DNA被用于鱼类和两栖类的特异性物种检测[11-12]. Thomsen等[13]于2012年首次报道了利用水样环境DNA高通量测序进行人工池塘两栖类和鱼类多物种检测,此后环境DNA宏条形码成为水生态系统生物多样性研究的热点,在海洋、湖泊、河流等生境得到应用[14-16]. 国内研究中,单一物种的环境DNA检测已针对鲢[17]、中华鲟[18]和长江江豚[19-20]在长江流域得到应用,已有通过水样环境DNA克隆测序进行的少数鱼类物种检测初步研究[21],水样环境DNA宏条形码研究已在海洋环境[22]开展,淡水生态系统环境DNA宏条形码研究在国内引发关注但鲜见报道.
长江中下游具有独特的江湖生态系统,分布有多种重要经济鱼类[23],也是濒危水生动物的重要栖息地[24],而长江中下游因流量大、泥沙含量高、物种资源衰退等因素对环境DNA检测提出了巨大挑战[20-21]. 该研究在长江中下游3个江段(新滩、安庆和芜湖)采集水样,利用通用分子标记建立长江水样环境DNA宏条形码物种检测体系,分析该方法在长江水生物种种类组成和资源量评估中的有效性,旨在探索高敏感度非入侵式的长江生物多样性监测新方法,以期为长江水生态监测体系建设提供技术支撑.
于2016年1月在长江中下游干流新滩、安庆、芜湖3个江段各设置1个采样点,每个采样点用全新2 L可密封广口瓶(Nikko,日本亚速旺公司)采集3个2 L表层水样,用0.45 μm聚醚砜滤膜(Pall,美国颇尔公司)真空抽滤,玻璃抽滤漏斗在每次使用前用10%次氯酸钠消毒液浸泡30 min并充分冲洗干净. 抽滤完成后,滤膜冷冻保存并尽快送回实验室. 用强力水样DNA提取试剂盒(Mobio,美国Mobio实验室)提取滤膜DNA,DNA溶液于-20 ℃下保存.
利用线粒体细胞色素B简并引物L14912-CYB (5′-TTCCTAGCCATACAYTAYAC-3′; Y=C或T)和H15149-CYB (5′-GGTGGCKCCTCAGAAGGACATTTG KCCYCA-3′; K=G或T, Y=C或T)[25]对环境DNA样本进行PCR扩增,扩增产物片段长度在285 bp左右. 该引物的目标位点是在脊椎动物中广泛存在的保守区域,用于河流水样环境DNA鱼类物种克隆测序检测[21,26]. 总体积50 μL的反应体系包括4 μL的10×PCR Buffer、1 μL的dNTP (10 mmol/L)、1 μL的正向引物(10 pmol/μL)、1 μL的反向引物(10 pmol/μL)、5 μL的DNA模板和36.5 μL的灭菌双蒸水. 反应条件:94 ℃预变性5 min;94 ℃变性1 min、50 ℃ 退火1 min、72 ℃延伸1.5 min (45个循环);72 ℃最后延伸7 min.
采用PCR产物进行文库构建,用Illumina Miseq平台进行高通量测序(深圳华大基因科技服务有限公司),下机数据经过过滤,滤除低质量的读长序列,得到有效序列. 使用软件FLASH v1.2.11进行序列拼接,利用重叠关系将双末端测序得到的成对有效序列组装成一条拼接序列,去除没有重叠关系的序列. 利用软件USEARCH v7.0.1090将拼接好的序列聚类为OTU(operational taxonomic unit,操作分类单元),获得OTU代表序列.
该研究利用美吉生物云生物信息分析云平台(www.majorbio.com)对测序序列进行生物信息学分析. 首先将OTU代表序列在Genbank的核酸序列数据库中进行物种注释,同时利用序列相似性在线检索工具BLAST进行人工校核,除去两端引物序列,核心目标片段100%匹配则被认定为该物种DNA. 然后利用分析平台的生物多样性云分析流程v3.0进行交互分析,对注释OTU进行筛选,统计不同分类水平的物种数量、OTU数量、序列数量,计算物种序列相对丰度.
2.1.1种类组成
长江中下游干流新滩、安庆、芜湖3个江段水样环境DNA宏条形码检测共得到有效拼接序列 334 623 条,平均序列长度为285 bp,聚类得到425个OTU,其中在数据库中达到匹配的OTU为54个(总序列数 162 398),共注释10目13科32种(见表1),其中鱼类20种、水生哺乳动物1种、鸟类4种、陆生哺乳动物7种. 鱼类包括鲤形目、鲇形目、鲈形目和鲱形目物种,其中鲤形目物种数最多(见图1);水生哺乳动物为国家一级保护动物长江江豚,是长江中唯一现存的水生哺乳动物,也是目前中国唯一的淡水豚类;鸟类包括家禽和野生物种;陆生哺乳动物包括偶蹄目牲畜、褐家鼠和人.
所有注释物种中,水生物种数占65.6%,其他为陆生动物. 以无脊椎动物为主要对象的环境DNA宏条形码研究报道了河流水样环境DNA包含水生和陆生生物群落的多样性信息[27]. 该研究以脊椎动物为主要对象,结果表明,长江水样环境DNA包含水生和陆生物种信息,通过环境DNA宏条形码可检测多类群物种. 该研究采样点位于人口密集的城镇江段,水样中检测到的陆生物种大多与人类生产生活息息相关,水样环境DNA物种组成体现出人类活动与河流生态系统的密切联系.
2.1.2序列相对丰度
在所有注释物种中,鱼类序列数占总序列数的78.5%,其中鲤形目序列数占鱼类总序列的96.2%,鲱形目占3.5%,鲇形目占0.2%,鲈形目占0.1%,各类目序列数见表2. 在20种鱼类物种中,序列数排序前5位的依次为鲢、团头鲂、青鱼、短颌鲚和草鱼,各物种序列数在鱼类中占比依次分别为52.5%、23.7%、17.0%、3.5%和1.3%. 水生哺乳动物长江江豚的序列数在所有注释水生物种中占比约为 1/4 000. 陆生物种中,序列相对丰度最高的是鸡,其序列数占陆生物种总序列数的96.0%,其次为人(占比为2.6%)、鸿雁(占比为0.8%)、牛(占比为0.6%)和绵羊(占比为0.2%). 综合种类数量和序列相对丰度,长江中下游水样环境DNA主要包含水生物种信息,鱼类是水生脊椎动物的主要类群,其中鲤形目占绝大多数.
序列相对丰度受物种丰度的影响,在一定程度上序列相对丰度体现了物种在环境中可能具有的丰度相互关系[28],但二者的相关性不太确定[13,29]. 从环境中物种DNA的释放、运输和降解,到环境样本的采集、目的片段的扩增和测序,检测结果中序列相对丰度受多种因素的影响,如环境DNA的来源、浓度分布、引物的扩增偏向性等[30]. 环境样本中特定物种的环境DNA检测浓度也不一定与其物种丰度具有显著相关关系[31]. 尽管在多变的自然环境中相关性较弱,物种环境DNA浓度依然体现了物种丰度的实质性变化,利用环境DNA宏条形码准确评估物种数量或生物量则需要对采样策略和检测分析技术提出更高的要求[15,32].
表1 长江中下游3个江段水样环境DNA宏条形码注释物种
注: +表示采样点水样中检测到物种环境DNA.
在OTU水平基于Bray-curtis距离算法的样本进行层级聚类(hierarchical clustering),根据不同的距离阈值可将样本划分为聚类小组,结果如图2所示. 由图2可见,同一采样点的样本首先聚类在一起,然后安庆采样点和芜湖采样点的样本聚类为一支,新滩采样点样本则单独聚类为一支. 相似性分析结果显示,新滩、安庆、芜湖3个采样点组间差异显著大于其组内差异(见图3)(P=0.001),置换多因素方差分析结果表明,分组对样本差异具有可信的解释度(R2=0.924 6,P=0.005). 综上,同一采样点样本差异度较小,表明在一定范围内环境DNA的分布与组成具有稳定性,而在河流中不同采样点间的显著差异可以为空间差异性研究提供基础. 河道间距较近的安庆和芜湖采样点聚类距离也较近,推测河流水样环境DNA的组成相似性可能与采样点间的距离有关. 不同江段物种组成的差异性会导致水样中环境DNA的组成差异,但在该研究中检测结果不排除其他因素的影响,环境DNA组成在河流中的空间差异性还需更多研究数据来证实.
图1 长江中下游水样环境DNA宏条形码检测种类组成Fig.1 The number of species detected from water samples in the middle and lower reaches of the Yangtze River using environmental DNA metabarcoding
表2 长江中下游水样环境DNA宏条形码检测物种注释序列数Table 2 The number of sequences for each order identified by environmental DNA metabarcoding of water samples in the middle and lower reaches of the Yangtze River
序号物种分类序列数序号物种分类序列数1鲤形目122 6326雁形目2792鸡形目33 2267鲇形目2763鲱形目4 4468鲈形目1694灵长目9239啮齿目415偶蹄目37410鲸目32
注: 1、2、3表示采样点样本序号. 下同.
图2 基于OTU水平的样本层级聚类树
Fig.2 Hierarchical clustering tree on OTU level
图3 基于OTU水平的组间距离盒型图Fig.3 Distance box plot on OTU level of each sample groups
2.3.1鱼类
该研究通过水样环境DNA检测到的20种鱼类均为在长江中下游传统鱼类资源调查中出现过的物种[33-35],包括构成重要渔业资源的经济鱼类草鱼、青鱼、鲢、鳙、鲤、鲫等. 环境DNA检测到的4个类目是长江中下游渔获物中的主要类目,在所有20种鱼类中,鲤形目种类数占总种类总数的60%,鲇形目占25%,鲈形目占10%,鲱形目占5%. 长江鱼类资源的特点是鲤科鱼类种类多[23,36],该次环境DNA检测结果也显示,鲤科种类占多数. 在已报道的长江中下游干流渔获物调查中,鲤科种类占52%~61%[33-35,37],该次环境DNA检测鲤科种类占鱼类种类总数的50%~60%. 在上述渔获物调查中,鲇形目、鲈形目和鲱形目分别占11%~22%、12%~19%和3%~6%,该次环境DNA检测结果显示,其分别占18%~23%、0~12%和0~6%. 水利部中国科学院水工程生态研究所于2017年3月在芜湖江段开展了6 d渔获物调查,记录每日渔民单船渔获物种类,各类目鱼类种类数从多到少依次为鲤形目、鲇形目、鲈形目和鲱形目. 图4显示了安庆和芜湖两江段的渔获物调查结果和该次环境DNA调查结果的各类目种类数量占比,从鱼类种类组成来看,环境DNA检测结果体现了渔获物的主要类目,且主要注释类目的种数比例与渔获物组成相似.
注: 采样次数=采样点数×采样天数;渔获物数据来源于文献[34-35,37]和调查记录. eDNA表示环境DNA.
图4 长江中下游两江段渔获物调查和环境DNA宏条形码检测中目水平鱼类种数占比
Fig.4 Proportion of fish species of each phylum investigated in traditional survey and environmental DNA metabarcoding
鲤形目鱼类不仅在种类数上占优势,在DNA序列数上也占据绝对优势,而鲤形目在长江中下游渔获物调查中的尾数和重量均居首位. 对于鲤形目、鲇形目、鲈形目和鲱形目4个主要类目,安庆江段4个类目资源量总和占渔获物总量的99%以上[33,37],芜湖江段4个类目尾数总和占渔获物总量的99%以上[34]. 但环境DNA序列相对丰度和渔获物资源量百分比之间存在较大差异,如鲤形目序列在3个江段中均占90%以上,而在渔获物中鲤形目资源量通常为渔获物总量的60%~80%. 以安庆江段为例,根据2003年3月—2010年3月[33]和2015年4—6月[37]的渔获物调查结果,发现鲤形目尾数和质量分别占渔获物总量的47%~72%和60%~77%,鲇形目尾数和质量分别占17%~46%和12%~32%,鲈形目尾数和质量分别占3%~12%和6%~11%,鲱形目尾数和质量分别占>1%~6%和>1%~4%. 该研究环境DNA检测安庆江段鲤形目序列占鱼类总序列的92%,鲱形目占7.7%,鲇形目和鲈形目均不足0.1%,与渔获物反映的各类目资源量相差较大. 一方面如2.1.2节所述序列相对丰度和物种资源量之间并不具有稳定的对应关系; 另一方面该研究每个江段单次调查获取的数据有限,不能满足准确的资源量统计分析需求,但序列相对丰度体现出了其中的优势类群鲤形目鱼类.
从调查取样数量来看,对比安庆和芜湖两江段的已有报道[33-35,37],将同一采样点同一天内的采样记为1次采样,该研究每个江段仅进行了1次采样,水样的采集通常在几分钟内完成,而渔获物报道中鱼类标本的采集在每个江段采样次数(采样点数×采样天数)为20~336次(部分数据见图4),每个采样点的总采样时间(单次采样时间×采样次数)最低为10 h,最高达 8 064 h,且多数调查在取样时同时使用多个网具,实际工作量更大. 从单次调查种类数量来看,在水利部中国科学院水工程生态研究所于2017年3月开展的芜湖江段6 d渔获物调查中,每天调查种类为9~21种,平均单次渔获物调查种类数为15.7种,网具数量4网,该研究芜湖江段仅一次环境DNA调查种类为17种,平行样本3个. 总的来说,该次环境DNA调查检测到的鱼类种数为各江段渔获物种数的31%~49%,而采样时间不足努力量最少渔获物调查的1%,是努力量最多渔获物调查的十万分之一. 可见,与传统渔获物调查相比,环境DNA调查方法在采样效率方面具有明显的优势. 对比长江中下游安庆江段4个不同年代的渔获物调查报道[33,35,37-38]中出现的物种,环境DNA调查安庆江段检测到17种鱼类,其中2种鱼类(团头鲂、红鳍原鲌)未在渔获物调查中出现,2种鱼类(达氏鲌、寡鳞飘鱼)仅在1次渔获物报道中出现,仅7种鱼类在4次渔获物报道中均有出现,渔获物调查时间跨度从短至长依次分别为3个月、12个月、6 a和8 a,而环境DNA取样仅用1 d完成,可见长时间渔获物调查结果仍然会受到捕捞方式的限制,甚至造成一些常见的经济鱼类(如草鱼、青鱼等)在渔获物调查中未采集到[35],环境DNA调查通过简单的取样方式,可以对传统调查结果进行补充. 在国际同类研究中,利用环境DNA宏条形码技术在淡水生态系统的流水生境中检测到的鱼类种数分别为8种[39]、12种[40]、16种[41]和19种[42],这些研究均表明环境DNA宏条形码和传统调查方法相结合能调查到更多的鱼类物种,且环境DNA调查方法比传统调查方法效率更高. 对于鱼类资源调查来说,环境DNA宏条形码需优化采样方案和检测体系,如增加采样点和采样次数、使用更高效的分子标记等,以检测到更多物种,并对其进行资源量评估.
注: 传统调查分布密度来源于农业农村部发布的2017年长江江豚科学考察结果[43].
图5 长江中下游干流3个采样点位置和长江江豚环境DNA序列相对丰度
Fig.5 Location of the three sampling sites in the middle and lower reaches of the Yangtze River mainstream and relative abundance of environmental DNA sequences of Neophocaena asiaeorientalis asiaeorientalis
2.3.2长江江豚
安庆采样点3个样本,以及新滩采样点1个样本,均检测到长江江豚DNA序列,且安庆采样点长江江豚DNA序列相对丰度(物种序列数在总测序序列数中占比)约是新滩采样点的36倍,芜湖采样点未检测到江豚DNA序列(见图5). 针对单一物种的环境DNA检测技术对濒危物种具有高敏感度[13],其为长江江豚监测提供了有力工具[19-20]. 该研究结果显示,以多物种为对象的水样环境DNA宏条形码对于濒危物种也具有一定的敏感度.
综合2006、2012和2017年3次长江江豚考察结果[43-45],发现安庆江段长期以来均处于长江干流长江江豚密度最高的江段,而中下游鄂州以上和安庆以下江段则密度较低. 该研究中安庆采样点位于皖河口下游约2 km处,皖河口是安庆市江豚自然保护区的重点水域,是长江江豚群体的重要栖息地[35],安庆采样点的长江江豚环境DNA检出率和DNA序列相对丰度在3个采样点中都是最高的,与传统调查中长江江豚密度最高的江段相一致. 与2017年江豚科学考察结果相比,新滩采样点处于长江江豚分布密度居中的江段,位于长江新螺段白鱀豚国家级自然保护区内靠近下游边界,该保护区也是长江江豚的聚集地;而芜湖采样点处于江豚分布密度最低的江段,距离上游铜陵淡水豚国家级自然保护区下游边界约60 km[46]. 该次环境DNA检测结果表明,长江中下游干流长江江豚环境DNA检出率和序列相对丰度与其密度分布相符合,但利用环境DNA宏条形码进行长江江豚群体分布和资源量的监测还需要针对性的开展进一步研究.
该研究体现了环境DNA宏条形码技术在长江水生态系统生物多样性研究领域的主要优势,采集水样将调查研究过程对目标物种和生态环境的影响降到极低,高效的采样需要较少的人力和时间付出,且宏条形码分析流程标准化,有利于不同样本间的比较,结果可体现大尺度区域群落结构,突出其中的优势类群,对濒危物种亦具有一定敏感度. 此外,高通量测序可以在同一样本中快速获得比克隆测序更多的序列作为物种注释的基础,在已报道的长江鱼类环境DNA研究[21]中,通过水样环境DNA克隆测序方法在24个采样点共检测到10种鱼类序列,每个采样点检测到1~3种鱼类,而该研究通过高通量测序在3个采样点共检测到20种鱼类,每个采样点检测到13~17种鱼类,比克隆测序方法效率大幅提高. 随着高通量测序技术和生物信息学分析方法的快速发展,环境DNA宏条形码将揭示更多水体中所包含的生物多样性信息. 该研究中环境DNA宏条形码方法在物种多样性检测方面也存在一定的局限性. 在所有测序得到的OTU中,成功注释物种的OTU比例低,约有87.3%的OTU无法在数据库中匹配到物种,这些未注释的OTU包含总测序序列数量的51.5%,数据库中收录的DNA序列成为环境DNA宏条形码物种注释的重要限制因素[29]. 环境DNA宏条形码技术极大程度上依赖于检测中所选择的分子标记,不同的分子标记可以针对不同的目标类群,同时分子标记也可能存在一定的扩增偏向性[30],该研究检测到的物种有限,后续研究可筛选更多通用分子标记,以获取足够的序列信息用于生物多样性研究. 在该研究中各物种序列相对丰度在不同采样点各不相同,同一采样点的样本则相当接近,表明水样环境DNA宏条形码检测在一定范围内具有稳定性和可重复性. 测序结果中存在某些物种的序列数量占比极大的现象,如在新滩采样点,鲢的序列数占总注释序列数的99.8%,而新滩采样点注释的物种总数最少,约为其他2个采样点的62%,安庆和芜湖采样点物种数量则较为接近,其原因可能是,新滩采样点水样中鲢的DNA浓度较高,在有限的测序数量下相对丰度较低的物种无法检出. 根据研究对象设计合理的采样方案、设置足够数量的采样点有助于获得更全面准确的检测结果.
近年来,利用环境DNA宏条形码技术针对淡水生态系统脊椎动物多样性的研究日益发展,但其中环境样本来源于大型河流的报道较少,尤其是长江中下游干流这种流量巨大的流水生境,对水样环境DNA物种检测来说存在相当的挑战. 河流是生物多样性信息传递的纽带[47],水样中包含的水陆物种信息可以用于景观尺度的生物多样性分析,以及人类活动对河流的影响研究. 该研究验证了在长江水生态系统中通过水样环境DNA宏条形码技术进行生物多样性检测分析的可行性,同时体现了通过对河流环境DNA检测进行景观尺度的水陆复合生态系统生物多样性综合分析的可能性. 技术方法的发展和环保理念的革新是未来生态保护不可阻挡的趋势,通过开发高效分子标记、发展测序分析技术,未来环境DNA所包含的物种条形码还将被进一步挖掘,环境样本中展现出的大尺度复合生态系统的丰富信息将为生物多样性研究开拓更广阔的发展前景.
a) 长江中下游干流新滩、安庆和芜湖3个江段的水样环境DNA宏条形码分析共检测到32个脊椎动物物种,其中包括20种鱼类,1种水生哺乳动物(长江江豚)和11种陆生动物,其中鱼类序列数占总注释序列数的78.5%,表明长江水样环境DNA主要包含水生物种信息,同时承载了陆生物种信息,通过水样环境DNA宏条形码可检测不同类群物种.
b) 环境DNA检测到的鱼类中,鲤形目种类占60%,鲇形目占25%,鲈形目占10%,鲱形目占5%,为长江中下游渔获物调查的主要类目,且各类目包含的种类数量比例与渔获物调查结果具有相似性. 在渔获物中尾数和质量均居首位的鲤形目在环境DNA调查中序列数也最多,但4个类目序列相对丰度与渔获物种资源量组成差异较大.
c) 环境DNA调查用约占传统渔获物调查几十至几百分之一的调查次数,百分之一至十万分之一的调查时间,检测到了31%~49%的鱼类种数,表明环境DNA调查效率更高,且环境DNA检测到了在渔获物调查中未出现的物种,可以对传统调查结果进行补充,将二者结合可获得更全面的调查结果.
d) 在新滩和安庆采样点均检测到濒危物种长江江豚环境DNA,表明环境DNA宏条形码对濒危物种具有一定敏感度. 在中下游长江江豚密度最高的安庆江段,长江江豚环境DNA检出率和序列相对丰度在3个采样点中都是最高的,而位于长江江豚密度最低的中下游下段的芜湖采样点,未检测到长江江豚环境DNA. 长江江豚环境DNA检出率和序列相对丰度可能主要受其分布密度的影响.
e) 长江中下游水样环境DNA宏条形码分析结果体现了调查区域的部分群落结构,显示出了优势类群,注释物种包含濒危物种. 环境DNA宏条形码在淡水生物多样性研究中具有高效率和无损伤的优势. 核酸序列数据库和分子标记体系是环境DNA宏条形码物种检测的主要限制因素,进行全面的物种多样性监测和资源量估算需要对采样方案和分析方法进行科学设计. 通过水样环境DNA,可进一步对河流水陆复合生态系统的生物多样性信息进行挖掘.