劉逸文， 張崇良??， 昝肖肖， 孫 銘， 徐賓鐸， 薛 瑩， 任一平,3

(1. 中國海洋大學(xué)水產(chǎn)學(xué)院, 山東 青島 266003； 2. 山東大學(xué)(威海)海洋學(xué)院， 山東 威海 264209；3.青島海洋科學(xué)與技術(shù)試點國家實驗室 海洋漁業(yè)科學(xué)與食物產(chǎn)出過程功能實驗室, 山東 青島 266237)

長蛇鯔(Sauridaelongate)屬暖溫性近海底層魚類，在中國渤海至南海均有分布，是中國重要的經(jīng)濟(jì)魚種之一。長蛇鯔的相關(guān)研究主要集中于年齡、生長、食性、遺傳學(xué)等方面[1-3]，有關(guān)其資源分布的研究較少，影響其分布的關(guān)鍵因素尚不明確。許多研究表明，在個體發(fā)育的不同階段，魚類的分布模式會發(fā)生變化[1, 4]，成體和幼體的分布與影響因子的關(guān)系可能存在差異。因此探究成體與幼體的環(huán)境適應(yīng)性差異，有助于解析魚類的分布模式及變化，識別和保護(hù)重要的魚類棲息地[5]，對于漁業(yè)資源的合理利用與保護(hù)有著重要意義。

當(dāng)前物種分布的研究越來越多地依賴于統(tǒng)計學(xué)方法，利用統(tǒng)計模型評估物種分布與影響因子的關(guān)系，其中，廣義可加模型(Generalized additive model，GAM)是研究物種分布常用的一種統(tǒng)計學(xué)模型。廣義可加模型是廣義線性模型的非參數(shù)化擴(kuò)展，其優(yōu)點是能夠處理高維數(shù)據(jù)中響應(yīng)變量與解釋變量間的非線性、非單調(diào)關(guān)系[6]。在漁業(yè)資源研究應(yīng)用中，水溫、鹽度、水深等自變量對于資源密度等因變量的影響往往不是線性的，因此，GAM模型的應(yīng)用有十分重要的意義[7]，已被越來越廣泛地應(yīng)用于漁業(yè)資源研究[5, 7-10]。

本文根據(jù)2016—2017年在山東南部海域進(jìn)行的漁業(yè)資源底拖網(wǎng)調(diào)查數(shù)據(jù)，利用GAM方法重點分析了秋季長蛇鯔的資源分布及其影響因素，并探討了影響長蛇鯔幼體與成體分布差異的原因，旨在更深入的了解長蛇鯔幼體和成體的生態(tài)習(xí)性和分布規(guī)律，揭示其生活史差異導(dǎo)致的環(huán)境適應(yīng)性差異，為長蛇鯔資源的合理利用和養(yǎng)護(hù)提供科學(xué)依據(jù)。

1 材料與方法

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

研究數(shù)據(jù)來源于2016年10月、2017年1、5、8月在山東南部進(jìn)行的4個航次(見表1)的漁業(yè)資源調(diào)查，調(diào)查范圍為119°30′E～124°E、35°0′N～36°45′N，共設(shè)置63個站位(見圖1)。用于樣品采集的調(diào)查船為單拖底拖網(wǎng)漁船，主機(jī)功率為220 kW，平均拖速為3 kn，拖網(wǎng)時間1 h，拖網(wǎng)時網(wǎng)口高度約7.53 m，寬約15 m。將漁獲物冷凍保存，帶回實驗室內(nèi)進(jìn)行樣品鑒定和生物學(xué)測量。調(diào)查及分析均按《海洋生物調(diào)查規(guī)范》(GB 12763.6—2007)[11]要求進(jìn)行。

圖1 山東南部近海漁業(yè)資源底拖網(wǎng)調(diào)查站位

