将生态系统信息纳入渔业管理科学建议——以海达瓜依太平洋鲱鱼为例的案例研究

本综述系统探讨了如何将生态系统信息整合到渔业管理的科学建议中，以加拿大不列颠哥伦比亚省海达瓜依太平洋鲱鱼（Clupea pallasii）为案例，重点分析了环境、生物及人为压力对其种群动态的影响。研究通过多模型方法（包括线性模型、广义可加模型GAM、梯度森林分析GF及动态因子分析DFA）识别了22个显著的压力-响应关系，揭示了捕食者丰度、上升流强度及太平洋年代际振荡（PDO）等关键压力因子。研究构建了风险分类表，提出需将生态系统考量纳入渔业管理策略，以促进更可持续的资源管理实践。

　　

摘要

本研究旨在推动生态系统信息与渔业管理科学建议的整合，聚焦于加拿大不列颠哥伦比亚省海达瓜依地区的太平洋鲱鱼（Clupea pallasii）案例。尽管近几十年来缺乏商业捕捞，海达瓜依鲱鱼的生产力仍持续低迷，表明非捕捞相关的物理和生物压力仍在产生影响。研究目标包括识别和总结假设通过机制影响鲱鱼的环境、生物及人为压力指标；采用多模型方法分析压力-响应关系；并为渔业管理者提供见解。单压力-响应分析揭示了22个显著的线性或非线性关系，主要涉及捕食者和猎物丰度。梯度森林分析还识别出捕食者丰度、鲱鱼商业捕捞量、上升流和太平洋年代际振荡（PDO）为重要压力因子。动态因子分析揭示了与鲱鱼捕捞量减少和捕食者丰度增加相对应的长期压力趋势。基于生态系统和种群动态的风险分类表显示对该种群的关注度增加。这些发现倡导将生态系统考量纳入渔业管理，以促进更有效和可持续的管理策略。

引言

管理机构日益需要实施生态系统方法渔业管理（EAFM）。例如，加拿大修订的《渔业法》要求渔业和海洋部（DFO）在提供总允许捕捞量（TAC）的科学建议时“考虑影响种群的环境条件”。该法还允许渔业、海洋和加拿大海岸警卫部部长在决定捕捞选项时考虑生态系统方法。因此，DFO已启动国家倡议，以制定和实施EAFM框架。在此倡议下的一个案例研究中，开发了生态系统考量，以纳入海达瓜依太平洋鲱鱼的科学建议。

作为小型中上层饵料鱼类，太平洋鲱鱼（以下简称鲱鱼）种群规模易受大幅波动影响。这种变异性可能源于补充量和自然死亡率在早期生活史阶段和成鱼期的波动。例如，估计的海达瓜依鲱鱼产卵生物量在1990年代中期降至近历史低点，在1990年代后期有所增加，但自2000年以来持续处于低位，多年生产力低下，阻碍了种群增长。自然死亡率（由种群评估模型估计）自1980年以来增加，而生物参数（如体况（通过体长-体重残差衡量）和年龄体重）从1980年代到约2010年下降。海达瓜依鲱鱼的这些趋势可能源于多个物理、生物和人为压力在不同生活史阶段的影响。

在冬末和春季，成鱼（3龄以上）从越冬区迁移到沿海水域产卵，将卵产在浅水大型植物上，包括潮间带和亚潮带水域的鳗草、岩藻、海带和藻类。鲱鱼卵孵化约需2-3周。新孵化的鲱鱼仔鱼依赖卵黄储备1-2周，之后以小型卵、桡足类和硅藻为食。幼鱼（0龄）通常第一个夏季在沿海区域度过，白天近底，夜间分散在上水层。0龄鲱鱼摄食桡足类、藤壶无节幼体和腺介幼体、幼形类、太平洋磷虾成体和幼体、端足类、蟹大眼幼体和蚤状幼体，以及其他无脊椎动物，如枝角类、介形类、昆虫、毛颚类和等足类。一般认为，大多数1龄及以上鲱鱼迁移到大陆架上的夏季摄食区，主要摄食磷虾和桡足类，尽管0龄鲱鱼在这些区域也有观察到。在秋季和冬季，未成熟和成熟鲱鱼使用近岸越冬/滞留区。成熟鲱鱼然后在早春移动到产卵区。到3龄时，90%的鲱鱼性成熟并加入成体产卵种群。

