基于混合效应模型的海州湾方氏云鳚的生长动态及其影响因素分析

赵天亚 ，杨晓改，李勋，薛莹 ，徐宾铎 ，纪毓鹏 ，任一平 ，张崇良 *

( 1. 中国海洋大学 水产学院，山东 青岛 266003；2. 青岛海洋科学与技术试点国家实验室 海洋渔业科学与食物产出过程功能实验室，山东 青岛 266237；3. 海州湾渔业生态系统教育部野外科学观测研究站，山东 青岛 266003；4. 中国农业出版社，北京 100125；5. 青岛西海岸新区海洋发展局，山东 青岛 266400)

1 引言

鱼类的生长特征是渔业资源评估的重要内容，可以反映鱼类群体的生存环境和资源状况。根据可用的数据材料可知，鱼类生长的研究方法主要有两种：基于长度的评估方法和基于年龄鉴定的评估方法。基于长度的方法如电子体长频率分析法（Electronic Length Frequency Analysis, ELEFAN）和MULTIFAN等适用于年龄数据难以获得、长度数据易于获得的数据有限的场景，但易受到采样随机性和算法稳健性的影响[1]。而鱼体硬组织材料如耳石和鳞片等可以通过年轮的形式记录鱼类的生长信息，其轮纹增量可以直接反映鱼类在该年份的生长速度，从而直观地反映鱼类的生长动态，为探讨影响鱼类生长的因素提供了重要材料。

本研究以方氏云鳚（Pholis fangi）为例研究鱼类的生长特征。方氏云鳚作为一种小型饵料鱼类，在食物网中维持关键的能流传递过程，在生态系统中发挥重要作用[2]。伴随着近海渔业资源的明显衰退，方氏云鳚等小型底栖鱼类逐渐成为海州湾及其邻近海域的主要优势种[2]。近年来，受到全球气候变化及人类活动的影响，许多鱼种在生理、行为及分布等方面上均出现多种适应性转变[3]。尽管已有研究估算了方氏云鳚的生长参数[4]，但在气候变化背景下，其生长特征可能发生变化，生长的动态及驱动因素目前尚缺乏探讨。

线性混合效应模型（Linear Mixed-Effects Model,LMEM）是一种应用广泛的统计学模型[5]，通过在一般线性模型中加入随机效应，可以有效解释数据中的嵌套关系对整体线性回归所产生的影响。该模型在渔业领域亦有广泛应用，如Ma等[6]通过线性混合效应模型研究了黄、渤海海域小黄鱼（Larimichthys polyactis）体长体重关系的时空异质性；Mu等[7]使用线性混合效应模型研究了黄海海域星康吉鳗（Conger myriaster）的生长动态。

为掌握方氏云鳚的生长动态及其对环境因子的响应，本研究基于2015-2019年秋季海州湾方氏云鳚耳石样本，利用线性混合效应模型研究方氏云鳚生长速率的年际变化，并分析了影响其生长动态的主要因素。本研究旨在加深对方氏云鳚生物学特性的认识，为资源评估和渔业管理提供科学支撑，为理解气候变化驱动下鱼类的生物学特征响应提供参考。

2 材料与方法

2.1 数据来源

本研究所用耳石样本全部来自于2015-2019年秋季于海州湾及其邻近海域进行的渔业资源底拖网调查，调查海域范围为34°20′～35°40′N，119°20′～121°10′E，站位设计采用分层随机采样的方法（图1），每个航次设置18个站位，在各区域中随机选取站位[8]。所用调查船为功率220 kW的单拖底拖网渔船，拖网速度为2～3 kn，每站拖网时间为1 h。调查网具网目为17 mm，网口扩张高度约为6 m，宽度约为25 m。将所采集样本带回实验室进行解剖，从每条方氏云鳚样本的左右耳石中随机选取一枚进行分析。样本的采集和分析参考《海洋调查规范》（GB/T 12763.6-2007）[9]。

图1 海州湾及邻近海域调查区域Fig. 1 Survey areas in Haizhou Bay and its adjacent waters

