体长-体重关系(length-weight relationship, LWR)是鱼类生活史研究的基础内容, 能够通过体长计算体重，进而求得种群生物量, 可以反映鱼类的生长模式(等速还是异速生长, 条件因子如何), 也可以间接反映鱼类栖息环境的质量^[1]。因此体长-体重关系的研究是鱼类资源评估模型中不可或缺的部分, 对资源评估和管理具有重要意义^[2-3]。

大量研究显示, 鱼类的体长-体重关系会受到季节、地理和环境条件等因子的影响, 因而导致其关系式中的参数值存在不同的异质性。例如, Moutopoulos等^[4]发现爱琴海中九带 (Serranus cabrilla)的体长-体重关系(W=aL^b)中的a值和b值存在季节性差异, a值在春季最大, 其次是冬季和夏季, 秋季最小, b值的规律与此相反; Li等^[5]研究得出, 渤海和黄海的小黄鱼(Larimichthys polyactis)体长-体重关系有水域和年际上的变化, 北黄海和渤海种群的b值在1982—1993年呈增长趋势, 而南黄海种群的b值在1986—1994年逐渐下降; 林斌和徐青^[6]在对大亚湾的带鱼(Trichiurus lepturus)、二长棘鲷(Paerargyrops edita)和蓝点马鲛(Scomberomorus niphonius)等鱼类的研究中发现, 同一种类的鱼在不同的生态环境下, 其体长-体重关系的b值也不同, 且比东海和黄渤海的值高。

在不同的异质性研究方面, 线性混合效应模型(linear mixed-effects model, LMEM)是一种应用极其广泛的统计学模型。其线性预测不仅包含固定效应, 也包含随时间和地点等因子变化产生的随机效应, 因此能够分析多样化的数据^[7], 而且已经广泛应用在很多领域中。Baayen等^[8]介绍了混合效应模型能估计各种相关因子产生的随机效应; Xu等^[7]用线性混合效应模型研究福建省杉木(Cunninghamia lanceolata)的冠幅, 得出了一个包含由地点引起的随机效应的模型; Ma等^[1]在中国北部沿海小黄鱼体长-体重关系的研究中使用了线性混合效应模型, 在模型中考虑了不同年份和海域下的随机效应。

黄(Lophius litulon)是中国渤海、黄海和东海的一种重要经济鱼类, 近年来其捕捞量逐年增多^[9-10]。随着海洋渔业资源衰退日渐严重, 对该经济种类资源评估的研究以及科学管理也越来越紧迫。而黄资源评估所需的基础生物学研究, 尤其是其时空差异方面的研究并不多。李忠炉等^[11]研究了黄海中南部黄的体长组成、体长-体重关系、性比和性腺成熟度等生活史特征的年际变化, 得出LWR的各参数存在年间差异。林龙山等^[12]描述了2002年东海区黄的体长组成、体重组成、体长-体重关系、繁殖习性和摄食习性等生物学特征, 徐开达等^[13]研究了东海北部、黄海南部黄的年龄与生长。而本文根据在山东近海进行的渔业资源底拖网调查数据, 对黄的体长-体重关系进行研究, 并运用线性混合效应模型分析了其时空差异, 以期为该种类的资源评估和渔业管理提供数据支撑。

黄样品取自在山东近海海域(图 1)进行的秋(2016年10月)、冬(2017年1月)、春(2017年5月)和夏(2017年8月)4个季节的底拖网调查。该调查按季节分4个航次, 每个航次的调查站位相同, 均设118个站位, 其中有116个站位在至少一个航次中采集到黄样品(表 1)。根据纬度和水深, 将调查海域分为6个水域(图 1), 其中A、B和C区位于36°48′N以北, D、E和F区位于36°48′N以南。A和D区处于水深30 m以浅的水域, B和E区处于30~50 m水深之间, C和F区处于水深50 m以深的水域。调查船为220 kW的单拖渔船, 底拖网为1380目× 4.67 cm。拖网时网口高度约7.53 m, 宽约15 m。参照《海洋调查规范》中的要求, 样品捕捞上来后, 先进行冷冻保存, 然后带回实验室进行生物学测定, 其中体长和体重分别测量到0.1 cm和0.01 g。

