基于混合效應(yīng)模型的海州灣方氏云鳚的生長動態(tài)及其影響因素分析

趙天亞 ，楊曉改，李勛，薛瑩 ，徐賓鐸 ，紀毓鵬 ，任一平 ，張崇良 *

( 1. 中國海洋大學 水產(chǎn)學院，山東 青島 266003；2. 青島海洋科學與技術(shù)試點國家實驗室 海洋漁業(yè)科學與食物產(chǎn)出過程功能實驗室，山東 青島 266237；3. 海州灣漁業(yè)生態(tài)系統(tǒng)教育部野外科學觀測研究站，山東 青島 266003；4. 中國農(nóng)業(yè)出版社，北京 100125；5. 青島西海岸新區(qū)海洋發(fā)展局，山東 青島 266400)

1 引言

魚類的生長特征是漁業(yè)資源評估的重要內(nèi)容，可以反映魚類群體的生存環(huán)境和資源狀況。根據(jù)可用的數(shù)據(jù)材料可知，魚類生長的研究方法主要有兩種：基于長度的評估方法和基于年齡鑒定的評估方法?；陂L度的方法如電子體長頻率分析法（Electronic Length Frequency Analysis, ELEFAN）和MULTIFAN等適用于年齡數(shù)據(jù)難以獲得、長度數(shù)據(jù)易于獲得的數(shù)據(jù)有限的場景，但易受到采樣隨機性和算法穩(wěn)健性的影響[1]。而魚體硬組織材料如耳石和鱗片等可以通過年輪的形式記錄魚類的生長信息，其輪紋增量可以直接反映魚類在該年份的生長速度，從而直觀地反映魚類的生長動態(tài)，為探討影響魚類生長的因素提供了重要材料。

本研究以方氏云鳚（Pholis fangi）為例研究魚類的生長特征。方氏云鳚作為一種小型餌料魚類，在食物網(wǎng)中維持關(guān)鍵的能流傳遞過程，在生態(tài)系統(tǒng)中發(fā)揮重要作用[2]。伴隨著近海漁業(yè)資源的明顯衰退，方氏云鳚等小型底棲魚類逐漸成為海州灣及其鄰近海域的主要優(yōu)勢種[2]。近年來，受到全球氣候變化及人類活動的影響，許多魚種在生理、行為及分布等方面上均出現(xiàn)多種適應(yīng)性轉(zhuǎn)變[3]。盡管已有研究估算了方氏云鳚的生長參數(shù)[4]，但在氣候變化背景下，其生長特征可能發(fā)生變化，生長的動態(tài)及驅(qū)動因素目前尚缺乏探討。

線性混合效應(yīng)模型（Linear Mixed-Effects Model,LMEM）是一種應(yīng)用廣泛的統(tǒng)計學模型[5]，通過在一般線性模型中加入隨機效應(yīng)，可以有效解釋數(shù)據(jù)中的嵌套關(guān)系對整體線性回歸所產(chǎn)生的影響。該模型在漁業(yè)領(lǐng)域亦有廣泛應(yīng)用，如Ma等[6]通過線性混合效應(yīng)模型研究了黃、渤海海域小黃魚（Larimichthys polyactis）體長體重關(guān)系的時空異質(zhì)性；Mu等[7]使用線性混合效應(yīng)模型研究了黃海海域星康吉鰻（Conger myriaster）的生長動態(tài)。

為掌握方氏云鳚的生長動態(tài)及其對環(huán)境因子的響應(yīng)，本研究基于2015-2019年秋季海州灣方氏云鳚耳石樣本，利用線性混合效應(yīng)模型研究方氏云鳚生長速率的年際變化，并分析了影響其生長動態(tài)的主要因素。本研究旨在加深對方氏云鳚生物學特性的認識，為資源評估和漁業(yè)管理提供科學支撐，為理解氣候變化驅(qū)動下魚類的生物學特征響應(yīng)提供參考。

2 材料與方法

2.1 數(shù)據(jù)來源