本研究共选取107枚方氏云鳚耳石样本，其中2015年共选取22枚样本，2016年共选取24枚，2017年共选取21枚，2018年共选取19枚，2019年共选取21枚。使用Olympus bx51型光学显微镜采集耳石样本图像并读取年龄，结合样本年龄和采集年份反推其每一轮纹的生成年份，再沿耳石长轴从耳石核心起始对每个年轮进行测距。读耳石过程遵守盲读和重复性原则[7]，年轮测距过程通过imageJ软件实现。本研究所选样本包含1～5龄个体，但5龄个体较少，第5年轮的耳石增量数据不具有代表性，因此在数据整理中将其剔除。经整理筛选后共得到海州湾方氏云鳚的241个耳石增量数据，其轮纹形成时间为2013-2018年。

2.2 因子筛选

本研究以耳石增量作为响应变量，以年龄、出生年份、底层海水温度、叶绿素含量和种群密度等为解释变量，构建回归模型。方氏云鳚为一种底栖小型饵料鱼类，营养级较低，因此本研究以海水叶绿素含量反映其食饵条件；同时以轮纹形成时间所对应年份的底拖网调查结果中的站位平均标准化数量反映其密度制约强度，其计算方式为：站位平均标准化数量=∑站位标准化数量/总站位数，单位为ind./站。本研究中叶绿素含量数据及底层温度数据提取自FVCOM（Finite-Volume Community Ocean Model）[10]。以方差膨胀因子（Variance Inflation Factor，VIF）作为因子多重共线性衡量指标对初始因子进行筛选，通常认为VIF＞10即存在多重共线性问题，应在建模前予以剔除。

2.3 线性混合效应模型

本研究使用线性混合效应模型研究海州湾方氏云鳚耳石增量的年际差异。线性混合效应模型包括固定效应和随机效应两部分，固定效应用于解释线性回归的整体期望，随机效应用于解释数据中嵌套关系对整体期望的影响。线性混合效应模型的一般形式为

式中，Yi为响应变量；α为模型截距；Xi×β为固定效应；Zi×bi为随机效应；εi为误差项。本研究根据加入模型的影响因素来源的不同，分别构建了内因模型和外因模型。

2.3.1 内因模型

内因模型考虑了影响生长的种群内在因素。其中年龄通常为生长速度的主要决定因素，因此将年龄作为固定效应加入模型，此外将个体编号、年份和世代作为随机效应加入模型，模型响应变量为耳石增量。使用小样本校正的赤池信息准则（AICc）作为模型选择标准[11]，AICc最小的模型为最优模型。使用极大似然估计法重新拟合最优模型，以得到耳石增量在年际间差异的无偏估计值[5]。

2.3.2 外因模型

外因模型考虑了外在环境因子对方氏云鳚的生长速度的影响，模型响应变量为耳石增量，解释变量包括底层温度、叶绿素含量和种群密度。模型中加入因子二次项，以反映环境因子与生长速度的非线性关系。由于方氏云鳚在幼体时期（小于1龄）生长速度显著大于成体时期（大于1龄），其生活史特征差异较大[12]，因此本研究以0龄时期所对应耳石增量（第1年轮）反映方氏云鳚幼体时期的生长特征，以1龄时期所对应耳石增量（第2年轮）反映方氏云鳚成体时期的生长特征，分别对方氏云鳚0龄时期及1龄时期的耳石增量对环境因素的响应进行建模分析。模型拟合过程采用逐步回归法，根据AICc最小准则分别选取最优模型。

本研究中数据分析、模型拟合及绘图过程均于R语言4.0.2中完成，其中LMEM的拟合由lme4软件包实现。

3 结果

3.1 耳石增量分析

方氏云鳚个体在不同年龄的耳石增量具有明显差异（图2），第1年轮的平均耳石增量为0.327 mm，显著高于高龄耳石增量。第1年轮的耳石增量在2013-2018年呈现先上升后下降的趋势，在2016年达到最大值（0.341±0.024）mm；第2至第4年轮的耳石增量均呈现不规则波动趋势。此外，方氏云鳚在不同年龄阶段的耳石增量波动范围变化较大，变异系数随年龄逐渐增大，表明个体间生长变异水平随年龄逐渐提高（表1）。

