程佩瑄， 孟凡生， 王业耀， 张铃松， 杨 琦， 薛 浩

1.中国地质大学(北京)水资源与环境学院， 北京 100083 2.中国环境科学研究院， 北京 100012 3.中国环境监测总站， 北京 100012

近年来，随着人类活动的影响，水环境污染物的输入导致水生生物群落的组成和结构发生了系统性的改变，进而影响河流水生生态系统的健康. 因此，为了保护生物群落结构的完整性以及保持生态系统的稳定性，了解生态系统所能承受的胁迫范围十分必要. 生态阈值是生态系统从一种状态快速转变为另一种状态的某个点或一段区间，导致这种转变的原因在于某个或多个关键生态因子较小的或者连续性的改变[1-3]. 推导并确定水环境理化因子的生态阈值，对保护群落结构的完整性及水生态系统的健康具有重要意义.

研究表明，生物参数与多种环境因子的关系为非线性，大量非线性统计方法如非参数突变点分析(nonparametric change-point analysis, nCPA)[4]、贝叶斯突变点分析(Bayesian changepoint)、分段回归分析(piecewise regression)等被应用于推导生态阈值研究中[5]. 然而，这些方法只能使用在聚集的群落水平上而不是物种水平上[6-8]. Baker等[9]提出了临界指示物种分析法(threshold indicator taxa analysis, TITAN)，可通过确定物种的变化点及其对环境因子的响应方向为生态阈值的确定提供依据. 目前针对生态阈值的研究多集中在整个流域，然而在不同地区或者不同生境条件下，水生生物的优势种类存在差异，其耐受性也有所不同，在区域间气候、地形地貌等自然环境因素的协同作用下及调查区域内环境变量取值范围等差异的影响下，各物种在不同区域内对环境因子的响应存在差异，进而导致同一物种在不同区域呈现不同的指示特征[10-12]. 将流域按照地形、水文等特征进行分区，并研究不同分区内环境因子的生态阈值对于流域的管理具有适用性.

底栖动物具有生活在水底、生活史相对较长、移动能力较弱、易于采集等特点，大型底栖动物种类对环境变化的敏感性不同，其群落结构变化能很好地反映环境污染状况[13-14]，因此经常被应用在河流健康的生物评价中. 松花江是我国东北地区的重要河流，近年来松花江流域处于轻度污染状态，主要污染指标为ρ(CODMn)和ρ(NH3-N)，亟需开展松花江流域的生态阈值研究，且流域内高程相差较大，为排除地形对底栖动物群落结构与水质指标之间响应关系的影响，该研究以海拔作为分区依据，采用TITAN法推导不同地形分区内底栖动物群落组成的水质指标的阈值，在此基础上探讨不同地形分区内水质指标对大型底栖动物群落的影响，并识别不同分区水质指标的指示物种，以期为松花江流域不同地形分区生态阈值的确定提供一定的数据支撑，并为不同分区底栖动物群落的保护提供依据.

1 研究方法

1.1 研究区概况与样品采集

松花江流域(41°42′N～51°38′N、119°52′E～129°31′E)位于我国东北地区，是我国七大流域之一，流域面积5.57×105km2. 松花江流域地势南北高、中间和东部低平，绝对高程相差极大，该地势地貌特征决定了松花江流域河流水系的流向和水文特征，进而影响到河流水生态系统的空间分布特征. 流域内多年平均降水量一般在500 mm左右，降水量分布整体呈现山丘区大、平原区小的特征[15-16].

图1 松花江流域采样点分布Fig.1 Distribution of sampling sites in Songhua River Basin

该研究于2016年7月、9月和2017年7月、9月以及2018年7月分别在松花江流域的典型河流断面设置97个采样点，并采集水样及生物样品. 为研究松花江流域不同地形分区的底栖动物群落保护生态阈值，利用ArcGIS 10.2软件中的hydrology工具，对30 m空间分辨率的海拔(DEM)进行子流域提取[17]，共提取257个子流域，同时进行坡度分析；结合坡度和海拔，将各子流域分为平原区(200 m以下)、丘陵区(200～500 m)和山区(500 m以上)(见图1)[18]，以此为基础将采样点划分为山区采样点、丘陵区采样点和平原区采样点.

1.2 样品采集

