环境DNA揭示东江流域水生昆虫生物多样性变化及影响因素

摘 要：以东江流域为例，通过环境DNA技术探究了水生昆虫群落的空间分布规律及其环境影响因子。结果显示：1）水生昆虫物种丰富度由高到低顺序为：双翅目（7科）＞蜉蝣目（6科）＞毛翅目（2科）=蜻蜓目（2科）=鞘翅目（2科）＞鳞翅目（1科）。双翅目（84.99%）和蜉蝣目（14.89%）共占总reads数99.88%，摇蚊科（69.20%）、四节蜉科（14.36%）及蚊科（14.26%）三个科合计占总reads数97.82%。2）物种丰富度指数、香农维纳指数、辛普森指数与功能丰富度指数（FRic）、功能离散度指数（FDiv）均与距河口距离呈正相关关系。功能均匀度指数与距河口距离呈负相关关系。3）海拔是影响物种丰富度最重要的环境因子（P＜0.01），COD是影响物种丰富度的重要环境因子（P＜0.05）。海拔、DO、COD是影响功能丰富度的重要环境因子（P＜0.05）。

关键词：东江流域；水生昆虫；物种多样性；功能多样性；环境DNA

中图分类号：S941.419

水生昆虫是河流生物多样性中重要生物类群［1］。水生昆虫对环境变化和人类干扰敏感，对干扰的反应范围广，我们可以通过其被干扰后的影响进行空间分析［2］；水生昆虫有较长的生活周期，能综合反映较长时间段内的水体质量状况［3］。水生昆虫的生态学性状（如摄食方式、大小、生存位置、栖息地改造、移动性等）在很大程度上影响其生态功能贡献［4-5］。盛天进等［6］研究认为物种多样性有助于确定未来河流群落恢复差距。物种多样性虽被广泛使用，但其仅考虑物种的丰度和多度，忽视了物种的功能性状［7］，近年来功能多样性概念得到重视，其是通过量化一个群落内的性状变异或多变量性状差异来理解群落和生态系统［8］。

环境DNA（Environmental DNA， eDNA）是指从土壤、水、冰川或沉积物等环境样品中提取的DNA，可揭示过去和现今的生物多样性信息［9］。与传统采样方法相比，该方法成本更低，样本采集受气候条件影响较小。近年使用环境DNA技术进行生物监测逐渐兴起。环境DNA最早应用于法国池塘中检测到北美牛蛙［10］。随着该方法研究的爆发式增长，不同样本类型的环境DNA在物种丰富度估算方面十分成功［9］。样本类型从淡水［11-13］、海洋［14-15］、陆地［16-18］到空气［19-20］。环境DNA技术的出现为水生生态系统生物多样性的研究开辟了一条新道路，具有极大的应用潜能。

目前，使用环境DNA技术对我国亚热带河流水生昆虫物种多样性及功能多样性变化研究仍处于空白。本研究选取我国典型亚热带河流东江水系为研究区域，旨在探究东江流域水生昆虫的群落多样性三个关键问题：一是东江流域水生昆虫的群落组成特征；二是流域内水生昆虫的空间格局分布规律；三是影响水生昆虫群落结构差异以及物种多样性的因素。

1 材料与方法

1.1 采样点的选择和样品的采集

东江为珠江流域三大水系之一，介于东经114°47′36″至115°52′36″，北纬24°20′30″至25°12′18″之间，位于珠江三角洲的东北端，地势东北高，西南低。干流全长 562 km，发源于江西省寻乌县桠髻钵山，流入广东省境内经龙川、东源、源城、紫金、惠阳、惠城、博罗至东莞市，最后注入狮子洋，经虎门出海，河道平均坡降0.388‰，流域面积35 340 km2。属亚热带季风湿润气候区［21］，河流年均径流量257亿m3。本研究采样于2021年10月，共计51个样点，其中干流23点位、支流浰江6点位、石马河6点位、西枝江16点位，监测点位（见封三图1）示。

环境DNA样品使用Rocker过滤装置、0.45 μm混合纤维素脂滤膜过滤采集的水样，每500 mL过滤一个滤膜，使用灭菌镊子从边缘夹起滤膜，将其放置于5 mL无菌离心管内，-20 ℃干冰保存，带回实验室后保存于-20 ℃冰箱，直到提取DNA。过滤用的杯子在每次使用后均用纯净水冲洗2～3遍。

