生物多样性度量是生态学重要问题之一, β多样性是除α多样性外的一个关键指标, 它可定义为不同时空群落间的物种成分变化。自从Whittaker提出β多样性概念后, β多样性已长期应用于陆地、水生动植物生态学研究和生物多样性保护规划中^[1-2], 因此, 有关β多样性度量方法已不断被关注和拓展^[3-4]。

在众多的β多样性测量方法中, 以平均非相似性的单指数测量法^[5]、多元排序的空间梯度长度法和方差分配法最常见^[6]。单指数测量法很难定量分辨出不同环境梯度和地理距离分别对β多样性的影响；梯度长度法则主要基于线性模型, 也难以反映生态学中最常见的2种非线性关系, 即物种成分非相似性随环境因子和空间距离增加呈曲线关系^[7-8], 以及物种周转速率与环境因子梯度间呈曲线关系^[9]；方差分配法尽管可以定量探索不同预测变量对β多样性的影响, 但只能对少于4个变量(组)进行分析, 因此, 许多影响因子需要合并为变量组^[10]。为解决这些问题, 广义非相似性模型(GDM, Generalised Dissimilarity Modelling)已得到开发并在不同领域得到大量应用^[11], 其实质是矩阵回归, 把各个预测变量通过样条函数转换为统一单位的生物成分非相似性, 于是, 可用不同预测变量、包括经纬度变量来预测区域β多样性变化, 因此, GDM结果可全面代表研究区域的β多样性格局^[12]。

β多样性格局生成的实质来自2个不同生态现象—地点间的物种替换和丰富度差异, 即β多样性可进一步分为周转和嵌套性2个成分^[13]。所谓周转, 指物种被替换(站位间有不同物种), 而嵌套性则必须具备物种丰富度差异的子集特性^[14], 即群落间不仅存在丰富度差异, 而且1个群落必须是另1个群落的子集(站位间无新物种), 或者说, 1个群落是另1个高丰富度群落非随机物种丢失的结果, 这种物种丢失的非随机过程可由任何因素产生^[15]。然而, 广泛应用的Whittaker β多样性无法分辨这2个成分^[16], 而区分这2个成分可为分析群落形成机制提供一个更强大的工具。如2个研究区域有相等的β多样性值, 似乎β多样性的格局是相似的, 但β多样性的形成机制可完全不同, 因1个区域可来自是高周转、低嵌套性, 而另1个区域则可以是低周转、高嵌套性的结果^[17]。此外, 区分这2个成分对了解生物地理学、生态学和规划保护问题至关重要, 因嵌套性地点间存在丰富度子集特性, 因此, 只保护丰富度最高地点即可, 而周转为主的区域则需要保护更多的地点^[15]。

在这2个β多样性成分中, 对嵌套性的认识实际上比β多样性概念还要早^[15]。嵌套性是二进制生态数据矩阵的属性, 可用于探索物种对生境片段化、对特定环境栖息地的适应性等^[18]。其定量方法研究较多, 早期主要采用矩阵温度及其改进方法, 随后又提出并建议使用NODF(基于重叠和下降填充的嵌套性)方法^[19]。然而, Baselga^[15]近年用只含嵌套性、无物种周转的β多样性案例, 指出了NODF方法不能准确反映Whittaker β多样性变化, 并提出了新的嵌套性计算方法。由于β多样性成分分解新方法刚出现不久, 目前知道这2者成分的相对重要性甚少^[2]。

讫今为止, 国内β多样性研究大多针对陆地植物, 对海洋生物定量研究很少。已有综述指出^[20], 国内针对浮游动物几乎只有α多样性研究, 探讨浮游动物β多样性与环境因子的定量关系目前也缺乏分析方法, 因而是一个薄弱环节。为此, 本文利用先前已报道的不同网目采集的群落数据和同步环境监测数据^[21], 进一步挖掘象山港浮游动物β多样性及其成分变化, 并探索环境因子的驱动作用。目标是探明象山港浮游动物：(1) β多样性的驱动因子及其梯度影响；(2) β多样性周转与嵌套性成分的时空格局特征；(3) β多样性成分与环境因子的相关性, 以期为象山港海洋牧场建设和渔业管理提供科学依据。

1 材料与方法 1.1 调查站位与采样