本研究所用耳石樣本全部來自于2015-2019年秋季于海州灣及其鄰近海域進行的漁業(yè)資源底拖網(wǎng)調(diào)查，調(diào)查海域范圍為34°20′～35°40′N，119°20′～121°10′E，站位設(shè)計采用分層隨機采樣的方法（圖1），每個航次設(shè)置18個站位，在各區(qū)域中隨機選取站位[8]。所用調(diào)查船為功率220 kW的單拖底拖網(wǎng)漁船，拖網(wǎng)速度為2～3 kn，每站拖網(wǎng)時間為1 h。調(diào)查網(wǎng)具網(wǎng)目為17 mm，網(wǎng)口擴張高度約為6 m，寬度約為25 m。將所采集樣本帶回實驗室進行解剖，從每條方氏云鳚樣本的左右耳石中隨機選取一枚進行分析。樣本的采集和分析參考《海洋調(diào)查規(guī)范》（GB/T 12763.6-2007）[9]。

圖1 海州灣及鄰近海域調(diào)查區(qū)域Fig. 1 Survey areas in Haizhou Bay and its adjacent waters

本研究共選取107枚方氏云鳚耳石樣本，其中2015年共選取22枚樣本，2016年共選取24枚，2017年共選取21枚，2018年共選取19枚，2019年共選取21枚。使用Olympus bx51型光學顯微鏡采集耳石樣本圖像并讀取年齡，結(jié)合樣本年齡和采集年份反推其每一輪紋的生成年份，再沿耳石長軸從耳石核心起始對每個年輪進行測距。讀耳石過程遵守盲讀和重復性原則[7]，年輪測距過程通過imageJ軟件實現(xiàn)。本研究所選樣本包含1～5齡個體，但5齡個體較少，第5年輪的耳石增量數(shù)據(jù)不具有代表性，因此在數(shù)據(jù)整理中將其剔除。經(jīng)整理篩選后共得到海州灣方氏云鳚的241個耳石增量數(shù)據(jù)，其輪紋形成時間為2013-2018年。

2.2 因子篩選

本研究以耳石增量作為響應(yīng)變量，以年齡、出生年份、底層海水溫度、葉綠素含量和種群密度等為解釋變量，構(gòu)建回歸模型。方氏云鳚為一種底棲小型餌料魚類，營養(yǎng)級較低，因此本研究以海水葉綠素含量反映其食餌條件；同時以輪紋形成時間所對應(yīng)年份的底拖網(wǎng)調(diào)查結(jié)果中的站位平均標準化數(shù)量反映其密度制約強度，其計算方式為：站位平均標準化數(shù)量=∑站位標準化數(shù)量/總站位數(shù)，單位為ind./站。本研究中葉綠素含量數(shù)據(jù)及底層溫度數(shù)據(jù)提取自FVCOM（Finite-Volume Community Ocean Model）[10]。以方差膨脹因子（Variance Inflation Factor，VIF）作為因子多重共線性衡量指標對初始因子進行篩選，通常認為VIF＞10即存在多重共線性問題，應(yīng)在建模前予以剔除。

2.3 線性混合效應(yīng)模型

本研究使用線性混合效應(yīng)模型研究海州灣方氏云鳚耳石增量的年際差異。線性混合效應(yīng)模型包括固定效應(yīng)和隨機效應(yīng)兩部分，固定效應(yīng)用于解釋線性回歸的整體期望，隨機效應(yīng)用于解釋數(shù)據(jù)中嵌套關(guān)系對整體期望的影響。線性混合效應(yīng)模型的一般形式為

式中，Yi為響應(yīng)變量；α為模型截距；Xi×β為固定效應(yīng)；Zi×bi為隨機效應(yīng)；εi為誤差項。本研究根據(jù)加入模型的影響因素來源的不同，分別構(gòu)建了內(nèi)因模型和外因模型。

2.3.1 內(nèi)因模型