现场使用多参数水质监测仪（YSI Incorporated，Yellow Springs，USA）进行水质原位测定，包括pH、溶解氧（DO，mg/L）、电导率（EC，μS/cm）、水温（Temp，℃）、溶解性总固体（TDS，mg/L）、氧化还原电位（ORP，mV）、盐度（SAL，PSU）、叶绿素a（CHLA，μg/L）。另采集1 L水样送到实验室，按照中国国家标准测定方法测定：氨氮（NH+4-N，mg/L）、总氮（TN，mg/L）、总磷（TP，mg/L）、磷酸盐（PHOS，mg/L）、高锰酸盐指数（COD，mg/L）、硝酸盐（NO-3-N，mg/L）、亚硝酸盐（NO-2-N，mg/L）等指标。现场还用手持GPS进行经纬度定位和海拔高度测量。

1.2 DNA提取和PCR扩增

DNA提取步骤按照ZymoBIOMICSTM DNA Miniprep Kit说明书进行。正向引物fwhF2（5'-GGDACWGGWTGAACWGTWTAYCCHCC-3'），反向引物EPTDr2n（5'-AAACAAATARDGGTATTCGDTY-3'），采用上述eDNA模板进行水生昆虫的PCR扩增。引物fwhF2-EPTDr2n扩增长度在150～315 bp之间。25 μL反应体系包括：4 μL模板DNA、正反向引物各1 μL、2 μL dNTPs、2.5 μL Buffer、0.125 μL Taq HS以及14.375 μL DEPC。PCR反应程序为：95 ℃预变性30 s; 35个循环包括： 95 ℃变性5 s， 54 ℃退火30 s， 72 ℃延伸10 s; 72 ℃最后延伸5 min。51个采样点各设置一个阴性对照，PCR产物使用2%琼脂凝胶电泳检测。COI区域所有隐性对照均无目的条带，表明无外源污染。

1.3 高通量测序和数据处理

将等量PCR产物汇集在一起，使用Illumina PE150平台进行高通量测序。对测序得到的序列信息进行拼接得到原始数据，经质控和去嵌合体得到可用于后续分析的有效数据。对有效数据进行去冗余、去单一序列和97%相似度的OTU（Operational Taxonomics Table，可操作分类单元）聚类分析，根据OTU聚类结果，对每个OTU的代表序列数据库比对作物种注释，得到对应的物种信息和基于物种的丰度分布情况，即OTU Distribution Table（以下简称 OTU Table）。

1.4 生物信息学分析

采用图图云平台（www.cloudtutu.com）的alpha多样性指数工具包计算α多样性指数，散点图工具包中的spearman方法完成散点图的制作。利用参考文献[22-23]进行功能性状的等级划分。本研究选择水生昆虫的游泳能力、吸附能力、形状、生活型、营养习性、最大体长以及生活史等7个功能性状的27个等级性状纳入分析（见表1）。大多数性状信息基于百度百科（https：//baike.baidu.com）及书籍［24-29］进行搜集。采用Rstudio软件FD功能包的dbFD函数计算功能多样性指数。使用SPSS 27.0.1对环境变量进行线性回归分析。保留ASL、pH、DO、ORP、SAL、CHLA、PHOS、COD、NO-2N等9个环境变量构建简约模型解释群落结构变化。采用凌波微课（http：//www.cloud.biomicroclass.com）DCA分析工具包先做DCA分析。本研究DCA分析后第一轴值为1.39，故选择RDA分析。RDA分析前，为避免环境因子间的共线性，对环境因子数据进行了log（x+1）转换。使用凌波微课RDA分析工具包对环境因子进行可视化。Mantel test分析对α多样性指数与功能多样性指数的相关关系进行检验并在图图云平台进行可视化。物种多样性、功能多样性与环境变量之间关系的网络图与显著影响的环境变量相关性分析图叠加热图的Mantel Test分析以及偏曼特尔检验均使用图图云平台计算制作。采样点图用Arcgis10.7软件制作。

2 结果

2.1 东江水生昆虫物种组成

东江水生昆虫共获得57 142 012条reads，经过生物信息学过滤，保留了6 435 709个reads。在所有采样点和重复中，共获得昆虫纲6目20科89属166种。其中，双翅目和蜉蝣目为两大优势类群，共占总reads数99.88%，优势十分明显。摇蚊科（69.20%）、四节蜉科（14.36%）及蚊科（14.26%）三个科合计占总reads数97.82%，东江流域水生昆虫分类情况见表2。

2.2 物种多样性与功能多样性空间分布