图2 海州湾方氏云鳚耳石增量随年龄的变化Fig. 2 Change in otolith increment of in Pholis fangi Haizhou Bay with different ages

表1 方氏云鳚不同年龄耳石增量的变异性Table 1 Variation of otolith increment of Pholis fangi at different ages

3.2 内因模型分析

通过调整随机效应的不同结构，共拟合了3个内因模型（表2）。其中模型M3的AICc值最低，为最优内因模型，包含个体编号和年份两个随机效应，共解释了85.9%的偏差。年龄固定效应表明，方氏云鳚的生长速度随着年龄的增加而显著减缓（图3a）。年份随机效应表明，方氏云鳚耳石增量在2016年达到最大值（0.260±0.004）mm，其生长速度在2013-2016年呈现逐步变快的趋势，在2016-2018年呈现波动趋势（图3b）。

图3 内因模型的年龄与年份效应Fig. 3 Age and Year effects of the intrinsic model

表2 内因模型的拟合Table 2 Fitting of the best intrinsic model

3.3 外因模型分析

多重共线性检验结果表明，底层温度、叶绿素含量和种群密度的VIF均小于10（表3）。因此3个环境因子均可作为初始解释变量。通过逐步回归法分别拟合最优0龄模型（Model 1）和1龄模型（Model 6），两个最优模型均只包含全年平均底温（BT）（表4），分别解释了22.9%和2.5%的偏差。本研究将全年平均底温分离为春季平均底温、夏季平均底温、秋季平均底温和冬季平均底温，并将4种底温分别作为初始解释变量加入模型，以评估不同季节的温度对其生长影响的强弱并改良模型的拟合效果。

分别将四季底温作为初始变量时，3个环境因子的VIF均小于10（表3）。将季节底温分别作为模型初始解释变量并重新拟合模型，最优模型筛选结果见表4。在0龄模型中，季节底温模型（Model 2～5）分别解释了24.7%～30.1%的偏差，其AICc值均小于Model 1；在1龄模型中，季节底温模型（Model 7～10）分别解释了0.9%～2.9%的偏差，其AICc值均大于Model 6。

表3 因子多重共线性检验结果Table 3 Results of multicollinearity test of factors

表4 最优外因模型筛选结果Table 4 Results of best extrinsic models screening

3.4 环境因素对耳石增量的影响

最优模型拟合结果表明，环境因素的加入可以优化0龄模型的拟合及解释效果（表4）。在0龄模型（Model 1～5）中，所加入的环境因子对0龄耳石增量均具有显著影响效应（p＜0.05）；而在1龄模型中（Model 6～10）中，环境因子对1龄耳石增量均无显著影响（p＞0.05）。整体而言，方氏云鳚0龄时期的生长速度受到环境因子的显著影响，而1龄时期的生长速度则主要受到自身年龄的制约。

在全年平均底温模型Model 1和Model 6中，底温每升高1℃，会导致0龄耳石增量增加14.49%，1龄耳石增量减少13.83%（图4），因此方氏云鳚在0龄及1龄时期的生长速度对底温的响应存在一定差异。

图4 全年平均底层温度对方氏云鳚耳石增量的影响Fig. 4 Effect of annual average bottom temperature on otolith increment of Pholis fangi

在季节底温模型Model 2～5中，四季平均底温对0龄耳石增量的影响效应存在一定差异，0龄耳石增量随着春季底温和冬季底温的升高而增加，随着夏季底温和秋季底温的升高而先增加后减少（图5）。种群密度对0龄耳石增量存在一定的负效应，耳石增量随着种群密度的上升而减少（图5）。叶绿素含量对0龄耳石增量的影响效应在模型间存在差异，在模型Model 2和Model 5中，0龄耳石增量随叶绿素含量的上升而减少；在模型Model 4中，0龄耳石增量随叶绿素含量的上升而增加（图5）。