图 1(Fig. 1)

图 1 山东近海黄的调查区域 A和D区处于水深30 m以浅的水域, B和E区的水深为30~50 m, C和F区处于水深50 m以深的水域. Fig.1 Survey map of Lophius litulon in Shandong coastal waters The depth of A and D areas is less than 30 m, the depth of B and E areas is 30-50 m, and the depth of C and F areas is greater than 50 m.

表 1(Tab.1)

表 1 山东近海黄在不同季节和水域的样本数量 Tab.1 Sample size of Lophius litulon in different seasons and regions in Shandong coastal waters 水域 region 秋季 autumn 冬季 winter 春季 spring 夏季 summer 总和 total A - 16 238 26 280 B 33 95 83 72 283 C 41 34 33 4 112 D 21 55 50 11 137 E 23 100 12 51 186 F 128 90 33 167 418 总和total 246 390 449 331 1416 注: A和D区处于水深30 m以浅的水域, B和E区的水深为30~ 50 m, C和F区处于水深50 m以深的水域, 详见图 1. Note: The depth of A and D areas is less than 30 m, the depth of B and E areas is 30~50 m, and the depth of C and F areas is deeper than 50 m, as shown in Fig. 1 for details.

表 1 山东近海黄在不同季节和水域的样本数量 Tab.1 Sample size of Lophius litulon in different seasons and regions in Shandong coastal waters

黄的体长-体重关系用幂指数方程表示:

式中, W为体重(g), L为体长(cm), a是条件因子。b是生长参数, 若b < 3, 则表示随着体长的增加, 鱼的肥满度下降, 身体趋于细长; 若b=3, 则为等速生长; 若b > 3, 则表示随着体长的增加, 鱼的肥满度升高, 体高和体宽也增大^[15]。

本文使用广义线性模型(generalized linear model, GLM)和9个LMEM模型来描述体长-体重关系(表 2), 其中这9个LMEM模型在参数a或b是否包含水域或季节产生的随机效应上各不相同。LMEM模型将水域和季节的影响作为参数a和b的随机效应, 以阐明体长和体重关系的时空差异^[7]。使用F检验评估时空因子对黄体长-体重关系的显著性。所有模型的构建过程都是用R软件的“lme4”程序包来完成^[16]。

表 2(Tab.2)