1.2.1水质指标

水样在河面以下约0.3 m处采集，水温(T)、溶解氧〔ρ(DO)〕、pH和电导率(EC)均使用手持式多参数水质分析仪(YSI Incorporated, Yellow Spring, Ohio, 美国)在野外实地测定. 其他参数测定所需水样存放在0.5或1.5 L的塑料瓶中，保温箱保存，寄回实验室，0～4 ℃下贮存，及时检测分析.ρ(CODMn)、ρ(NH3-N)、ρ(TN)、ρ(TP)等指标的测定均参照《水和废水监测分析方法》(第4版)[19].

1.2.2底栖动物

共采集148个底栖动物样品：对于可涉水河流，样品采集采用D形抄网，每个采样点采集10个30 cm×30 cm的样方，这10个样方涵盖了采样点附近100 m河段的主要生境，最终混合成一个样品；对于不可涉水河流，使用1/16 m2彼得逊采泥器进行采集，每个样点采集3次混合成一个样品. 每个采样点所采集的样品经过60目(0.25 mm)尼龙筛网现场筛洗，剩余物置于白色瓷盘中，将所有大型底栖动物逐一挑出，并用95%的酒精保存后带回实验室鉴定. 大型底栖动物的鉴定在解剖镜和显微镜下进行，所有样品均鉴定至尽可能低的分类单元. 每个采样点所采到的大型底栖动物按不同种类准确地统计数量并用电子天平称量，称量前用滤纸吸去表面固定液. 根据每个样点的采样面积，最终换算出各样点每种大型底栖动物的密度(ind./m2)[20-21].

1.3 分析方法

为研究不同分区水质指标的分布特征，利用单因素方差(one-way ANOVA)分析不同分区之间水质指标的差异. 该研究利用临界指示物种法(TITAN)确定松花江流域底栖动物群落的水质指标的生态保护阈值. TITAN法结合了指示物种分析法和非参数突变点分析法(nonparametric change-point analysis, nCPA)[22]，其原理是对群落中全部物种的水质指标突变点进行比较，当有多个物种在某一较小的浓度变化范围内同时发生相似响应时，该浓度范围即为群落的响应阈值. 在进行数据分析前，排除出现频次小于5的物种. 通过TITAN法找出每个物种的最佳变化点，以得到最大指示值得分(IndVal)，所有物种都会根据在变化点两侧的相对丰度和频率区分为负响应物种和正响应物种[9,23-24]. 同时，对250个观察值样本进行500次自举重抽样，并得到各类群变化点的不确定性(uncertainty)、纯度(purity)和可靠度(reliability)，其中，不确定性指的是突变点分布与自举重抽样所得数据集分布的相异程度，表征从抽样数据集中得到突变点的可能性，纯度和可靠度表征了在抽样过程中有效性和方向性的一致程度；然后以不确定性(p<0.05)、纯度(purity≥0.95)、可靠度(reliability≥0.95)为依据，验证阈值和对应指示物种的可靠性[6].

单因素方差在SPSS 22.0软件中完成，TITAN在R v3.4.2软件中完成，采用mvpart软件包以及Baker和King编写程序[9].

2 结果与分析

2.1 不同分区水质指标和底栖动物群落特征

2.1.1水质指标特征

单因素方差分析表明，各水质指标在山区、丘陵区和平原区之间差异显著(P<0.05)(见表1)，各水质参数的平均值呈现出山区最好、平原区最差的趋势. 松花江流域水质状况整体表现为山区和丘陵区优于平原区，原因可能在于，随着社会经济的快速发展，平原区土地开发利用活动强于丘陵区和山区，导致平原区的水质差于丘陵区和山区.

2.1.2底栖动物群落结构特征

从采集的148个底栖动物样品中共鉴定出大型底栖动物216种，隶属于3门6纲22目61科. 其中，寡毛纲7种，蛭纲8种，腹足纲22种，双壳纲5种，昆虫纲167种，甲壳纲7种. 松花江流域不同地形分区间物种种类数有较大差异，丘陵区物种数(149种)大于平原区(110种)和山区(97种)，主要差异体现在昆虫纲种类的分布上，丘陵区昆虫纲种类数最多，为113种，而平原区仅出现78种，山区共发现86种. 山区出现频率较高的物种主要有亚美蜉属(Ameletus)、弯握蜉属(Drunella)和长跗摇蚊属(Tanytarsus)，出现频率分别为36.84%、26.32%和31.58%；丘陵区出现频率大于20%的物种主要有四节蜉属(Baetis)和壳粗腹摇蚊属(Conchapelopia)，出现频率分别为26.39%和29.17%；平原区出现频率大于20%的物种主要有短沟蜷属(Semisulcospira)和秀丽白虾(Exopalaemonmodestus)，分别为21.2%和26.96%.