調(diào)查發(fā)現(xiàn)，本研究海域長蛇鯔由于繁殖和洄游[1]等的影響，分布具有明顯的季節(jié)性變化，在春、夏和冬季航次捕獲率很低(見表1)。長蛇鯔的繁殖期為5～7月[1]，春、夏航次均無長蛇鯔幼體的分布，而夏季長蛇鯔雖資源量較高但呈現(xiàn)集群分布，僅在小部分站位出現(xiàn)。冬季長蛇鯔群體進(jìn)行越冬洄游，于調(diào)查海域分布較少，相對資源量和出現(xiàn)頻率均為最低。秋季長蛇鯔的相對資源量和出現(xiàn)頻率均為年內(nèi)最高，滿足模型的數(shù)據(jù)量要求，因此本文僅選用2016年10月(秋)數(shù)據(jù)進(jìn)行分析。

本研究中涉及的環(huán)境因子包括水溫、鹽度、水深、葉綠素a含量、離岸距離和底質(zhì)類型。其中，拖網(wǎng)時同步使用CTD儀(CTD75M/1167)測量各個站位的水溫、鹽度、水深、葉綠素a含量等環(huán)境因子。離岸距離根據(jù)站位經(jīng)緯度計算得來，底質(zhì)類型包括砂、砂-粉砂-黏土、淤泥質(zhì)粉砂、粉砂質(zhì)黏土四類。

本研究同時考慮生物因子(餌料生物豐度)[12]對長蛇鯔分布的影響。生物因子的選擇基于已發(fā)表的長蛇鯔攝食生態(tài)研究[2]，包括戴氏赤蝦(Metapenaeopsisdalei)、短鰭魚銜(Callionymussagitta)、六絲鈍尾蝦虎魚(Amblychaeturichthyshexanema)、槍烏賊(Loligospp.)、鳀(Engraulisjaponicus)和細(xì)條天竺鯛(Apogonlineatus)。

表1 山東南部近海長蛇鯔的相對資源量

注：平均值±標(biāo)準(zhǔn)誤差。Note: Mean ± standard error.

1.2 數(shù)據(jù)處理

將各站位長蛇鯔及餌料生物的實際漁獲量標(biāo)準(zhǔn)化為拖速3.0 kn與拖網(wǎng)時間1.0 h的相對資源量(g/h)，作為資源量指標(biāo)。根據(jù)本次調(diào)查的生物學(xué)測量數(shù)據(jù)，以全長197.5 mm作為50%性成熟體長，并根據(jù)50%性成熟體長將長蛇鯔劃分為幼體和成體[13]。

在構(gòu)建模型前對環(huán)境因子進(jìn)行Pearson相關(guān)性分析，以剔除相關(guān)性較高的因子對模型擬合的影響。結(jié)果顯示，表溫和底溫，表鹽和底鹽，水深、底溫和離岸距離間存在顯著的相關(guān)性，可能對模型擬合造成影響。根據(jù)長蛇鯔底棲的生態(tài)習(xí)性，初步篩選出5個因子：水深、底層鹽度、底層水溫、表層葉綠素A含量、底質(zhì)類型。

根據(jù)長蛇鯔與餌料生物的空間重疊指數(shù)(Schoener indices, SIs)[14]篩選生物因子，公式如下：

(1)

式中：Zi(x)和Zi(y)分別表示x、y物種在i站位的相對資源量；n為站位個數(shù)。結(jié)果如表2所示，僅有槍烏賊的重疊指數(shù)較高，為0.41，其它餌料生物的重疊指數(shù)范圍為0.01～0.18。本文剔除重疊指數(shù)小于0.1的物種，篩選出槍烏賊、細(xì)條天竺鯛、短鰭和六絲鈍尾蝦虎魚作為生物因子，將其相對資源量進(jìn)行對數(shù)化處理作為豐度指標(biāo)加入到模型中。

表2 餌料生物與長蛇鯔的空間重疊指數(shù)

1.3 GAM模型

以長蛇鯔相對資源量作為響應(yīng)變量，利用GAM模型分析響應(yīng)變量與解釋變量(環(huán)境因子和生物因子)之間的關(guān)系。GAM模型的表達(dá)式如下[6]：