存活假设与影响鱼类早期生活的过程相关。海洋输送、中尺度海洋学特征、保留和集中猎物的海流以及为表层水域提供营养物质的上升流事件与仔鱼存活和补充量相关。已识别对鲱鱼重要的一个物理压力是上升流的时间、持续时间和强度，它直接影响鲱鱼可获得的主要和次级生产。此外，强上升流可能导致仔鱼保留减少，从而降低仔鱼存活率。其他物理压力包括温度、盐度、海平面、河流排放和埃克曼输送。例如，鲱鱼补充量和存活与水温负相关；尽管夏季鲱鱼分布似乎与温度呈穹顶形关系。温度影响卵孵化时间，并与盐度一起影响仔鱼卵黄利用率。盐度也可能影响产卵时间。

影响鲱鱼补充量和存活的生物压力包括浮游动物猎物可获得性、捕食者的存在（特别是太平洋无须鳕Merluccius productus）以及竞争（如与太平洋沙丁鱼Sardinops sagax）。浮游动物物种组成从较大尺寸和脂质丰富的北方和亚北极桡足类物种转变为较小尺寸和脂质贫乏的南方桡足类物种，影响了鲱鱼可获得食物的质量。浮游动物生物总量与鲱鱼生长相关，浮游植物和浮游动物生产与鱼类生产力和产量相关。产卵鲱鱼与随后幼鱼猎物可获得性之间的时间或匹配-失配在决定秋季0龄鲱鱼丰度方面起着至关重要的作用。鲱鱼捕食者和竞争者（如太平洋无须鳕和太平洋沙丁鱼）丰度的增加通常与水温升高相关，可能部分解释了观察到的海表温度（SST）与鲱鱼补充量和存活之间的负相关关系。

物理和生物压力可以在种群评估模型之外进行定性和定量总结。这些物理和生物压力可以直接或间接整合到种群评估模型中。尽管物理压力尚未明确纳入海达瓜依鲱鱼种群评估，但在该种群管理策略评估（MSE）中使用的操作模型开发中测试了时变自然死亡率。时变自然死亡率的基本原因尚未识别且未与物理或生物压力联系起来，但相信反映了生态系统内的潜在变化。目前，生态系统考量以两种方式隐含地纳入鲱鱼科学建议：（1）通过在闭环模拟中模拟鲱鱼自然死亡率的潜在变化；（2）允许在种群评估中时变自然死亡率，直接影响1年预测。然而，当前方法未识别死亡率变化的机制，因此无法预测死亡率的未来趋势。定性和定量总结及分析有助于解决这一知识差距。

本研究的目标是开发相关的生态系统信息，将其纳入渔业管理科学建议，作为海达瓜依鲱鱼评估的一部分。为实现此目标，研究设定了以下目标：（1）识别假设通过机制影响海达瓜依鲱鱼不同生活史阶段的环境和生物压力；（2）识别这些压力的指标；（3）总结压力和潜在海达瓜依鲱鱼响应指标的状态和趋势；（4）检查压力-响应关系；（5）开发一种标准化格式，以有效沟通生态系统信息和风险给管理者。

方法

研究区域

在不列颠哥伦比亚省，鲱鱼被单独管理为两个次要和五个主要种群评估区域，其中之一是海达瓜依。海达瓜依鲱鱼的主要洄游种群通常在4月在海达瓜依东侧产卵，并根据标记研究，夏季移动到赫卡特海峡（HS）摄食。本研究包括了来自HS和夏洛特皇后湾（QCS）区域的生态系统信息，以及影响该区域的大气和海洋学压力的盆地尺度指标。

压力-结果映射

驱动因素（如大气压力、温室气体排放增加、海鲜需求和自然变异性）导致物理、生物和人为压力。通过系统文献映射，识别了15个物理、生物和人为压力，与21个饵料鱼类结果相关，包括分布、生长和生活史、迁移、性能、生产力、繁殖和存活等方面，可能通过机制影响饵料鱼类。

在当前研究中，扩展了文献覆盖范围（1931-2020年）。识别的压力被映射到不同生活史阶段的结果，产生了卵、仔鱼、幼鱼和成鱼鲱鱼的压力-结果矩阵。压力-结果矩阵使用R（v. 4.1.0）和ForageFishLitReview包组装。

生态系统压力指标

