基于混合效應模型的海州灣方氏云鳚的生長動態及其影響因素分析

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

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

1 引言

魚類的生長特征是漁業資源評估的重要內容，可以反映魚類群體的生存環境和資源狀況。根據可用的數據材料可知，魚類生長的研究方法主要有兩種：基于長度的評估方法和基于年齡鑒定的評估方法?；陂L度的方法如電子體長頻率分析法（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]，數據有限的資源評估方法成為了目前漁業資源評估的主要手段。本研究強調應當重視傳統評估方法的精確性，在不同的場景應用不同的方法，充分利用不同方法各自的優勢，以更好地應對未來氣候變化對魚類群體以至于生態系統的影響。