(2)

式中：g(μ)為連接因變量的連接函數(shù)；xi為解釋變量；α表示模型中的截距；ε為隨機(jī)誤差項，fi(xi)為解釋變量的非參數(shù)平滑函數(shù)，通過樣條平滑得到。

采用向前逐步篩選的方法構(gòu)建模型，首先將初步篩選的9個因子逐個分別帶入模型，根據(jù)AIC(Akaike information criterion)準(zhǔn)則[15]，選擇AIC最小的單因子模型，逐漸加入其他因子，篩選出AIC值最小的雙因子模型，循序漸進(jìn)直至得到一個最小AIC的模型，即為擬合效果最好的模型[16]。

AIC的計算方法如下[15]：

(3)

模型篩選過程中，考慮到水深與底溫之間顯著的相關(guān)性，對底溫進(jìn)行殘差處理[17]，以排除解釋變量間的共線性影響：

res=residual[lm(sbt～depth)] 。

(4)

式中：sbt代表底溫；depth代表水深；lm表示線性回歸。在底溫和水深同時存在的模型中，以res底溫殘差值代替底溫進(jìn)行模型的擬合。

本研究對成體、幼體的分布分別進(jìn)行了建模和檢驗，其方法是首先擬合了分類模型M1，分別模擬成體和幼體相對資源量與影響因子的關(guān)系。其中平滑函數(shù)擬合了對成體、幼體不同的回歸曲線，即與影響因子的不同關(guān)系。此外以相同數(shù)據(jù)集擬合集合模型M2，在同一站位中成體、幼體相對資源量相加進(jìn)行模型擬合，即反映了總資源量與影響因子的共同關(guān)系。

利用ANOVA方差分析，對M1、M2兩個嵌套模型的差異進(jìn)行檢驗，根據(jù)卡方檢驗結(jié)果驗證兩者間是否存在顯著差異，以檢驗成體和幼體與影響因子關(guān)系的差異顯著性。

1.4 模型檢驗

利用交叉驗證(Cross validation)檢驗所得模型的預(yù)測效果，選取80%作為訓(xùn)練數(shù)據(jù)構(gòu)建模型，剩余的20%作為驗證數(shù)據(jù)，與模型的預(yù)測結(jié)果相對照。交叉驗證過程運算1 000次，通過計算預(yù)測值與實際值間的均方根誤差(Root Mean Squared Error，RMSE)、線性關(guān)系截距斜率和決定系數(shù)R2等4個參數(shù)判斷模型的預(yù)測效果。RMSE用來衡量模擬值和真值之間的偏差，決定系數(shù)R2表示可由預(yù)測變量變異解釋的預(yù)測值線性變異的比例，截距描述模型的偏差，表示觀測與模擬之間系統(tǒng)恒定的差異，斜率描述模型的一致性，表示觀測值和實際值之間的差異，解釋為低估或高估[18-19]。

1.5 分布預(yù)測

本研究使用兩階段的GAM模型[12]進(jìn)行長蛇鯔資源分布的預(yù)測，一階段僅根據(jù)環(huán)境因子構(gòu)建餌料生物(Xj)的分布模型：

(5)

將一階段的結(jié)果，即餌料生物豐度作為變量加入二階段GAM模型，預(yù)測長蛇鯔的相對資源量(Y)分布：

(6)

以上兩式中：xj為環(huán)境因子；Xj為餌料生物j的豐度；Y為長蛇鯔的相對資源量。

用于預(yù)測的環(huán)境數(shù)據(jù)提取自FVCOM(Finite-Volume Community Ocean Model)模型，F(xiàn)VCOM是一種無結(jié)構(gòu)網(wǎng)格、有限體積、自由表面三維原始方程的海洋環(huán)流模型[20]。本文根據(jù)山東南部海域的FVCOM模型，提取33 825個數(shù)據(jù)點的環(huán)境數(shù)據(jù)(水溫、鹽度、水深)預(yù)測長蛇鯔的資源分布，并繪制相對資源量分布圖。