图2 松花江流域各分区底栖动物相对密度和相对生物量Fig.2 The relative abundance and biomass of macroinvertebrate in different regions of Songhua River Basin

表1 松花江流域丘陵山区和平原区水质指标特征

松花江流域山区、丘陵区和平原区的大型底栖动物类群组成差异较大(见图2)，其中，昆虫纲相对密度均最高，分别为59.75%、81.47%和62.68%，昆虫纲中摇蚊科的相对密度分别为11.21%、15.77%和34.74%；甲壳纲的相对密度在山区为38.67%，在丘陵区和平原区分别为1.84%和9.42%；寡毛纲的相对密度在平原区为20.35%，远高于山区(1.14%)和丘陵区(5.39%)；腹足纲在丘陵区的相对密度大于山区和平原区，为10.89%；蛭纲相对密度在3个分区内的差异不大. 从相对生物量来看，昆虫纲的相对生物量在丘陵区最高，为67.39%，其次为山区(60.17%)，而在平原区仅为11.59%；双壳纲的相对生物量在平原区最高，为57.66%，而在丘陵区仅为19.73%；腹足纲的相对生物量在平原区远高于山区，分别为19.39%和0.47%；甲壳纲的相对生物量在山区最高(31.42%)，远高于丘陵区(4.02%)和平原区(9.61%)，其他生物类群的相对生物量在3个分区间差异较小.

2.2 不同分区水质指标的阈值

利用TITAN法分别计算得到松花江流域山区、丘陵区和平原区水质指标的生态阈值(见表2、图3)，根据计算结果可知，除ρ(DO)和ρ(CODMn)以外，松花江流域其他水质指标的负响应阈值(即响应变量随环境因子梯度增加而减少)均表现为山区<丘陵区<平原区，而ρ(DO)的负响应阈值表现为山区>丘陵区>平原区，ρ(CODMn)的负响应阈值(5.46 mg/L)在丘陵区最高，在山区最低(4.01 mg/L)；除ρ(DO)以外，其他水质指标的正响应阈值(即响应变量随环境因子梯度增加而增加)均呈现出山区<丘陵区<平原区的趋势，而ρ(DO)的变化趋势则与之相反. 个别水质指标〔如EC、ρ(NH3-N)、ρ(TN)和ρ(TP)〕的正响应阈值在不同分区内变化较明显：在山区、丘陵区和平原区，EC的正响应阈值分别为197.80、314.10和326.00 μS/cm，ρ(NH3-N)的正响应阈值分别为0.23、0.42、0.68 mg/L，ρ(TN)的正响应阈值分别为0.94、2.30、4.64 mg/L，ρ(TP)的正响应阈值分别为46.88 μg/L、213.11 mg/L和257.00 mg/L.

根据水质指标阈值的变化范围，在山区分别有78%、33%、66%、71%、66%和74%的采样点超过了EC、ρ(DO)、ρ(CODMn)、ρ(NH3-N)、ρ(TN)和ρ(TP)的负响应阈值，分别有32%、76%、37%、32%、29%和28%的采样点超过各水质指标的正响应阈值；在丘陵区分别有80%、8%、56%、78%、88%和86%的采样点超过了EC、ρ(DO)、ρ(CODMn)、ρ(NH3-N)、ρ(TN)和ρ(TP)的负响应阈值，分别有10%、74%、22%、37%、14%和8%的采样点超过了正响应阈值；而在平原区分别有92%、9%、76%、88%、76%和93%的采样点超过EC、ρ(DO)、ρ(CODMn)、ρ(NH3-N)、ρ(TN)和ρ(TP)的负响应阈值，其中6%、51%、40%、6%、6%和6%的采样点超过其正响应阈值.

表2 松花江流域山区水质指标生态阈值