图5 环境因子对方氏云鳚0龄耳石增量的影响Fig. 5 Effect of environmental factors on otolith increment of Pholis fangi at age 0

4 讨论

本研究基于耳石增量数据，利用混合效应模型研究了海州湾方氏云鳚在2013-2018年生长速度的年际差异及其生长对外界环境因子的响应。结果表明，海州湾方氏云鳚存在一定的生长变异，同时该群体在2013-2018年内的生长速度整体呈现变快的趋势；方氏云鳚在0龄时期的生长速度受环境因子的影响较大，受到底层温度、食饵条件和密度制约的调控；而年满1龄后其生长速度受环境因子变化的影响较小，主要受到自身年龄的限制。

在2013-2018年内，方氏云鳚的生长速度具有明显的变化，2016-2018年的生长速度相较于2013-2015年更快，这与王琨等[4]使用ELEFAN方法所评估的结果较为一致。方氏云鳚生长变快的趋势可能与目前的人类活动及气候变化的大环境具有密切联系。在适温范围内，温度的适当提升可以促进鱼类的新陈代谢和摄食，进而加快鱼类的生长速度，而Gao等[13]的研究表明，2016-2018年黄海南部出现较大程度的海洋变暖，这与本研究中2016-2018年方氏云鳚生长速度较快的现象较为一致。此外，海州湾海域长期的过度捕捞也会使鱼类产生适应性转变[14]，外界压力可能使方氏云鳚选择加快生长速度的策略以快速到达性成熟。值得注意的是，方氏云鳚生长速度的年际差异主要来源于0龄时期生长速度的差异，而年满1龄后的生长速度波动不大，这可能与其生活史特征有关。方氏云鳚为一种小型短周期鱼类，在0龄时期生长速度极快，而性成熟后的生长速度大幅降低，资料显示方氏云鳚年满1龄时的体长即可达到100 mm[15]。通常而言，鱼类在幼体时期的生长和生理状态对外界环境较为敏感，因此其生长速度容易因外界环境因素的变化而产生较大程度的变化；而成体生长缓慢，其生长速度较为稳定，对环境的适应能力较强，其生长过程不易受到外界环境因素的影响。

温度对方氏云鳚0龄时期的生长速度具有显著的正效应，对1龄时期的生长速度具有一定程度的负效应。方氏云鳚的新生群体在孵化后进入为期1年的索饵期，此时会广泛分布于黄、渤海沿岸海域[15]，而黄海沿岸海域在夏、秋季节的底层水温为18～22℃[16]。本研究中Model 1～5的结果表明，温度的升高可以促进方氏云鳚0龄时期的生长，直至21.5℃左右时生长速度达到最快，符合其生活史特征，因此较高的温度可以促进方氏云鳚0龄群体的生长，且该温度可能为0龄幼体的最适生长温度，在其生长过程中起到关键作用。底层温度对方氏云鳚1龄时期的生长速度具有一定程度的负效应，这可能与其冷温习性有关。年满1龄后的性成熟个体通常在产卵后游向深水区，并在夏、秋季节集中于青岛冷水团区域[17]，在冬季重新返回产卵场进行产卵繁殖，其生活史特征与幼体阶段具有较大差异。相较于0龄时期，方氏云鳚年满1龄后的生长速度较为稳定，环境因素的变化对其耳石增量的影响较小，其生长速度的变化主要来自于自身年龄限制。