利用Surfer13軟件繪制長蛇鯔相對資源量分布圖，利用R(version 3.4.4)中的mgcv包[21]完成模型的擬合和檢驗過程。

2 結(jié)果與分析

2.1 長蛇鯔資源分布特征

長蛇鯔相對資源量的分布在成體、幼體間存在明顯差異。秋季幼體的分布范圍為119°30′E～122°15′E、35°N～36°30′N，主要分布在30 m等深線附近及以淺水域；成體分布范圍為119°30′E～122°30′E、35°N～36°45′N，較幼體分散，分布水深范圍更廣，大致存在海州灣水域及調(diào)查水域中部兩個分布中心(見圖2)。

(左：2016年10月成體；右：2016年10月幼體。 Left: Adult in October 2016； Right: Juvenile in October 2016.)

2.2 GAM模型

根據(jù)AIC準(zhǔn)則篩選出的最優(yōu)模型為：

log(Y+1)=α+s(log(lsp+1))+s(res)+s(depth)+s(sbs)+ε。

(7)

式中：Y為長蛇鯔的相對資源量，對Y進(jìn)行對數(shù)化處理得log(lsp+1)作為響應(yīng)變量；lsp為槍烏賊相對資源量；res為殘差處理后的底層水溫；depth為水深；sbs為底層鹽度。區(qū)分和不區(qū)分成體與幼體的兩模型M1和M2對應(yīng)的AIC值分別為523.17、545.07，累計偏差解釋率分別為78.3%、69.3%(見表3)。方差分析結(jié)果表明，M1、M2兩模型差異極顯著(P<0.01)，即成體、幼體與影響因子的關(guān)系有極顯著差異。

方差分析表明，對長蛇鯔成體相對資源量有顯著影響的因子包括槍烏賊豐度、底層水溫和水深，對長蛇鯔幼體相對資源量有顯著影響的因子包括烏賊豐度、底層鹽度和水深，對長蛇鯔整體資源量影響顯著的因子為底層水溫和水深和底層鹽度(見表3)。

表3 GAM擬合結(jié)果的方差分析

Note: lsp = relative abundance ofLoligospp., sbt = sea bottom temperature, sbs = sea bottom salinity.

GAM模型表明(見圖3)，長蛇鯔成體的相對資源量隨槍烏賊豐度和底層水溫的增加均呈現(xiàn)先上升后下降的趨勢，幼體相對資源量隨槍烏賊豐度和底層水溫的增加均呈現(xiàn)上升趨勢；成體、幼體的相對資源量隨水深增加的變化有所差異，成體在25～35 m的區(qū)間內(nèi)相對資源量急劇下降，隨后緩慢上升進(jìn)而下降，在45～50 m水域附近出現(xiàn)小的峰值，幼體相對資源量隨水深的增加始終保持穩(wěn)定下降的趨勢；成體、幼體的相對資源量隨鹽度的變化趨勢明顯不同，成體的變化趨勢不明顯，幼體的相對資源量隨鹽度的增加呈現(xiàn)穩(wěn)定上升的趨勢。長蛇鯔整體的相對資源量隨槍烏賊豐度和底層水溫的增加均呈現(xiàn)先上升后下降的趨勢；隨著水深的增加，相對資源量呈現(xiàn)總體下降的趨勢，其中，在35～50 m出現(xiàn)小的峰值；相對資源量隨鹽度的增加呈現(xiàn)上升趨勢。

圖3 M1和M2中解釋變量對長蛇鯔相對資源量的影響

2.3 長蛇鯔相對資源量的預(yù)測

交叉驗證的結(jié)果表明(見表4)，M2模型有較低的均方根誤差和較高的決定系數(shù)，兩模型線性回歸平均的截距和斜率差異不大，M2更近于無偏估計，且預(yù)測值的波動范圍更小(見圖4)。因此，不區(qū)分成體幼體的M2模型具有更好的預(yù)測效果。

表4 模型交叉驗證結(jié)果