內(nèi)因模型考慮了影響生長的種群內(nèi)在因素。其中年齡通常為生長速度的主要決定因素，因此將年齡作為固定效應(yīng)加入模型，此外將個體編號、年份和世代作為隨機效應(yīng)加入模型，模型響應(yīng)變量為耳石增量。使用小樣本校正的赤池信息準則（AICc）作為模型選擇標準[11]，AICc最小的模型為最優(yōu)模型。使用極大似然估計法重新擬合最優(yōu)模型，以得到耳石增量在年際間差異的無偏估計值[5]。

2.3.2 外因模型

外因模型考慮了外在環(huán)境因子對方氏云鳚的生長速度的影響，模型響應(yīng)變量為耳石增量，解釋變量包括底層溫度、葉綠素含量和種群密度。模型中加入因子二次項，以反映環(huán)境因子與生長速度的非線性關(guān)系。由于方氏云鳚在幼體時期（小于1齡）生長速度顯著大于成體時期（大于1齡），其生活史特征差異較大[12]，因此本研究以0齡時期所對應(yīng)耳石增量（第1年輪）反映方氏云鳚幼體時期的生長特征，以1齡時期所對應(yīng)耳石增量（第2年輪）反映方氏云鳚成體時期的生長特征，分別對方氏云鳚0齡時期及1齡時期的耳石增量對環(huán)境因素的響應(yīng)進行建模分析。模型擬合過程采用逐步回歸法，根據(jù)AICc最小準則分別選取最優(yōu)模型。

本研究中數(shù)據(jù)分析、模型擬合及繪圖過程均于R語言4.0.2中完成，其中LMEM的擬合由lme4軟件包實現(xiàn)。

3 結(jié)果

3.1 耳石增量分析

方氏云鳚個體在不同年齡的耳石增量具有明顯差異（圖2），第1年輪的平均耳石增量為0.327 mm，顯著高于高齡耳石增量。第1年輪的耳石增量在2013-2018年呈現(xiàn)先上升后下降的趨勢，在2016年達到最大值（0.341±0.024）mm；第2至第4年輪的耳石增量均呈現(xiàn)不規(guī)則波動趨勢。此外，方氏云鳚在不同年齡階段的耳石增量波動范圍變化較大，變異系數(shù)隨年齡逐漸增大，表明個體間生長變異水平隨年齡逐漸提高（表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 內(nèi)因模型分析

通過調(diào)整隨機效應(yīng)的不同結(jié)構(gòu)，共擬合了3個內(nèi)因模型（表2）。其中模型M3的AICc值最低，為最優(yōu)內(nèi)因模型，包含個體編號和年份兩個隨機效應(yīng)，共解釋了85.9%的偏差。年齡固定效應(yīng)表明，方氏云鳚的生長速度隨著年齡的增加而顯著減緩（圖3a）。年份隨機效應(yīng)表明，方氏云鳚耳石增量在2016年達到最大值（0.260±0.004）mm，其生長速度在2013-2016年呈現(xiàn)逐步變快的趨勢，在2016-2018年呈現(xiàn)波動趨勢（圖3b）。

圖3 內(nèi)因模型的年齡與年份效應(yīng)Fig. 3 Age and Year effects of the intrinsic model

表2 內(nèi)因模型的擬合Table 2 Fitting of the best intrinsic model

3.3 外因模型分析

多重共線性檢驗結(jié)果表明，底層溫度、葉綠素含量和種群密度的VIF均小于10（表3）。因此3個環(huán)境因子均可作為初始解釋變量。通過逐步回歸法分別擬合最優(yōu)0齡模型（Model 1）和1齡模型（Model 6），兩個最優(yōu)模型均只包含全年平均底溫（BT）（表4），分別解釋了22.9%和2.5%的偏差。本研究將全年平均底溫分離為春季平均底溫、夏季平均底溫、秋季平均底溫和冬季平均底溫，并將4種底溫分別作為初始解釋變量加入模型，以評估不同季節(jié)的溫度對其生長影響的強弱并改良模型的擬合效果。