图3 松花江流域底栖动物负响应物种(z-)和正响应物种(z+)总指示得分〔Sum(z)〕对水质指标突变点的响应曲线Fig.3 Titan Sum(z-) and Sum(z+) values corresponding to all candidate change points along the water quality parameters gradient for macroinveterbrate in Songhua River Basin

2.3 指示物种

经过TITAN法分析筛选出不同分区各水质指标的指示物种(见表3、图4～6)，水质指标阈值的指示物种在不同分区具有差异. 对山区和丘陵区EC负响应阈值有较好指示作用的物种主要为蜉蝣目，对正响应阈值指示较好的物种主要为摇蚊科和舌蛭属；平原区EC的正响应物种中，其他各指示物种的突变点都在负响应阈值和正响应阈值之间，只有水丝蚓属的突变点(367.9 μS/cm)大于正响应阈值，说明此时只有水丝蚓属的密度发生突变，而其他大部分物种已达到耐受极限，群落结构对EC不再产生明显响应.ρ(DO)在3个分区的正响应指示物种主要为蜉蝣目，四节蜉科一种、颤蚓属和水丝蚓属分别对山区、丘陵区和平原区的ρ(DO)负响应阈值有较强的指示作用.ρ(CODMn)在山区得到的指示物种中，短丝蜉科一种和埃蜉属分别是具有较好指示作用的负响应物种和正响应物种；ρ(CODMn)在丘陵区得到的最佳负响应和正响应指示物种分别为条纹角石蚕幼虫和萝卜螺属；ρ(CODMn)在平原区仅获得摇蚊属这1种正响应物种，戴春蜓则对负响应阈值有最强的指示作用.ρ(NH3-N)在山区和丘陵区的最佳负响应指示物种均为钩虾属，最佳正响应物种分别为长跗摇蚊属和萝卜螺属；ρ(NH3-N)在平原区获得的最佳负响应物种和正响应物种分别为白虾属和细蜉属；河花蜉属和条纹角石蚕幼虫分别是山区ρ(TN)负响应阈值和正响应阈值的最佳指示物种；蜉蝣科一种及舌蛭属分别对丘陵区负响应阈值和正响应阈值具有较强的指示作用，而舌蛭属在丘陵区仅出现了5次，说明区域罕见种对ρ(TN)阈值的确定有一定影响；在平原区，短沟蜷属和水丝蚓属分别对ρ(TN)负响应阈值和正响应阈值具有较强的指示作用.ρ(TP)在山区筛选得到的最佳负响应和正响应指示物种分别是短丝蜉科一种和长跗摇蚊属；在丘陵区动蜉属和摇蚊属分别是对ρ(TP)负响应阈值和正响应阈值指示性最强的物种；在平原区对ρ(TP)负响应阈值和正响应阈值指示性最强的物种分别是短沟蜷属和水丝蚓属.

相同的指示物种在不同的地形分区内可能具有相反的指示方向，萝卜螺属在丘陵区为ρ(NH3-N)的正响应指示物种，在平原区则转变为负响应指示物种，突变点分别为0.45和0.52 mg/L. 短沟蜷属在丘陵区为ρ(TN)和ρ(TP)的正响应物种，突变点分别为1.16 mg/L和36.30 μg/L，而在平原区则转变为负响应物种，突变点分别为0.84 mg/L和38.00 μg/L.

2.4 讨论

2.4.1不同分区底栖动物水质指标的阈值分析

根据TITAN法分析结果可以发现，松花江流域水质指标的负响应阈值中除ρ(DO)和ρ(CODMn)以外，其余均表现为山区<丘陵区<平原区，ρ(DO)的负响应阈值变化趋势则相反，ρ(CODMn)的负响应阈值在丘陵区为5.87 mg/L，明显大于山区(4.62 mg/L)和平原区(5.75 mg/L)；水质指标的正响应阈值中除ρ(DO)以外，其余均呈现出山区<丘陵区<平原区的趋势，ρ(DO)的变化趋势则与之相反. 松花江流域不同地形分区水质指标生态阈值的空间差异特征说明地形因素影响了阈值的推导结果，平原区土地利用强度和人类活动较强，物种对污染物的耐受度增强，导致其保护阈值的升高，山区和丘陵区由于地形因素的限制人类活动强度较弱，物种较为敏感，因此保护阈值较低. 而ρ(CODMn)的负响应阈值在丘陵区高于山区和平原区的原因可能是，松花江流域丘陵区的部分河流，如汤旺河、梧桐河、嫩江右岸小兴安岭地区等，由于流域内土壤腐殖质含量偏高而造成水中高锰酸盐指数本底值偏高[25]. 当河流的某些水质指标值持续较高时，水生生物群落在这种长期的环境压力下逐渐趋于稳定，对环境干扰的耐受程度有所提高，提高了其响应阈值[22,26].