(虛線表示1 000次線性回歸的平均值，實線代表模擬值與真值的無偏估計。The dotted lines represent the mean of 1000 times linear regressions. The solid lines represent the unbiased prediction between simulated values and true values.)

圖4 模型1 000次交叉驗證的散點圖和擬合曲線

Fig.4 Scatter plot and fitting curves for results of 1 000 times cross validation

基于預(yù)測的槍烏賊分布和提取自FVCOM的環(huán)境因子數(shù)據(jù)，本文根據(jù)已構(gòu)建的M1模型，繪制了山東南部近海秋季長蛇鯔成體和幼體的資源分布圖(見圖5)。長蛇鯔的分布模式與實際觀測接近，成體存在海州灣和調(diào)查海域中部兩個分布中心，在深水區(qū)的分布較少；幼體的分布呈現(xiàn)明顯的隨水深分層的現(xiàn)象，主要分布于30 m等深線以淺的區(qū)域，近岸相對資源量高，在深水區(qū)無分布。

圖5 山東南部近海秋季長蛇鯔成體和幼體的相對資源量分布圖

3 討論

許多魚類的棲息分布隨著生活史不同階段發(fā)生變化[22-23]，而本文作為一個研究案例，揭示了長蛇鯔成體和幼體不同的分布模式。相對資源量分布圖中長蛇鯔成體主要分布于60 m以淺海域，在123°E以東的海域幾乎無分布(見圖5)?？赡苡捎谶\動能力發(fā)達(dá)，長蛇鯔成體有著廣泛的分布，在海州灣及整個沿岸帶形成了第一個高值區(qū)，在調(diào)查海域中部形成了第二個高值區(qū)。幼體分布范圍相對狹小，呈現(xiàn)近岸高遠(yuǎn)岸低的主要趨勢(見圖5)，這也表明近岸海域是長蛇鯔的重要的產(chǎn)卵場和育幼場，這與以往的研究結(jié)果相一致[24]。

本文分析了生活史差異對模型擬合的影響，兩模型存在極顯著差異(P<0.01)，表明成體與幼體的分布隨影響因子的變化趨勢存在極顯著差異。不區(qū)分成體、幼體的模型M2會忽略成體、幼體間存在的部分差異(見圖3)。因此，在模型擬合的過程中，需要考慮個體不同生長階段的差異，以增加模型的準(zhǔn)確性。本研究也發(fā)現(xiàn)，相比于M2模型來說，考慮到生活史差異的M1模型有較大的預(yù)測誤差，這也說明隨著模型復(fù)雜度的增加，模型的預(yù)測效果可能下降，預(yù)測誤差可能增加[25]。因此，模型構(gòu)建應(yīng)考慮到擬合效果和預(yù)測效果間的權(quán)衡問題，避免過度擬合的出現(xiàn)。

長蛇鯔的相對資源量與生物因子和環(huán)境因子(水溫、水深、鹽度)顯著相關(guān)。其中餌料生物通常表現(xiàn)為正效應(yīng)，本研究中隨槍烏賊豐度的增加，長蛇鯔成體和幼體的相對資源量均呈現(xiàn)上升的趨勢，但在槍烏賊豐度達(dá)到e8(g/h)后，成體的相對資源量出現(xiàn)下降的趨勢，這可能說明在較高餌料密度條件下其它因素限制了長蛇鯔成體的資源量。水溫是影響魚類生長、發(fā)育和分布的主要環(huán)境因子之一，有研究表明長蛇鯔的適溫區(qū)間為20.9～28.4 ℃[26]。本研究中，海域底溫范圍為7.1～22.3 ℃，成體相對資源量隨水溫增加呈先上升后下降趨勢，推測成體的最適水溫低于22.3 ℃，幼體始終為上升趨勢，推測幼體的最適水溫高于22.3 ℃，即成體幼體的最適水溫存在差異。Xue等的研究表明，底溫在>21.8 ℃的區(qū)間內(nèi)開始引起長蛇鯔相對資源量下降[12]，本研究結(jié)果與其結(jié)論相近。水深對長蛇鯔相對資源量表現(xiàn)為負(fù)效應(yīng)，隨水深的增加，成體與幼體相對資源量均表現(xiàn)為下降趨勢。研究發(fā)現(xiàn)，在55 m以深水域無長蛇鯔的分布，這可能是由于調(diào)查季節(jié)受到南黃海東側(cè)冷水團(tuán)中心的影響，底層水溫低于10 ℃[27]，不適合長蛇鯔棲息分布。由于本研究中水深與底層水溫有較強(qiáng)的負(fù)相關(guān)性，這種分布模式可能是水溫與水深共同作用的結(jié)果，對成體和幼體的分布均表現(xiàn)出限制作用。另外長蛇鯔相對資源量也受到鹽度的影響，鹽度的變化對成體相對資源量無顯著影響，但對幼體影響極顯著，隨底層鹽度的上升，幼體相對資源量有明顯增加，這與Xue等的研究結(jié)果相似[12]。綜上所述，本研究認(rèn)為成體和幼體間的主要差異表現(xiàn)在對水溫和鹽度的適應(yīng)性，幼體可能對鹽度的變化更敏感，而成體受水溫的影響可能較大。