对物种多样性空间格局统计分布（见封三图2-a、2-b、2-c），结果表明：物种丰富度指数范围在3～15之间；香农指数范围在0.03～1.75之间；辛普森指数范围在0.01～0.79之间。对东江流域水生昆虫物种多样性分析后发现，物种丰富度、香农指数、辛普森指数均呈现干流上游及支流高，干流中下游低的分布现象。

对功能多样性空间格局统计分布，见图2-d、2-e、2-f，结果表明：功能丰富度指数范围在0.63～12.06之间；功能均匀度指数范围为0～0.88；功能离散度指数范围在0.52～0.99之间。对东江流域水生昆虫功能多样性分析后发现，功能丰富度指数在干流上游及支流分布明显高于干流中下游，见图2-d；功能均匀度指数在干流上下游及支流分布较为均匀，见图2-e；功能离散度指数较高，且均匀分布在各个点位，见图2-f。

环境DNA还揭示了东江流域水生昆虫物种多样性及功能多样性随距河口距离变化的趋势。结果显示，物种丰富度指数（见图3-a）、Shannon指数（见图3-b）、Simpson指数（见图3-c）与FRic指数（见图3-d）、FDiv指数（见图3-f）均与距河口距离呈正相关关系。FEve指数（见图3-e）与距河口距离呈负相关关系。

2.3 水生昆虫生物多样性变化确定因素

通过方差膨胀因子去除共线性较高（VIF值＞10）的环境因子EC、Temp、TDS、NH+4-N、TN、TP、NO-3-N等引起东江流域水生昆虫群落结构差异的影响因素，保留ASL、pH、DO、ORP、SAL、CHLA、PHOS、COD、NO-2-N等9个环境变量。冗余分析（RDA分析）结果显示，东江流域水生昆虫群落结构同时受ASL、COD（P＜0.01）和PHOS（P＜0.05），见图4。不同环境因子对各轴的贡献也不同。NO-2-N、COD、pH、CHLA、ASL与轴1相关，PHOS、SAL、DO、ORP与轴2相关。其中，NO-2-N、COD与轴1正相关，PHOS、SAL与轴2正相关。第一轴和第二轴分别可以解释采样点与环境关系的78.76%和9.94%。

跟着图5，沿着轴2正方向，PHOS、SAL逐渐增大、DO、ORP逐渐减小，蛾蚋科、摇蚊科、蚊科、螟蛾科、四节蜉科、水蝇科、龙虱科、蜉蝣科逐渐增多，其它水生昆虫逐渐减少。

使用偏曼特尔分析可以排除环境因子之间自相关的干扰，如图6（见封三）及表3所示，ASL是影响物种丰富度最重要的环境因子（P＜0.01），COD是影响物种丰富度的重要环境因子（P＜0.05）。 ASL、DO、COD是影响功能丰富度的重要环境因子（P＜0.05）。

3 讨论

河流作为一个完整的生态系统，具有一定程度净化污染物的能力［30］。物种丰富度指数的分布情况可以说明东江干流中下游及西枝江下游采样点位于海拔较低处，且受生活污染的排放及土地利用格局变化的影响较大。东江源海拔较高，河岸为山林，人口较少，主要种植脐橙和蜜橘，属农业区，河流受干扰相对较少，水质较好。东江干流中游城镇人口多，受生活污染物的排放及有少量航运、渔业活动的影响较大，河岸遍布竹林；东江下游多为平地，气候条件适宜，更适合人类生存与城镇建设，因此下游城镇化水平比上、中游明显提高，工业化和城市化进程加速了经济发展的同时也引起了一系列环境问题，大量工厂和货运码头，航运频繁［31-32］。东江下游西枝江是东江的第二大一级支流，发源于惠州市惠东县与河源市紫金县交界处，流经惠东县、惠阳区、惠城区后汇入东江［33］。西枝江上游区域有自然保护区白盆珠水库，居民少，污染排放少，水质情况较好；而中游和下游区域人口相对集中，工业较为发达，点源污染占到污染源负荷的85%以上，西枝江下游部分水域处于东江干流剑潭库区内，水生态状况直接影响东江干流水质［34］。石马河为东江的一级支流，源于深圳市宝安区龙华镇大脑壳山，上接观澜河，宝安区工业企业分布密集，工业废水排放较多，剧烈的人类活动使石马河受污染严重［35］。浰江是东江流域一级支流，张辉等［36］研究表明浰江水体优势种较多，优势种的多样化也反映了浰江水体生境的多样性，整体而言，浰江受人类活动干扰小，河网结构相对简单，对生态变化抗性较弱。