分別將四季底溫作為初始變量時，3個環(huán)境因子的VIF均小于10（表3）。將季節(jié)底溫分別作為模型初始解釋變量并重新擬合模型，最優(yōu)模型篩選結(jié)果見表4。在0齡模型中，季節(jié)底溫模型（Model 2～5）分別解釋了24.7%～30.1%的偏差，其AICc值均小于Model 1；在1齡模型中，季節(jié)底溫模型（Model 7～10）分別解釋了0.9%～2.9%的偏差，其AICc值均大于Model 6。

表3 因子多重共線性檢驗結(jié)果Table 3 Results of multicollinearity test of factors

表4 最優(yōu)外因模型篩選結(jié)果Table 4 Results of best extrinsic models screening

3.4 環(huán)境因素對耳石增量的影響

最優(yōu)模型擬合結(jié)果表明，環(huán)境因素的加入可以優(yōu)化0齡模型的擬合及解釋效果（表4）。在0齡模型（Model 1～5）中，所加入的環(huán)境因子對0齡耳石增量均具有顯著影響效應(yīng)（p＜0.05）；而在1齡模型中（Model 6～10）中，環(huán)境因子對1齡耳石增量均無顯著影響（p＞0.05）。整體而言，方氏云鳚0齡時期的生長速度受到環(huán)境因子的顯著影響，而1齡時期的生長速度則主要受到自身年齡的制約。

在全年平均底溫模型Model 1和Model 6中，底溫每升高1℃，會導致0齡耳石增量增加14.49%，1齡耳石增量減少13.83%（圖4），因此方氏云鳚在0齡及1齡時期的生長速度對底溫的響應(yīng)存在一定差異。

圖4 全年平均底層溫度對方氏云鳚耳石增量的影響Fig. 4 Effect of annual average bottom temperature on otolith increment of Pholis fangi

在季節(jié)底溫模型Model 2～5中，四季平均底溫對0齡耳石增量的影響效應(yīng)存在一定差異，0齡耳石增量隨著春季底溫和冬季底溫的升高而增加，隨著夏季底溫和秋季底溫的升高而先增加后減少（圖5）。種群密度對0齡耳石增量存在一定的負效應(yīng)，耳石增量隨著種群密度的上升而減少（圖5）。葉綠素含量對0齡耳石增量的影響效應(yīng)在模型間存在差異，在模型Model 2和Model 5中，0齡耳石增量隨葉綠素含量的上升而減少；在模型Model 4中，0齡耳石增量隨葉綠素含量的上升而增加（圖5）。

圖5 環(huán)境因子對方氏云鳚0齡耳石增量的影響Fig. 5 Effect of environmental factors on otolith increment of Pholis fangi at age 0

4 討論

本研究基于耳石增量數(shù)據(jù)，利用混合效應(yīng)模型研究了海州灣方氏云鳚在2013-2018年生長速度的年際差異及其生長對外界環(huán)境因子的響應(yīng)。結(jié)果表明，海州灣方氏云鳚存在一定的生長變異，同時該群體在2013-2018年內(nèi)的生長速度整體呈現(xiàn)變快的趨勢；方氏云鳚在0齡時期的生長速度受環(huán)境因子的影響較大，受到底層溫度、食餌條件和密度制約的調(diào)控；而年滿1齡后其生長速度受環(huán)境因子變化的影響較小，主要受到自身年齡的限制。