有关调查站位设置细节、室内样品处理、种类鉴定方法等见先前报道的群落比较研究^[21]。采样时间为2010年3月、7月、9月和12月。按象山港空间水文特征, 把狭长海湾划分为6个区(Ⅰ—Ⅵ), 其中, Ⅰ区为杭州湾冲淡水影响区, 并受区域围垦影响；Ⅱ区是与东海水交换的主要区域；Ⅲ区为狭湾入口, 受港口水交换影响, 还受码头运输和大桥建设影响；Ⅳ区为狭湾区, 其潮汐混合强度较Ⅲ区更弱, 此外还受乌沙山电厂影响；Ⅴ区主要受养殖、国华电厂、颜公河影响；Ⅵ区则受养殖、国华电厂、围垦、凫溪河影响。

在各区划分经纬度网格(每格4.7 km × 5.6 km), 各月分别在各区随机取得一个经纬度网格(去除围垦区、养殖区等受限区域), 并在网格内任意取1个站位, 以获得更多含不同经纬度的空间环境信息(图 1)。采样时, 同时用浅水Ⅰ、Ⅱ型浮游生物网(网目分别为505 μm和160 μm)由底到表垂直拖网2次。全年共得24个站位样品。

1.2 环境因子监测

现场用水质仪(Horiba U-5000, 日本)同步测定表层水0.5 m处的温度(SWT)、盐度(SLT)、溶解氧(DO)和pH；透明度(Trs)用圆盘法、水深(Dep)用绳长法、流速(CV)用智能流速仪(LGY-Ⅱ, 南京)测定。水样在室内用电导率仪测定电导率(CDT), 悬浮颗粒物(SPM)按重量法测定^[22], 叶绿素a (Chla, Chlorophyll a)荧光法测定^[23]。

1.3 β多样性分析 1.3.1 GDM分析

GDM建模时, 以站位间所有物种Bray-Curtis非相似性作为响应变量, 以多项式样条函数(polynomial spline function)变换的环境变量差值作为预测变量, 再用指数函数连接响应变量和预测变量, 变换后得^{[11, 24]}：

式中, i, j代表不同站位, d_ij为站位间群落成分非相似性；b为回归常数；k为预测变量(k=1, 2, …, p), x_ki、x_kj分别为i, j站位的k变量监测值。f_k(x_k)代表经样条拟合的转换函数, f_k(x_k)=a×I－spline, 其中, a为常数, I-spline是样条函数, 该样条函数形式为平滑单调增加、以符合随站位间生态距离越远、群落成分差异性越大的生态特点^[7]。为防止模型过度拟合, 采用各个预测变量分布的0、50%和100%分位数作为样条函数分段拟合节点(knot), 并使用gdm工具包缺省划分的3个样条函数^[11]。有关变量的样条函数转换细节参照文献^[25]。

GDM拟合后, 由于预测变量已通过样条函数统一转换为无量纲的生物成分非相似性单位, 故不同预测变量对β多样性的相对贡献大小, 可通过软件输出的偏响应图(Partial response graph)的曲线高低幅度直接进行比较；而各预测变量在图中的形状变化可指示预测变量不同梯度对β多样性影响^{[8, 11]}。GDM运行后, 同样也会自动去除对GDM无贡献的预测变量, 此时, 该变量的偏响应图不会输出^[26]。

1.3.2 β多样性的成分分解

采用betapart工具包, 把总β多样性分解为周转和嵌套性2个成分。当考虑2个配对站位时, 总β多样性计算采用广泛使用的Sørensen非相似性指数(β_sor), 该测量代表了真正的物种空间周转、丰富度引起的嵌套性差异^[5]。当仅有空间周转而无嵌套性差异时, 可用Simpson非相似性指数(β_sim)进行测量^[27], 即：。式中a代表 2个站位中共有物种数, b、c分别为仅在第1、2个站位出现的物种数, min(b, c)代表取b、c中最小值。

不难计算, 当2个站位有相同物种数时, 因丰富度相同, 嵌套性不发生, 即b＝c时, β_sor＝β_sim, 说明β多样性差异完全由空间周转产生；当2个站位有不相同物种数时, 因b≠c, 则β_sor≠β_sim。很明显, 由于β_sor和β_sim在无嵌套性时是相等的, 因此, 它们的差值则代表了嵌套性(β_nes)成分的测量, 即：β_nes=β_sor-β_sim。当计算多站位情况时, 同样可按照该公式推导计算, 但需要考虑多站位间物种的共有信息, 各成分详细推导过程见文献^[15]。

1.4 数据分析 1.4.1 物种丰度数据处理