压力指标的选择基于压力-结果矩阵、定向文献回顾和专家知识，考虑了鲱鱼按季节和生活史阶段的空间分布。此过程确保满足指标选择的关键标准和考虑因素，包括理论基础和历史数据的可用性。DFO的年度太平洋海洋状况报告是定位时间序列数据和研究人员获取专家知识和额外时间序列数据的主要资源。

压力-结果矩阵中识别的所有物理、生物和人为压力的指标均被包括。物理压力指标包括区域（如SST）和盆地尺度（如PDO和北太平洋环流振荡NPGO）的指标，以及季节到年度基础。所有盆地尺度大气型指标均按年度计算。盆地尺度海洋学或大气指标（如PDO和NPGO）可导致区域海洋条件的变化，如海洋温度、盐度、溶解氧（DO）、海流、营养物质和输送；这些海洋条件变量在研究区域内计算。为考虑区域海洋条件对鲱鱼不同生活史阶段的潜在影响，使用研究区域内可用的月度测量获取季节性SST。还使用了来自区域海洋建模系统（ROMS）BC大陆边缘模型的盐度和DO的季节模型估计作为研究区域内的区域指标。

生物压力指标包括猎物场的变化，如总浮游动物和群落浮游动物生物量和丰度估计（及其变异性），以及浮游植物动态（对鲱鱼有潜在间接影响），如叶绿素a（Chl-a）和春季水华日期。还包括捕食者（如底层鱼类和海洋哺乳动物）和竞争者（底层鱼类和其他浮游动物食者的早期生活史阶段）的指标，如太平洋无须鳕、太平洋刺狗鲨、太平洋大比目鱼、座头鲸和斯特勒海狮的生物量或丰度，以及面向对象海洋生态系统模拟器（OSMOSE）模型估计的捕食者消耗的鲱鱼生物量。最后，人为压力指标包括海达瓜依鲱鱼的商业上岸量和研究区域内商业渔业兼捕鲱鱼量。所有生物压力指标均针对研究区域计算。

鲱鱼响应指标

在压力-结果文献映射中，鱼类结果包括：分布和栖息地使用；生长和生活史（如质量或大小、体况、摄食）；存活和死亡率；性能（如应激、能量、游泳）；迁移（如产卵迁移）；繁殖（如成熟、产卵、卵丰度、卵大小或重量、繁殖力）；和生产力（如生物量或丰度、补充量）。鲱鱼响应指标源自DFO种群评估计划和生态系统建模的可用数据。种群评估指标包括3龄体重（Weight_Age3）、补充量丰度的中位数和可信区间（记为R_Med和R_CI）、产卵生物量的中位数和可信区间（记为SB_Med和SB_CI）以及雄性和雌性体况（即双对数转换体长-体重回归残差的均值和方差，记为Mean_LWResid和SE_LWResid）。生态系统模型指标包括OSMOSE估计的海达瓜依鲱鱼捕捞、自然、捕食和饥饿死亡率，以及平均体型。这些为压力-结果矩阵中的一些主要结果类别提供了指标，特别是：生长和生活史（质量或大小、体况）；存活和死亡率；和生产力（生物量或丰度、补充量）。

状态和趋势

大多数鲱鱼响应指标时间序列跨度从1950年至今，因此，压力指标也为此时间段组装（如可用）。指标时间序列的状态和趋势在标准化为均值为0、标准差为1后进行评估和描述。时间序列状态被分类为阳性、中性或阴性，如果时间序列最后5年的均值大于、 within 或小于时间序列均值的一个标准差。状态可能随时间序列长度变化，因此为最大化信息，使用最长的可用时间序列确定状态。时间序列趋势被分类为增加或减少，如果在序列最后5年观察到显著（p < 0.05）线性趋势（正或负）。如果没有明显趋势，则分类为中性。

分析

应用多模型方法确定压力-响应关系的函数形式，并识别线性和非线性关系。分析包括：（1）线性和非线性压力-响应关系；（2）多变量梯度森林（GF）分析；（3）多变量动态因子分析（DFA）。用于分析的R代码发布在GitHub上。对于所有分析，通过选择单个代表性指标减少了高度相关指标的数量，其中压力或响应指标之间的相关性超过0.80。例如，在压力中，ENSO_ONI与ENSO_MEI相关，因此仅ENSO_MEI被包括在分析中。此外，桡足类丰度和生物量时间序列的均值和标准误相关；因此，包括丰度和生物量的均值，排除标准误。冬季和春季SST相关；因此，仅春季、夏季和秋季SST被包括在分析中。平均南方桡足类生物量和丰度时间序列相关，因此仅生物量被包括在分析中。座头鲸和斯特勒海狮丰度相关；它们在单压力-响应模型中保持分开，但在多变量分析中相加。在响应指标中，均值与鲱鱼补充量丰度和产卵生物量估计的变异性高度相关，因此仅包括这些时间序列的均值。最后，仅考虑雌性体长-体重残差的均值和标准误（即Mean_LWResid_Female和SE_LWResid_Female），排除雄性对应项，以避免冗余。