在2013-2018年內(nèi)，方氏云鳚的生長速度具有明顯的變化，2016-2018年的生長速度相較于2013-2015年更快，這與王琨等[4]使用ELEFAN方法所評估的結(jié)果較為一致。方氏云鳚生長變快的趨勢可能與目前的人類活動及氣候變化的大環(huán)境具有密切聯(lián)系。在適溫范圍內(nèi)，溫度的適當提升可以促進魚類的新陳代謝和攝食，進而加快魚類的生長速度，而Gao等[13]的研究表明，2016-2018年黃海南部出現(xiàn)較大程度的海洋變暖，這與本研究中2016-2018年方氏云鳚生長速度較快的現(xiàn)象較為一致。此外，海州灣海域長期的過度捕撈也會使魚類產(chǎn)生適應(yīng)性轉(zhuǎn)變[14]，外界壓力可能使方氏云鳚選擇加快生長速度的策略以快速到達性成熟。值得注意的是，方氏云鳚生長速度的年際差異主要來源于0齡時期生長速度的差異，而年滿1齡后的生長速度波動不大，這可能與其生活史特征有關(guān)。方氏云鳚為一種小型短周期魚類，在0齡時期生長速度極快，而性成熟后的生長速度大幅降低，資料顯示方氏云鳚年滿1齡時的體長即可達到100 mm[15]。通常而言，魚類在幼體時期的生長和生理狀態(tài)對外界環(huán)境較為敏感，因此其生長速度容易因外界環(huán)境因素的變化而產(chǎn)生較大程度的變化；而成體生長緩慢，其生長速度較為穩(wěn)定，對環(huán)境的適應(yīng)能力較強，其生長過程不易受到外界環(huán)境因素的影響。

溫度對方氏云鳚0齡時期的生長速度具有顯著的正效應(yīng)，對1齡時期的生長速度具有一定程度的負效應(yīng)。方氏云鳚的新生群體在孵化后進入為期1年的索餌期，此時會廣泛分布于黃、渤海沿岸海域[15]，而黃海沿岸海域在夏、秋季節(jié)的底層水溫為18～22℃[16]。本研究中Model 1～5的結(jié)果表明，溫度的升高可以促進方氏云鳚0齡時期的生長，直至21.5℃左右時生長速度達到最快，符合其生活史特征，因此較高的溫度可以促進方氏云鳚0齡群體的生長，且該溫度可能為0齡幼體的最適生長溫度，在其生長過程中起到關(guān)鍵作用。底層溫度對方氏云鳚1齡時期的生長速度具有一定程度的負效應(yīng)，這可能與其冷溫習性有關(guān)。年滿1齡后的性成熟個體通常在產(chǎn)卵后游向深水區(qū)，并在夏、秋季節(jié)集中于青島冷水團區(qū)域[17]，在冬季重新返回產(chǎn)卵場進行產(chǎn)卵繁殖，其生活史特征與幼體階段具有較大差異。相較于0齡時期，方氏云鳚年滿1齡后的生長速度較為穩(wěn)定，環(huán)境因素的變化對其耳石增量的影響較小，其生長速度的變化主要來自于自身年齡限制。

密度制約也是影響方氏云鳚0齡時期生長速度的一個重要因素。密度制約是一個綜合性效應(yīng)，通常會與魚類棲息地的理化環(huán)境及捕食與競爭等生物過程密切相關(guān)[18]，普遍被認為是魚類補充群體的一個重要調(diào)節(jié)手段[19]。方氏云鳚為小型短周期魚類，其生活史策略通常傾向于r選擇，傾向于將更多的能量用于繁殖過程，因此冬季孵化后的新生群體將存在較大的密度制約效應(yīng)[20]，主要表現(xiàn)為對食餌等資源的競爭，其結(jié)果為隨種群密度的上升，方氏云鳚幼魚將具有逐漸降低的生長速度及逐漸升高的幼體死亡率。此外，當0齡幼魚經(jīng)歷為期1年的索餌和發(fā)育后，將于冬季與高齡魚群共同開始向產(chǎn)卵場集中，進行產(chǎn)卵繁殖過程。在此過程中，個體的生長過程因能量分配和個體間競爭而受到限制，同時方氏云鳚的集群習性會進一步擴大密度制約效應(yīng)對其生長速度的影響。因此，方氏云鳚在0齡時期會經(jīng)歷較大程度的資源競爭，其生長速度受到密度制約較為顯著的調(diào)控。此外，許多魚類都有不同程度的生長變異現(xiàn)象[21]，生長變異的程度也可以反映個體間對資源的競爭強度，如藍紋鱸（Gramma loreto）和 前 唇 霓 虹 蝦 虎 魚（Elacatinus prochilos）等魚類群體會因資源競爭而在個體間產(chǎn)生較為明顯的體長差異[22]。方氏云鳚群體內(nèi)部較大的競爭強度導致個體間資源利用程度的均勻性較低，從而使該群體產(chǎn)生較大程度的生長變異現(xiàn)象。