Ⅰ型、Ⅱ型拖网样品的2重复浮游动物丰度分别取平均。对Ⅰ型和Ⅱ型网中的相同物种, 取高丰度数据作为该物种有效采样数据。如象山港主要优势种针刺拟哲水蚤(Paracalanus aculeatus)在Ⅰ型网中的丰度仅为0.34个/m³, 而在Ⅱ型网高达1169个/m^3[21], 显然, 使用Ⅰ型网或两网平均值数据, 都不能代表该种在象山港中的丰度, 而用Ⅱ型网的高丰度数据更具代表性。

1.4.2 统计分析

GDM建模时, 使用浮游动物的丰度数据。经纬度采用UTM投影(WGS 1984 UTM Zone 51N), 用以探讨地理距离上的β多样性差异, 在GDM拟合时, 配套的经纬度作为1个预测变量^[11]。

为探讨β多样性的周转与嵌套性成分变化, 物种丰度数据需经P/A(Presence/Absence)转换。由于象山港以桡足类和中小型浮游动物为主^[21], 为探索不同类群可能存在的β多样性成分上差异, 还进一步把浮游动物分为桡足类、幼体类, 以及大型浮游动物类(毛颚类、水母类、端足类、涟虫类、樱虾类、磷虾类、糠虾类、多毛类、仔鱼类)3个类群分别进行探讨。β多样性成分与环境因子间的Pearson相关系数显著性经Mantel permutation排列法检验, 排列999次。以上所有统计分析经R语言(版本2.15.3) 编程, 主要加载包为vegan、gdm^[11]和betapart^[28]。

2 结果 2.1 环境因子监测结果

不同环境因子的时空监测结果列于表 1。由表 1可知, 不同月份盐度(SLT)监测值变化不大。流速(CV)差异与采样时的潮汐有关, 监测点的水深(Dep)除与潮汐变化有关外, 还与采样点地形有关。悬浮颗粒物(SPM)9月明显低于其他各月, 故透明度(Trs)也以9月最大, 且明显高于其他月份。叶绿素a(Chla)以9月最高, 12月最低；电导率(CDT)以12月最低, 而pH低值主要发生在7月和9月。象山港不同季节水温(SWT)变化明显, 以9月最高, 3月最低；溶解氧(DO)变化趋势与水温相反。从空间上分析, 港口到港底, 盐度、流速、pH基本呈下降态势, 而叶绿素a、透明度、水温增加。水深在港口区较浅, 悬浮颗粒物主要分布在港口区与狭湾区, 电导率与溶解氧不同区域有所波动。

2.2 象山港不同浮游动物类群时空差异

由表 2可知, 桡足类丰富度(种类数)4个月份总体相似, 大型浮游动物丰富度以7月和9月较高, 幼体类也以7月和9月较多。空间上, 桡足类丰富度以港口区最高, 港底区最低；大型浮游动物在港口区和狭湾区相似, 都高于港底区；而幼体类空间分布较相似。

2.3 浮游动物β多样性的GDM拟合及预测变量贡献

据GDM输出结果, Null模型的总偏差量为29.40, 而GDM预测变量模型产生的偏差量为7.29, 故预测变量可解释浮游动物β多样性比例为75.2% (即100%-7.29/29.4×100%), 未能解释比例占24.8%。

在构建GDM模型的11个变量中, 未检出流速对β多样性的影响(表 3)。从各变量相对贡献百分比上看, 列于前3位变量分别是表层水温、溶解氧和水深, 累计贡献63.9%。其他变量累计占36.1%, 其中, 经纬度产生的地理距离对β多样性有贡献, 比例占6.3%, 几乎与叶绿素a和pH的相对贡献比例相当, 也高于盐度的相对贡献, 指示在象山港β多样性研究时, 需要考虑采样点间空间距离产生的影响。

2.4 浮游动物β多样性的预测变量梯度影响

预测变量对β多样性的梯度影响见图 2。各变量的梯度影响大体可分为4种类型, 第1种是变量仅在高梯度区有影响, 如来自经纬度的地理距离(Geo)、pH和盐度, 分别在约高于25 km、7.8和25时会促进β多样性增加；第2种是仅在低梯度区有影响, 如水温和叶绿素a分别在约低于22 ℃和0.5 μg/L时, β多样性才随梯度值的增加而增大；第3种是几乎随着变量值的增加一直增加, 如溶解氧、电导率、透明度和水深；第4种是影响极低, 如悬浮颗粒物含量几乎不会对β多样性产生影响。

2.5 象山港浮游动物β多样性成分的时空特征