密度制约也是影响方氏云鳚0龄时期生长速度的一个重要因素。密度制约是一个综合性效应，通常会与鱼类栖息地的理化环境及捕食与竞争等生物过程密切相关[18]，普遍被认为是鱼类补充群体的一个重要调节手段[19]。方氏云鳚为小型短周期鱼类，其生活史策略通常倾向于r选择，倾向于将更多的能量用于繁殖过程，因此冬季孵化后的新生群体将存在较大的密度制约效应[20]，主要表现为对食饵等资源的竞争，其结果为随种群密度的上升，方氏云鳚幼鱼将具有逐渐降低的生长速度及逐渐升高的幼体死亡率。此外，当0龄幼鱼经历为期1年的索饵和发育后，将于冬季与高龄鱼群共同开始向产卵场集中，进行产卵繁殖过程。在此过程中，个体的生长过程因能量分配和个体间竞争而受到限制，同时方氏云鳚的集群习性会进一步扩大密度制约效应对其生长速度的影响。因此，方氏云鳚在0龄时期会经历较大程度的资源竞争，其生长速度受到密度制约较为显著的调控。此外，许多鱼类都有不同程度的生长变异现象[21]，生长变异的程度也可以反映个体间对资源的竞争强度，如蓝纹鲈（Gramma loreto）和 前 唇 霓 虹 虾 虎 鱼（Elacatinus prochilos）等鱼类群体会因资源竞争而在个体间产生较为明显的体长差异[22]。方氏云鳚群体内部较大的竞争强度导致个体间资源利用程度的均匀性较低，从而使该群体产生较大程度的生长变异现象。

叶绿素含量对方氏云鳚0龄时期生长速度的影响效应在模型间存在一定差异，其原因为温度与叶绿素含量间可能存在一定的关系，在模型中加入不同季节的温度作为解释变量时，会导致叶绿素含量对耳石增量的解释效果出现差异。通常认为食饵条件与鱼类的生长速度间存在正相关关系，而叶绿素含量可以反映浮游生物的丰度变化，进而在一定程度上反映方氏云鳚的食饵条件，但浮游生物丰度可能受其他诸多外界因素如温度或其他低营养级鱼类捕食的影响[23-24]，从而导致使用叶绿素含量反映方氏云鳚食饵条件可能存在一定的局限性。此外，有研究指出营养级联效应可能会影响鱼类的捕食者群体动态或饵料丰度，并结合密度制约效应间接调控鱼类的生长[25-26]。其具体影响机制还有待进一步探讨。

年龄与生长是渔业资源评估与管理的重要研究内容，相关研究不仅可以掌握鱼类的生物学特性及其资源现状，还可以为资源评估模型如单位补充量渔获量（Yield Per Recruitment, YPR）模型[27-28]、生态系统模型如OSMOSE等提供基本的生物学参数[29]。耳石作为一种可以记录鱼类生长信息的材料，在国内外已被广泛应用于对年龄与生长的研究[30-32]。耳石和鱼类体细胞之间存在着一定的异速生长关系，因此通常情况下，耳石的宽度和鱼类的体长之间存在一定的线性关系[33-34]，其年轮宽度或增量对鱼类的生长动态具有足够的代表性。需要注意的是，适用于数据有限场景的资源评估方法虽然可以借助少量且易于获得的数据对鱼类种群进行评估，但其过程大多依赖于算法求解，其结果可能受到算法随机性的影响[1]，而年轮信息的真实性可以保证评估结果的稳健。另一方面，回归模型可以较为直观地分析耳石增量的年际差异，但研究耳石增量对环境因素的响应时，回归模型只能用于分析耳石增量与诸多因素间在统计学上的关系，并不能给出在生物学上的原因及内在机制，因此在应用过程中需提前了解所研究鱼种的生物学特征，以保证研究结果在生物学上的合理性。

随着气候变化、过度捕捞和生境破坏等多重因素的影响，许多鱼类的生物学特征发生重要变化。本研究结果表明，海州湾方氏云鳚具有较为明显的生长变快的现象，相关研究显示黄海海域中如鳀鱼（Engraulis japonicus）等鱼种也存在生长变快的现象[35]，因此掌握鱼类对环境因素的响应是应对未来气候变化的重要措施。目前世界绝大多数渔业属于数据缺乏渔业[36]，数据有限的资源评估方法成为了目前渔业资源评估的主要手段。本研究强调应当重视传统评估方法的精确性，在不同的场景应用不同的方法，充分利用不同方法各自的优势，以更好地应对未来气候变化对鱼类群体以至于生态系统的影响。