表 2 山东近海黄体长-体重关系模型的拟合效果指标(AIC)和交叉验证的评价指标(RMSE和R2) Tab.2 Fitting effect index (AIC) and cross validation evaluation indices (RMSE and R2) of length-weight relationship models of Lophius litulon in Shandong coastal waters 模型缩写 model abbreviations 模型model 对数转化log-transformed AIC RMSE R2 GLM W=a×L^b ln(W)=ln(a)+b×ln(L) 784 1.348 0.002 T.I W=[a×exp(ReT.I)]×L^b ln(W)= [ln(a)+ReT.I]+b×ln(L) 758 0.308 0.898 T.S W=a×L^(b+ReT.S) ln(W)=ln(a)+(b+ReT.S)×ln(L) 760 0.304 0.899 T.I & S W=[a×exp(ReT.I)]×L^(b+ReT.S) ln(W)=[ln(a)+ReT.I]+(b+ReT.S)×ln(L) 747 0.318 0.891 R.I W=[a×exp(ReR.I)]×L^b ln(W)=[ln(a)+ReR.I]+b×ln(L) 679 0.31 0.894 R.S W=a×L^(b+ReR.S) ln(W)=ln(a)+(b+ReR.S)×ln(L) 706 0.304 0.901 R.I & S W=[a×exp(ReR.I)]×L^(b+ReR.S) ln(W)=[ln(a)+ReR.I]+(b+ReR.S×ln(L) 468 0.29 0.906 R & T.I W=[a×exp(ReR.I)×exp(ReT.I)]×L^b ln(W)=[ln(a)+ReR.I+ReT.I]+b×ln(L) 666 0.31 0.898 R & T.S W=a×L^(b+ReR.S+ReT.S) ln(W)=ln(a)+(b+ReR.S+ReT.S)×ln(L) 695 0.299 0.903 R & T.I & S W=[a×exp(ReR.I)×exp(ReT.I)]× L^(b+ReR.S+ReT.S) ln(W)=[ln(a)+ReR.I+ReT.I]+ (b+ReR.S+ReT.S)×ln(L) 439 0.278 0.913 注:第一列为第二、三列模型的缩写.T.I、T.S、T.I & S、R.I、R.S、R.I & S、R & T.I、R & T.S和R & T.I & S分别为截距ln(a)的季节随机效应、斜率b的季节随机效应、截距ln(a)与斜率b的季节随机效应、截距ln(a)的水域随机效应、斜率b的水域随机效应、截距ln(a)与斜率b的水域随机效应、截距ln(a)的季节与水域随机效应、斜率b的季节与水域随机效应，以及截距ln(a)与斜率b的季节与水域随机效应。 Note: The first column shows the abbreviations of models detailed in the second and third columns. T.I, T.S, T.I & S, R.I, R.S, R.I & S, R & T.I, R & T.S和R & T.I & S indicate random effects on intercept [ln(a)] from seasons, random effects on slope(b) from seasons, random effects on intercept [ln(a)] and slope(b) from seasons, random effects on intercept [ln(a)] from regions, random effects on slope(b) from regions, random effects on intercept [ln(a)] and slope(b) from regions, random effects on intercept [ln(a)] from seasons and regions, random effects on slope(b) from seasons and regions, and random effects on intercept [ln(a)] and slope(b) from seasons and regions, respectively.

表 2 山东近海黄体长-体重关系模型的拟合效果指标(AIC)和交叉验证的评价指标(RMSE和R2) Tab.2 Fitting effect index (AIC) and cross validation evaluation indices (RMSE and R2) of length-weight relationship models of Lophius litulon in Shandong coastal waters

利用赤池信息量准则(Akaike Information criterion, AIC)比较10个模型的拟合程度, 其值越小, 则表示模型的拟合效果越好^[17]。为了评估模型的预测性能, 本研究应用交叉验证, 随机抽取80%的数据构建模型, 然后用建立的模型为剩余20%的数据进行预测, 将预测值和真实值进行对比检验, 重复100次, 最后得出该模型预测性能的评估结果^[18]。用均方根误差RMSE (root mean square error)和决定系数(R²)判断最适合的模型。均方根误差越接近于0, 决定系数越接近于1, 则表示模型的预测值与样品的真实值越吻合, 即预测性能越好。

2.1 黄的体长和体重分布

山东近海此次调查所得黄样品的体长范围为1.1~48.3 cm(图 2), 平均体长为22.67 cm, 优势体长组为20~30 cm。体重范围为4.19~3478.67 g, 平均体重为339.54 g, 优势体重范围为0~500 g。大部分春季采集的黄个体较大(图 3), 其次是冬季和秋季个体, 夏季调查到的黄个体最小。分布在水深30 m以浅的水域黄个体较大, 其次是分布在水深30 m和50 m之间水域的个体, 而水深50 m以深的水域个体最小。同等水深的海域, 北部的个体较大, 即大部分黄在A、B和C区分别比在D、E和F区的个体偏大。

图 1(Fig. 1)

图 2 山东近海黄体长和体重分布 Fig.2 Body length and weight distributions of Lophius litulon in Shandong coastal waters

图 1(Fig. 1)

图 3 山东近海黄在不同季节和水域的体长和体重分布 Fig.3 Seasonal and regional variations of body length and body weight of Lophius litulon in Shandong coastal waters

在10个模型中, LMEM(R & T.I & S)的AIC值最小(表 2), 为439, 即同时具有水域和季节对两个参数的随机效应的线性混合效应模型对黄体长-体重关系的拟合效果最好。交叉验证的结果也表明此模型的RMSE最接近于0, R²最接近于1, 预测效果最好(表 2)。在这个模型中, a和b的固定效应估计值分别是0.0477和2.7729。包含季节和水域随机效应的a估计值范围为0.0089~ 0.9384, b估计值范围为1.8773~3.2813(图 4 A)。

图 1(Fig. 1)

图 4 LMEM(R & T.I & S)中a和b的季节和水域差异 LMEM(R & T.I & S)为同时具有季节和水域对参数a和b的随机效应的线性混合效应模型. A:具有季节和水域随机效应的a值和b值; B: a和b的时间变化; C: a和b的空间变化. Fig.4 Variations of a and b among seasons and regions from the LMEM(R & T.I & S) LMEM (R & T.I & S) is a linear mixed-effects model with random effects on parameter a and b from seasons and regions. A: estimates of parameter a and b with random effects from seasons and regions; B: temporal variations of a and b; C: spatial variations of a and b.

LMEM(R & T.I & S)模型的结果显示, 黄体长-体重关系具有明显的时空差异。F检验结果表明, LMEM(R & T.I & S)和GLM两个模型差异极显著(P < 0.01), 即季节和水域对黄体长-体重关系具有极显著影响。在LMEM(R & T.I & S)模型中, 从不同季节来看, a值在春季最大, 其次是秋季和冬季, 而夏季最小; b值在夏季最大, 其次是冬季和秋季, 而春季最小。从不同水域来看, a值的大小顺序是B > C > D > A > F > E, b值与此相反。整体来看, 纬度高的a值较大, b值较小, 而a、b值与水深的变化没有明显规律。a值和b值的水域差异范围分别为0.3435和0.9637, 而a值和b值的季节差异范围分别为0.0974和0.4402, 说明水域差异比季节差异大(图 4)。

在体长0~30 cm的黄中(图 5), 体重的季节和水域差异较小; 在体长大于30 cm时, 体重的时空差异较大。黄的体重增长率在春季与其他季节明显不同, 而在夏季和冬季却十分相似。B区黄的体重增长率与其他地区也明显不同, 而A区和F区的曲线几乎重合。在30~50 cm的体长下, 各季节的体重增长率大小顺序是:夏 > 冬 > 秋 > 春; 各水域的体重增长率大小顺序是: E > F > A > D > C > B。该结果与b值的时空差异趋势是一致的。

图 1(Fig. 1)

图 5 山东近海黄在不同季节(A)和水域(B)的体长-体重关系差异 Fig.5 Variations of length-weight relationships for Lophius litulon among seasons (A) and regions (B) in Shandong coastal waters

3.1 黄体长-体重关系的时空差异

根据LMEM(R & T.I & S)的结果, 由4个季节和6个水域得出的24个a的估计值(0.0089~0.9384)中, 有13个a值的范围为0.0089~0.0404, 在大多数鱼类的共同范围(0.001~0.05)之内^[3], 也在Fishbase^[19]中黄的a值范围(0.00804~0.04119)之内。本研究中估计值偏高的a值出现在B、C区的各个季节, 以及春季的几乎所有水域(只有E区的a值在范围内, 但也偏高, 达到0.0404)。根据相同体长的鱼体重越大, 则表示处于更好的生长条件之中, 所以条件因子可以用来比较鱼的生长条件、肥满度或健康程度^[20]。上述结果可能是因为黄存在北黄海和南黄海群系上的差别, 北黄海水温较低, 黄长得慢, 但是个体更大, 且主要分布于水深超过 40~50 m 的海域^{[11-12, 21]}, 反映出B、C区的黄具有更好的生长条件, 所以a值较大; 黄主要是在春季进行产卵^[22-23], 此时个体的肥满度更大, 所以可能造成黄表现出较大的a值。另外, 在李忠炉等^[11]对于黄海中南部1985年、2000年、2005年和2009年秋季黄的研究中, a的范围为0.0219~0.6524, 也出现了较大的a值。

Fishbase^[20]中黄的b值为 2.96(2.75~3.17), Carlander^[24]和 Froese^[3]的研究显示大部分鱼类b值通常分别在 2.5~3.5 和 2.7~3.4 的范围之内, 与 本研究得出的b估计值(1.8773~3.2813)大体上一致。异速生长参数b的固定值(2.7729)小于3, 表示黄为负异速生长, 随着体长的增加, 黄的肥满度降低, 身体朝细长方向发展。本研究小于2.5的b值出现于所有季节的B区和春季的C区, 且春季所有水域的b值均小于3, 春季B区的b值(1.8773)甚至小于2, 反映出B区和春季的黄体型偏于细长。

在本研究中, a值在总体上表现出北部大于南部。在水深方面, a、b值没有呈现出明显的规律, 这可能是C和D区样品的数量较少和体长范围较小造成的。本研究共有1416个样品, 体长范围为1.1~48.3 cm, 而C和D区分别只有112和137个样品, 体长范围分别为15.0~42.5 cm和7.0~ 39.5 cm。Froese^[3]的研究显示, 小的样品范围和少的数据会导致体长-体重关系参数值的结果变化性较大。因此, C和D区可能出现了异常的结果, 对整体结果的分析产生了干扰。

根据AIC值比较和交叉验证的结果可知, 任何一个混合效应模型的拟合效果和预测能力, 均优于简单的广义线性模型。而包含对参数a和b的水域和季节随机效应的LMEM(R & T.I & S)模型拟合效果最好, 所以用于研究黄体长-体重关系的时空差异。Petrakis等^[25]、Bolognini等^[26]、Fatih等^[27]和Wilson等^[28]的研究显示, 鱼类的体长-体重关系会随水域和时间而变化, 本研究中的线性混合效应模型也表现出这种差异。

由不同的数据来源建立不同的模型, 然后对这些模型的参数估计值进行比较, 这是异质性研究的传统方法。例如, Liu等^[29]为了探求海州湾小黄鱼的体长-体重关系随时间的变化, 根据在5月、7月、9月和12月采集的数据, 构建了4个LWR模型; 管卫兵等^[30]由2006—2008年的调查数据构建了各年的LWR模型, 以研究长江口水域刀鲚(Coilia nasus)溯河群体体长-体重关系的年际变化; 牛玉娟等^[31]根据伊犁河三条支流采集的数据, 分别构建了一个LWR模型, 研究了各支流中新疆裸重唇鱼(Gymnodiptychus dybowskii)的体长-体重关系。这些单独构建起来的模型得出的参数估计值, 需要根据差异性检验进一步比较, 而随机效应的使用可以非常有效、合理地进行异质性研究。最常用的随机效应应用方法是混合效应模型, 在反映生物学时空差异的研究中更为普遍, 是一种将固定效应和随机效应结合起来的特殊统计模型^[32-35]。在该模型中, 固定效应能反映所有样品的总体特点, 而随机效应则因样品不同的来源, 呈现出拟合参数值与参数平均值(即固定效应)的差异^[36-38]。随机效应的估计值能反映多种数据来源下体长-体重关系的参数分布, 也可以避免统计分析中的过度参数化问题^{[34, 39-40]}。

混合效应模型可以将所有数据来源产生的效应包含进一个模型之中。许多研究使用某一特定水域或时间的LWR来表现鱼类体长-体重关系的时空差异^[41-43], 而混合效应模型能把水域和时间的随机效应在单个模型中体现出来, 从而能更快捷、有效地估计LWR的时空差异, 进一步证实了此模型在数据来源异质性研究中的优势。由于本研究只考虑了季节和水域的随机效应, 存在一定的不足, 所以在未来的研究中, 应该充分考虑性别、年份、环境因子和生长阶段等因素产生对黄生长的影响, 更全面地研究黄生活史特征的变化。