利用TITAN法分析松花江水质指标的阈值，当水质指标超过负响应阈值时，底栖动物的群落结构和密度开始受到一定的干扰，表现为敏感物种密度的降低，而当水质指标超过正响应阈值时，部分耐污种也达到了耐受极限，开始减少，底栖动物群落的结构将发生明显改变[1]. 因此，可将负响应阈值作为触发底栖动物群落发生变化的最低值，正响应阈值为底栖动物群落的耐受极限值[10]. 以此为依据，松花江流域内50%以上的采样点水质指标值都超过了其负响应阈值，超出正响应阈值的采样点比例在6%～40%之间，说明流域虽受到了一定程度的人为干扰，但干扰程度并不严重.

表3 松花江流域TITAN法得到的水质指标阈值指示物种

续表3

注： 圆圈大小表示指示物种响应强度的大小. 实心圆圈表示负响应物种，空心圆圈表示正响应物种.图4 松花江流域山区TITAN法得到的指示物种Fig.4 Significant indicator taxa derived from TITAN in mountainous zone in Songhua River Basin

注： 圆圈大小表示指示物种响应强度的大小. 实心圆圈表示负响应物种，空心圆圈表示正响应物种.图5 松花江流域丘陵区TITAN法得到的指示物种Fig.5 Significant indicator taxa derived from TITAN in hilly zone in Songhua River Basin

注： 圆圈大小表示指示物种响应强度的大小. 实心圆圈表示负响应物种，空心圆圈表示正响应物种.图6 松花江流域平原区TITAN 分析得到的指示物种Fig.6 Significant indicator taxa derived from TITAN in plain zone in Songhua River Basin

目前，大多数研究主要聚焦在推导营养盐的阈值方面，利用TITAN法推导其他水质指标的生态阈值的研究还不多，各研究所得到的结果也不尽相同. Smith等[27]利用HBI、EPT等底栖生物参数通过nCPA计算得到ρ(TN)和ρ(TP)的阈值分别为 0.41～1.2 mg/L和9～70 μg/L，阈值的下限均低于笔者得到的负响应阈值，上限则在负响应阈值和正响应阈值范围之内. Kail等[28]以欧洲河流为研究对象，得到基于大型底栖动物的ρ(NH3-N)阈值为0.12 mg/L，低于笔者得到的ρ(NH3-N)阈值；Sundermann等[6]利用TITAN得到的ρ(NH3-N)负响应阈值为0.368 mg/L，正响应阈值则为0.704 mg/L，ρ(TN)的负响应阈值为1.224 mg/L，正响应阈值为3.51 mg/L，均高于笔者在松花江流域山区和丘陵区得到的阈值，但与平原区的阈值相似. Baker等[9,29]利用TITAN得到Everglades地区ρ(TP)的负响应阈值为14.4 μg/L，正响应阈值为30.3 μg/L，均低于笔者所得ρ(TP)的阈值. 张莉等[30]利用LOWESS拟合曲线分析了浑太河流域3个生态区基于底栖硅藻完整性指数的ρ(NH3-N)和ρ(TP)的阈值，得到ρ(NH3-N)的阈值为0.13、0.30、1.98 mg/L，ρ(TP)的阈值分别为0.04、0.06、0.20 mg/L，所得3个生态区的阈值变化特征与笔者所得结果相似，都体现出山区等人为活动较少的地区阈值较低以及平原区等城市化程度较高、人为干扰较大的地区阈值较高的趋势. 可见，不同研究之间阈值的差异可能来源于研究区环境因子和生物类群分布的差异，此外，地形等环境因素的不同所导致的土地利用类型差异也会影响阈值推导结果.

2.4.2不同分区水体理化指标阈值的指示物种分析