與傳統(tǒng)回歸模型相比，GAM模型可以綜合多個變量，定量描述資源與影響因素間的關(guān)系，評估各變量的重要程度，并可對資源密度等進(jìn)行預(yù)測[28-29]，本研究中建模過程可以為相關(guān)研究提供有益經(jīng)驗。例如模型的回歸方程存在多個解釋變量時，變量間可能具有不同程度的共線性，強(qiáng)共線性會使模型的擬合中出現(xiàn)方差膨脹而增大誤差，因此在模型構(gòu)建前應(yīng)加以排除[30-31]。本研究數(shù)據(jù)中，底層水溫和水深間具有較強(qiáng)的相關(guān)性，考慮到水深是一個綜合因子，為保留更多原始數(shù)據(jù)信息，本文將水深也加入分析，在模型擬合過程中，對底層水溫進(jìn)行殘差處理，以res值代替水溫進(jìn)行擬合，去除相關(guān)性的影響。這種處理得到的結(jié)果可以理解為水深保持一定時，相對資源量隨底層水溫的相對變化趨勢，但res不代表實際的數(shù)值，分析長蛇鯔的最適水溫區(qū)間等信息存在一定困難。

GAM模型研究要求高準(zhǔn)確度、可信度以及全面的調(diào)查數(shù)據(jù)，以盡可能包含影響物種分布的所有因素。本文考慮的環(huán)境因子包括水溫、鹽度、水深、葉綠素a含量、底質(zhì)類型，但未能包含所有重要環(huán)境因子，如溶解氧等。此外，長蛇鯔作為一種底層魚類，底質(zhì)類型可能是影響其分布的原因之一，但在模型的擬合過程中，底質(zhì)類型并沒有表現(xiàn)出對相對資源量顯著地影響，在最優(yōu)模型中予以排除。這可能是由于調(diào)查海域的底質(zhì)類型均一性較強(qiáng)[32]，對長蛇鯔的棲息分布影響不大。此外本文選定的底質(zhì)類型分類方式較為粗糙，只根據(jù)物理特性將底質(zhì)分為砂質(zhì)、黏土類和礫石類，而忽略了其生物因素(如植被覆蓋率等)。另一方面，調(diào)查海域是長蛇鯔秋季的重要棲息地，在特定時期長蛇鯔相對資源量大、分布廣泛，是魚類群落中重要的優(yōu)勢種[12]，確有必要進(jìn)行物種分布研究。因此在今后的研究中，需進(jìn)一步優(yōu)化研究的時間、空間尺度，覆蓋長蛇鯔的分布區(qū)域，在更大的空間尺度上研究長蛇鯔不同季節(jié)的分布模式，并獲取更全面的環(huán)境數(shù)據(jù)，增加模型的準(zhǔn)確性和預(yù)測效果。