象山港浮游动物β多样性及其周转、嵌套性成分时空差异列于表 4。对总浮游动物, 总β多样性在7月份最低, 主要原因是其周转最低, 指示该月份物种替换较少, 即7月份群落成分比其他月份更稳定；对比各月份的周转与嵌套性结果, 周转均比嵌套性值高出2倍以上, 显示各月份不同地点间浮游动物总体上都以物种替换为主。空间上, 从港口区到港底区总β多样性呈下降趋势, 但周转却增加, 显然, 嵌套性值从港口到港底区下降是导致空间上总β多样性下降的主因。

据表 4还可知, 3个类群浮游动物的β多样性及成分时空变化也明显不同。首先, 对桡足类, 不同时空上的总β多样性主要由周转引起, 嵌套性较小。在7月份, β多样性较低, 而在其他月份更相似, 主要原因是桡足类周转较低；从港口区到港底区β多样性只略有下降, 说明桡足类在象山港周转和嵌套性都较稳定。其次, 大型浮游动物在7月份总β多样性也最低, 其他月份间相近, 但7月份β多样性中的周转与嵌套性值比较接近, 反映站位间的物种嵌套性主要发生在7月；空间上, 港口区到港底区大型浮游动物总β多样性增加, 这主要由物种周转增加、以及嵌套性值仅在港口区和狭湾区较高引起。最后, 对幼体类, 除7月份总β多样性最低、周转与嵌套性相似外, 其他月份的嵌套性都高于周转, 指示其他月份幼体类β多样性中的嵌套性比周转起更重要作用；空间上, 港口区幼体类的β多样性比港底区高, 主要原因是周转和嵌套性在港口区都高于港底区, 而且嵌套性都高于周转。

2.6 不同类群浮游动物β多样性成分与环境因子相关性

不同类群浮游动物的β多样性成分与环境因子间的相关性、及其Mantel检验见表 5。对桡足类, 经纬度、水深、水温、溶解氧与周转存在显著相关性, 但未检出各因子对嵌套性有显著相关。对大型浮游动物, 经纬度、盐度与物种周转有显著相关性, 水温和溶解氧则与嵌套性显著相关。在幼体类, 仅pH值与周转显著相关, 而流速、水温、电导率、溶解氧与嵌套性显著相关。对总浮游动物而言, 影响浮游动物周转的显著相关因子分别是经纬度、水深、pH、水温、电导率和溶解氧, 而影响嵌套性的显著相关因子是流速、水温和溶解氧。

3 讨论 3.1 驱动象山港浮游动物β多样性变化的主要环境因子

据本研究结果, GDM根据样条函数的分段曲线拟合功能, 不仅可分辨出不同环境因子、地理距离分别产生的β多样性影响, 而且还能识别各变量梯度的影响速率。经GDM拟合, 10个预测变量解释了浮游动物β多样性的75.2%, 该解释比例高于大多数水生生物、陆地生物β多样性研究结果^{[24, 29]}。这说明监测的预测变量可较好地反映象山港浮游动物β多样性变化, 也进一步证实：由于大多数海洋生物调查的物种数据偏少, 往往很难建立代表性的物种模型, 但GDM使用生物成分非相似性, 故对数据点很少时特别有用, 并可广泛用于物种分布分析^[30]。

在监测的变量中, 水温、溶解氧、水深是驱动β多样性变化3个最主要因素。除水深与浮游动物空间容纳量有关外, 由于溶解氧含量往往伴随季节水温变化呈极强的负相关^[31], 说明溶解氧对β多样性的驱动强烈受季节水温影响, 且当水温约在22℃以下时, β多样性随水温上升而快速增加, 该水温上升条件明显发生于3月到7月份。据已有研究, 象山港水温从2月到8月为上升周期, 随水温增加浮游动物种类会快速增加^[32]。从整个象山港浮游动物优势种更替看, 3月到7月种类更替率高达77.8%, 而高温期的7月到9月更替率极低(16.7%)^[21]。这些结果指示, 随着水温上升, β多样性随着浮游动物更替率增加而增加, 到高温期, 伴随溶解氧的降低, β多样性几乎不再增加。