葉綠素含量對方氏云鳚0齡時期生長速度的影響效應(yīng)在模型間存在一定差異，其原因為溫度與葉綠素含量間可能存在一定的關(guān)系，在模型中加入不同季節(jié)的溫度作為解釋變量時，會導致葉綠素含量對耳石增量的解釋效果出現(xiàn)差異。通常認為食餌條件與魚類的生長速度間存在正相關(guān)關(guān)系，而葉綠素含量可以反映浮游生物的豐度變化，進而在一定程度上反映方氏云鳚的食餌條件，但浮游生物豐度可能受其他諸多外界因素如溫度或其他低營養(yǎng)級魚類捕食的影響[23-24]，從而導致使用葉綠素含量反映方氏云鳚食餌條件可能存在一定的局限性。此外，有研究指出營養(yǎng)級聯(lián)效應(yīng)可能會影響魚類的捕食者群體動態(tài)或餌料豐度，并結(jié)合密度制約效應(yīng)間接調(diào)控魚類的生長[25-26]。其具體影響機制還有待進一步探討。

年齡與生長是漁業(yè)資源評估與管理的重要研究內(nèi)容，相關(guān)研究不僅可以掌握魚類的生物學特性及其資源現(xiàn)狀，還可以為資源評估模型如單位補充量漁獲量（Yield Per Recruitment, YPR）模型[27-28]、生態(tài)系統(tǒng)模型如OSMOSE等提供基本的生物學參數(shù)[29]。耳石作為一種可以記錄魚類生長信息的材料，在國內(nèi)外已被廣泛應(yīng)用于對年齡與生長的研究[30-32]。耳石和魚類體細胞之間存在著一定的異速生長關(guān)系，因此通常情況下，耳石的寬度和魚類的體長之間存在一定的線性關(guān)系[33-34]，其年輪寬度或增量對魚類的生長動態(tài)具有足夠的代表性。需要注意的是，適用于數(shù)據(jù)有限場景的資源評估方法雖然可以借助少量且易于獲得的數(shù)據(jù)對魚類種群進行評估，但其過程大多依賴于算法求解，其結(jié)果可能受到算法隨機性的影響[1]，而年輪信息的真實性可以保證評估結(jié)果的穩(wěn)健。另一方面，回歸模型可以較為直觀地分析耳石增量的年際差異，但研究耳石增量對環(huán)境因素的響應(yīng)時，回歸模型只能用于分析耳石增量與諸多因素間在統(tǒng)計學上的關(guān)系，并不能給出在生物學上的原因及內(nèi)在機制，因此在應(yīng)用過程中需提前了解所研究魚種的生物學特征，以保證研究結(jié)果在生物學上的合理性。

隨著氣候變化、過度捕撈和生境破壞等多重因素的影響，許多魚類的生物學特征發(fā)生重要變化。本研究結(jié)果表明，海州灣方氏云鳚具有較為明顯的生長變快的現(xiàn)象，相關(guān)研究顯示黃海海域中如鳀魚（Engraulis japonicus）等魚種也存在生長變快的現(xiàn)象[35]，因此掌握魚類對環(huán)境因素的響應(yīng)是應(yīng)對未來氣候變化的重要措施。目前世界絕大多數(shù)漁業(yè)屬于數(shù)據(jù)缺乏漁業(yè)[36]，數(shù)據(jù)有限的資源評估方法成為了目前漁業(yè)資源評估的主要手段。本研究強調(diào)應(yīng)當重視傳統(tǒng)評估方法的精確性，在不同的場景應(yīng)用不同的方法，充分利用不同方法各自的優(yōu)勢，以更好地應(yīng)對未來氣候變化對魚類群體以至于生態(tài)系統(tǒng)的影響。