该研究中不同分区水质指标生态阈值的指示物种存在一定差异. 大型底栖动物群落结构在不同地形分区的分布特征是影响水质指标阈值指示物种的主要原因，不同分区的自然地理状况是造成大型底栖动物群落结构分布差异的主要因素，地形的差异导致山区、丘陵区和平原区土地利用类型和人口密度的不同，进而影响了水质状况和生境状况，这反过来会放大区域底栖动物结构的差异[31]. 松花江流域山区溪流水质状况较好，生境适宜大部分底栖动物生存，因此在该研究中，蜉蝣目和摇蚊科在山区采样点出现频率较高，如亚美蜉属(Ameletus)、弯握蜉属(Drunella)和长跗摇蚊属(Tanytarsus)，在山区的出现频率分别为36.84%、26.32%和31.58%，最终被筛选为山区ρ(TP)的指示物种. 丘陵区水质指标的负响应物种以EPT为主，正响应物种则以腹足纲和摇蚊科为主，主要原因在于昆虫纲(EPT和摇蚊科)和腹足纲是丘陵区的主要底栖动物种类，相对密度分别为81.47%和10.89%[32]. 平原区底质多为泥沙，水流相对较缓，水质状况在3个分区中最差，寡毛纲在平原区的相对密度为20.35%，摇蚊科的相对密度在昆虫纲中为34.74%，可见平原区水环境更适宜寡毛纲和摇蚊科这种环境适应能力强的小个体物种生存[31-32]，更容易被TITAN法筛选为水质指标的正响应物种. 相同的指示物种在不同分区可能具有相反的指示方向. 萝卜螺属在丘陵区为ρ(NH3-N)的正响应指示物种，而在平原区则转变为负响应指示物种；短沟蜷属在丘陵区为ρ(TN)和ρ(TP)的正响应物种，而在平原区则转变为负响应物种. 这说明生境条件的差异影响了底栖动物的耐受程度，短沟蜷属和萝卜螺属均为一般耐污种[33]，在丘陵区由于人为干扰有限，其丰度可在一定范围内随着污染物浓度的升高而增加，而在平原区，人为干扰较大，超出了其耐污值，从而导致物种的指示方向发生改变. 松花江底栖动物群落和水质状况具有空间异质性，采取统一的保护阈值不足以起到保护底栖动物群落的目的，因此在地形分区内研究水质指标阈值，可以更有针对性地保护河流的水生态系统.

3 结论

a) 松花江流域丘陵山区和平原区的水质指标特征和底栖动物群落结构具有较为明显的差异性. 流域水质状况以山区为最好，其次为丘陵区，平原区最差；底栖动物群落相对密度和相对生物量具有空间差异，寡毛纲的相对密度和相对生物量均呈现出平原区>丘陵区>山区的趋势，而甲壳纲则刚好相反，昆虫纲的相对密度和相对生物量均呈现丘陵区>山区>平原区的趋势.

b) 松花江流域水质指标在不同地形分区内的阈值不同，负响应阈值的变化趋势除ρ(DO)和ρ(CODMn)外，其余均表现为山区<丘陵区<平原区，ρ(DO)则刚好相反，ρ(CODMn)则在丘陵区出现最高阈值，山区为最低；正响应阈值除ρ(DO)以外，其余均呈现出山区<丘陵区<平原区的趋势，ρ(DO)的变化趋势则与之相反.

c) 将TITAN法所得的负响应阈值作为触发底栖动物群落发生变化的最低值，正响应阈值为底栖动物群落的耐受极限值. 松花江流域内50%以上的样点水质指标值都超过了其负响应阈值，超出正响应阈值的采样点比例在6%～40%之间，说明流域受到一定的干扰，但干扰程度不严重.

d) 不同分区的自然地理状况、栖境状况和水质状况是造成大型底栖动物群落结构分布差异的主要因素，进而影响了水质指标阈值的确定和指示物种的筛选结果. 同一物种在不同地形条件下可能呈现相反的指示特征，如萝卜螺属在丘陵区为ρ(NH3-N)的正响应指示物种，在平原区则转变为负响应指示物种；短沟蜷属在丘陵区为ρ(TN)和ρ(TP)的正响应物种，而在平原区则转变为负响应物种.