除上述3个主要因素外, 盐度往往是驱动浮游动物空间分布的另一个重要因素^[20], 然而, 本研究发现盐度对β多样性影响较低, 主要原因是象山港盐度梯度周年变化不大；再者, 由于象山港水交换率低, 外海高盐水所携带的浮游动物影响区域有限, 因高盐水充其量也只能到达象山港约1/3处^[33]。本研究还发现, 往往作为站位描述的经纬度(地理距离)却对浮游动物β多样性有影响, 其影响比例高于盐度。当站位间距离大约超过25 km时, β多样性影响开始增大, 事实上, 这反映了生物成分非相似性随地理距离增加而增大特点^{[1, 7]}, 这种非相似性差异显然不是人类活动产生, 因此, 在探讨人类活动对β多样性影响的大尺度研究时, 地理距离尺度本身存在的β多样性差异必须关注, 否则所识别的驱动因子也会不同。如对小型底栖动物线虫的4个尺度大小研究时发现, 在宏观尺度上(>1500 km), 物种周转由食物资源的数量驱动, 而在小尺度上(m), 则由食物资源的质量控制^[4]。此外, 由地理距离尺度产生的β多样性差异甚至还会放大人类活动的影响因子解释, 如在澳大利亚沿岸大型海藻β多样性研究中, 当添加经纬度变量时, 对GDM模型解释量增加了不到1%^[8]；而在墨西哥, 对节肢动物分布的贡献则以地理距离最重要^[9], 显然, 此时若不考虑地理距离影响, 则该影响会计算到人类活动的影响中, 从而放大人类活动的影响。

3.2 不同类群浮游动物周转、嵌套性与环境因子关系

本研究结果显示, 象山港浮游动物β多样性以物种周转为主。对大型浮游动物, 经纬度差异和盐度是驱动其周转变化的显著影响因子。在象山港, 港口与港底区水环境很不同, 早被认为是2个不同的水体, 各自以不同的浮游动物优势种在港中部交汇混合^[32], 因此, 经纬度引起的地理距离相距越远物种差异越大、周转越高。而盐度越高周转加快, 这明显与外海高盐水携带大型浮游动物、从而增加物种丰富度有关。如象山港肥胖箭虫(Sagitta enflata)的出现可指示其来源于暖水性外海种^[33], 主要大型浮游动物-百陶箭虫(Sagitta bedoti)也主要分布于盐度梯度较大、外海高盐水与沿岸水的交汇区, 双生水母(Diphyes chamissonis)也同样与表层盐度有较强相关性^[34]。对幼体类, 物种周转则显著相关于pH值。从监测结果看, 象山港pH值以狭湾区和港底区最低, 且发生在7月和9月份, 主要原因是6—9月份为当地丰水期, 由于上游化工厂多, 径流会带入酸化水, 而酸化水会影响浮游动物幼体发育, 甚至死亡^[35], 从而降低浮游动物的丰富度。而对桡足类, 经纬度、水深、水温和溶解氧均与周转有显著相关性, 由于象山港主要类群是桡足类, 丰度几乎占79%以上(不包括桡足类幼体)^[21], 因此, 这些因素的解释类似于总β多样性, 这已在3.1节中论述。

从象山港浮游动物的嵌套性结果看, 除桡足类嵌套性极低外, 大型浮游动物、幼体类都存在不同程度的嵌套性。从概念上看, 嵌套性指不同站位间的浮游动物种类发生丢失, 而不是指种类发生更替, 因此, 这明显与种类的扩散分布状态有关。对大型浮游动物而言, 嵌套性主要发生在港口区和狭湾区, 驱动因子主要是水温和溶解氧。原因是港口区易受外海影响, 水交换率高溶解氧较大, 通过水交换也会带来东海区大型浮游动物, 但水交换影响通常只侵入到象山港约1/3处^[33], 因此, 从港口到港中部, 随着水动力的减弱, 主要大型浮游动物的毛颚类、水母类、糠虾类种类数分布几乎呈递减趋势^[21]。而这些优势类群在象山港海域主要发生在7—9月的高温期, 此时, 随着台湾暖流靠近长江口, 带来了毛颚类^[32]、水母类^[34]、糠虾类、磷虾类、十足类, 以及仔鱼等象山港优势类群^[21]。在幼体类中, 嵌套性在时空上大多高于周转, 空间上也以港口区嵌套性更高, 驱动幼体类嵌套性的主要因子除与大型浮游动物相似的水温、溶解氧影响外, 还有流速和电导率。众所周知, 幼体类的运动能力比大型浮游动物弱的多, 流速更易改变幼体类的扩散分布, 这是驱动幼体类嵌套性显著变化的一个影响因素, 其中, 港口区流速受潮汐影响更大, 而港底区水动力则很弱, 使得港口区幼体类分布受影响更大。而电导率驱动浮游动物嵌套性变化目前不清楚, 据报道, 电导率与水体中可溶性营养物有关^[36], 因此, 有关电导率与幼体嵌套性关系是否与水体中营养物分布有关有待深入研究。