单压力-响应关系

为探索指标之间的潜在关系，使用四种类型的模型检查单压力-响应关系：线性模型、带自相关的线性模型（LMAC）、广义可加模型（GAM）和广义可加混合模型（GAMM）。如果残差偏离正态分布，则对数据进行对数转换（ln + 1）。模型选择程序遵循概述的程序。GAMM和GAM使用对数似然比检验进行比较；如果发现时间自相关（由值 < 0.05指示），则将GAMM与LMAC进行比较。使用校正小样本大小的Akaike信息准则（AICc）和估计自由度比较GAMM与LMAC或GAM与线性模型。如所述，“当以下情况时，选择非线性模型 over 线性模型：（i）估计自由度 ≥ 2，（ii）AICc值差异 ≥ 2，（iii）广义交叉验证得分较低（GCV）。”在最佳拟合模型中，仅那些p值 < 0.05且R2 > 0.3的模型被考虑在结果中。使用最大5个节点的薄板样条拟合GAM。GAM在R中使用“mgcv”包运行。非线性压力-响应关系中的阈值被识别为压力值，在该值处二阶导数的符号发生变化。为此，通过自举残差（允许自相关）获取阈值可能发生的压力值范围，并计算平滑函数和二阶导数的95% CI。阈值被定义为该范围内二阶导数最不同于零的值。

梯度森林分析

多变量GF分析用于识别沿物理、生物和人为压力梯度的多个鲱鱼响应中的阈值。GF分析整合了不同响应上的个体随机森林，提供了随机森林的功能，同时捕获多个压力-响应关系。随机森林本质上是回归树技术，并构建决策规则以基于单个压力x递归地将数据划分为连续更小的组。森林中的每棵树使用压力的随机子集构建，最小化它们之间相关性的潜在影响。在回归树中，选择分裂值s以最小总不纯度，即关于组均值的平方偏差之和。回归树的第一次分裂s通常接近阈值的值，响应对此特别敏感，并且超出该值响应值转移到另一个水平。沿压力x的分裂s的重要性值从每棵树聚合，并通过核密度估计估计密度曲线。每个响应指标y的拟合优度R2表示随机森林解释的方差比例。运行梯度森林1000次以获得R2的均值和标准差。然后使用具有最高整体性能的运行导出GF分析的所有统计量。计算累积鲱鱼对压力的响应，并且对于那些将R2增加≥0.01的压力，基于压力值范围识别阈值。

分析和结果在两个时间段之间进行比较：1983-2018年，有5个响应和25个压力；和1995-2018年，有5个响应和36个压力。选择较长的时间段以最大化时间序列长度，同时保持大多数压力时间序列。后者时间段最大化压力时间序列数量，同时保持足够的年份进行有意义的分析。分析使用R包gradientForest执行。GF结果与单压力-响应模型结果进行比较。

动态因子分析

为识别趋势并探索多个鲱鱼响应和压力指标之间假设的机制关系，使用了多变量DFA。对于与GF分析相同的两个时间段（1983-2018和1995-2018），使用DFA检查鲱鱼响应（模型1.1和1.2，其中“0.1”和“0.2”代表两个时间段）和压力指标（模型3.1和3.2）的趋势，并将鲱鱼响应指标建模为与某些压力指标相关的多个状态过程的线性组合。DFA识别了最简约的基础状态过程数m及其权重（即Z矩阵表示这些状态过程如何线性组合），使得鲱鱼响应变量可以由大量压力指标最好地解释。对与GF分析相同的时间序列进行分析，并对两个时间段1983-2018（模型1.1、2.1、3.1、4.1和5.1）和1995-2018（模型1.2、3.2、4.2和5.2）的数据进行归一化。如果不正态分布，则转换时间序列。对具有正值的