趙天亞，劉逸文，王 琨，徐賓鐸，管清華，薛 瑩,4，任一平,4，張崇良??

(1.中國海洋大學(xué)水產(chǎn)學(xué)院，山東 青島 266003；2.海州灣漁業(yè)生態(tài)系統(tǒng)教育部野外科學(xué)觀測研究站，山東 青島 266003；3.諸城市密州街道農(nóng)業(yè)農(nóng)村服務(wù)中心，山東 諸城 262200; 4.青島海洋科學(xué)與技術(shù)試點(diǎn)國家實(shí)驗(yàn)室 海洋漁業(yè)科學(xué)與食物產(chǎn)出過程功能實(shí)驗(yàn)室，山東 青島 266237)

魚類的生長特征是漁業(yè)資源評估的重要內(nèi)容，是研究漁業(yè)資源群體數(shù)量變動(dòng)和合理利用漁業(yè)資源的基礎(chǔ)。傳統(tǒng)的漁業(yè)資源評估通常假設(shè)魚類生長是勻質(zhì)的，將固定的生長參數(shù)應(yīng)用于資源評估模型并給出管理建議。而文獻(xiàn)[1]指出，魚類的生長特征并不是勻質(zhì)的，而是存在時(shí)間和空間上的異質(zhì)性，這種異質(zhì)性可能會(huì)在一定程度上影響漁業(yè)管理決策的合理性，增加漁業(yè)管理的風(fēng)險(xiǎn)[2]。

短吻紅舌鰨(Cynoglossusjoyneri)為亞熱帶及暖溫帶淺海底層魚類，具有活動(dòng)范圍小、洄游距離短、圍捕率高、環(huán)境適應(yīng)能力強(qiáng)等特性[3]，是重要的漁業(yè)資源種類。近年來，由于過度捕撈、生境破壞及氣候變化等原因，中國近海漁業(yè)資源嚴(yán)重衰退[4]，很多海域的漁業(yè)資源群落結(jié)構(gòu)出現(xiàn)顯著變化[5-7]。目前，短吻紅舌鰨逐漸成為山東近海海域的優(yōu)勢種[8]，在食物網(wǎng)中的影響逐漸增大[9]，在海洋生態(tài)系統(tǒng)中的地位和作用不斷提高。近年來，對短吻紅舌鰨的研究主要涉及資源量及分布[10-11]、產(chǎn)卵場分布[3,12]和攝食[13-14]等方面，而對其生長的研究較為缺乏，更缺乏對其生長特征的時(shí)空異質(zhì)性的深入認(rèn)識(shí)。

本研究基于2016—2017年萊州灣及2013—2020年海州灣的漁業(yè)資源調(diào)查數(shù)據(jù)，使用Bootstrap重抽樣的ELEFAN方法，估算了這兩個(gè)海域的短吻紅舌鰨群體的VBGF參數(shù)，比較了其生長參數(shù)的空間異質(zhì)性，并結(jié)合2013—2020年海州灣漁業(yè)調(diào)查數(shù)據(jù)，研究了海州灣短吻紅舌鰨群體生長參數(shù)的年際變化。本研究旨在探討短吻紅舌鰨生長特征的時(shí)空異質(zhì)性，為相關(guān)研究提供生物學(xué)數(shù)據(jù)支撐，并為漁業(yè)資源評估和漁業(yè)管理提供參考。

1 材料與方法

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

本研究于2013—2020年在海州灣及鄰近海域進(jìn)行了漁業(yè)資源調(diào)查，共獲得短吻紅舌鰨樣本1 790尾。調(diào)查海域范圍為119°20′E—121°10′E，34°20′N—35°40′N，每年調(diào)查2次，分為春(4—5月)和秋(9—10月)兩個(gè)季節(jié)。站位設(shè)計(jì)采用分層隨機(jī)采樣的方法，每個(gè)航次設(shè)置18個(gè)站位，在各區(qū)域中隨機(jī)選取對應(yīng)數(shù)量的站位[19]。所用調(diào)查船為功率220 kW的單拖底拖網(wǎng)漁船，拖網(wǎng)速度為2～3 kn，每站拖網(wǎng)時(shí)間為1 h。調(diào)查網(wǎng)具的網(wǎng)目大小為17 mm，網(wǎng)口擴(kuò)張時(shí)網(wǎng)口高度約為6 m，寬度約為25 m。

為研究短吻紅舌鰨生長參數(shù)的空間異質(zhì)性，本研究于2016年10月、2017年1、5和8月進(jìn)行了萊州灣漁業(yè)資源底拖網(wǎng)調(diào)查，共獲得樣本1 874尾。調(diào)查海域經(jīng)緯度范圍為119°05′E—120°00′E，37°12′N—37°40′N。調(diào)查船為220 kW的單拖底拖網(wǎng)漁船，平均拖速為3 kn，拖網(wǎng)時(shí)間為1 h。調(diào)查網(wǎng)具為1 380目、每網(wǎng)目大小為4.67 cm的底拖網(wǎng)，網(wǎng)口擴(kuò)張時(shí)網(wǎng)口高度約為7.53 m，寬度約為15 m。樣品的采集和分析參考《海洋調(diào)查規(guī)范》(GB/T 12763.6—2007)[20]。

圖1 山東近海漁業(yè)資源底拖網(wǎng)調(diào)查區(qū)域

1.2 ELEFAN方法

本研究假設(shè)短吻紅舌鰨的生長符合von Bertalanffy生長方程(von Bertalanffy growth function，VBGF)：

Lt=Linf{1-exp[-K(t-t0)]}。

式中:Linf為生物體的極限體長，也稱漸進(jìn)體長或理論最大體長；K為生長參數(shù)，衡量體長趨向于極限體長的相對速率；t0為生物體的理論生長初始年齡，即理論上體長為0時(shí)的年齡。

本研究使用ELEFAN方法分別估計(jì)了海州灣2014—2019年和萊州灣2016—2017年短吻紅舌鰨的極限體長Linf和生長參數(shù)K。ELEFAN方法可以利用連續(xù)時(shí)間序列的體長頻率數(shù)據(jù)來估算魚類的VBGF參數(shù)，其主要流程為[15]：

(1)輸入基于時(shí)間序列的體長頻率數(shù)據(jù)，并按照一定的組距對體長數(shù)據(jù)進(jìn)行分組，繪制體長頻率分布圖，圖中的峰值對應(yīng)體長頻率密度較大的點(diǎn)。

(2)對體長數(shù)據(jù)進(jìn)行數(shù)據(jù)重構(gòu)。以一定的項(xiàng)數(shù)計(jì)算移動(dòng)平均值，高于相對應(yīng)的移動(dòng)平均值的點(diǎn)為波峰，低于相對應(yīng)的移動(dòng)平均值的點(diǎn)為波谷。得到一系列的波峰與波谷后，計(jì)算可擬合波峰和(Available sum of peaks，ASP)。

(3)輸入并調(diào)整Linf和K的值，觀察不同參數(shù)組合所描述的生長曲線與體長數(shù)據(jù)中波峰分布的擬合情況，計(jì)算曲線所穿過的波峰的相對總和。曲線每到達(dá)波峰一次加1，每到達(dá)波谷一次減1，所得的總和為解釋波峰和(Explained sum of peaks，ESP)。擬合效果以擬合優(yōu)度指數(shù)(Rn)衡量：Rn=10(ESP/ASP)/10。選取Rn最大的一組為最優(yōu)解，此時(shí)所輸入的VBGF參數(shù)為最優(yōu)參數(shù)。

本研究以10 mm組距對原始體長數(shù)據(jù)進(jìn)行分組，以移動(dòng)平均步長5進(jìn)行數(shù)據(jù)重構(gòu)，優(yōu)化過程選擇遺傳算法(Genetic algorithm，GA)，以提高結(jié)果的穩(wěn)健性[21]。由于采樣過程和樣本數(shù)據(jù)的獲取具有一定的隨機(jī)性，本研究使用結(jié)合了Bootstrap重抽樣的ELEFAN方法來估算生長參數(shù)，以評估參數(shù)估計(jì)的不確定性[22]。Bootstrap重抽樣是指從原始樣本中重新抽取新樣本，新樣本的數(shù)據(jù)量與原始樣本相同且抽取方式為放回抽樣，重復(fù)此過程1 000次，即可得到1 000個(gè)新樣本，將每個(gè)抽取的新樣本用ELEFAN方法估算一次VBGF參數(shù)，即可得到1 000組VBGF參數(shù)估算結(jié)果，從而獲得其概率分布。最后，利用Pauly經(jīng)驗(yàn)公式[23]逆向估算理論生長初始年齡t0，并繪制出生長曲線。

1.3 時(shí)空異質(zhì)性分析

對短吻紅舌鰨體長頻率進(jìn)行統(tǒng)計(jì)，分析不同海域、不同時(shí)間其體長分布范圍、優(yōu)勢體長組和體長平均值等，使用Kolmogorov-Smirnov檢驗(yàn)比較不同樣本間體長分布的差異。

本研究通過比較2016—2017年海州灣群體和萊州灣群體的生長特征，反映短吻紅舌鰨生長參數(shù)的空間變化，并通過研究海州灣短吻紅舌鰨群體在2014—2019年的生長特征來反映生長參數(shù)的時(shí)間變化。由于來自海州灣漁業(yè)調(diào)查的樣本中每年只存在2個(gè)月份的體長數(shù)據(jù)，ELEFAN方法的擬合效果較差[17]，因此本研究采用移動(dòng)平均的方法研究VBGF參數(shù)隨時(shí)間變化的變動(dòng)趨勢，以相鄰3年的調(diào)查數(shù)據(jù)作為其中間年份對應(yīng)的樣本，再對其進(jìn)行Bootstrap重抽樣和運(yùn)行ELEFAN，如海州灣群體2014年的VBGF參數(shù)是使用2013—2015年的數(shù)據(jù)估算得出。利用Z檢驗(yàn)進(jìn)行生長參數(shù)的差異顯著性檢驗(yàn)。

本研究過程中的數(shù)據(jù)分析及圖形構(gòu)建采用R 4.0.2環(huán)境下TropFishR[24]軟件包完成。

2 結(jié)果

2.1 短吻紅舌鰨的體長頻率分布

對2016—2017年海州灣短吻紅舌鰨樣本的體長頻率統(tǒng)計(jì)分析結(jié)果表明，海州灣群體的體長分布范圍為43～274 mm，體長平均值為(123.9±39.2)mm，優(yōu)勢體長組為133～177 mm，占總樣本量的38%。在2013—2019年，海州灣短吻紅舌鰨體長均值的最大值出現(xiàn)于2015年，優(yōu)勢體長的最大值出現(xiàn)于2016年；而體長均值和優(yōu)勢體長的最小值均出現(xiàn)于2017年(見表1)。體長均值與優(yōu)勢體長在2014—2019年整體呈波動(dòng)變化(見表1)。

表1 海州灣短吻紅舌鰨樣本體長分布的年際變化

萊州灣群體的體長分布范圍為40～225 mm，體長平均值為(119.4±28.9)mm，稍小于海州灣群體(見圖2)。其優(yōu)勢體長組為114～150 mm，占總樣本量的43.5%。對萊州灣和海州灣2016—2017年兩個(gè)群體體長分布情況的Kolmogorov-Smirnov檢驗(yàn)結(jié)果表明其具有顯著差異(P<0.05)。

圖2 2016—2017年海州灣和萊州灣短吻紅舌鰨樣本體長頻率分布

2.2 短吻紅舌鰨生長參數(shù)的空間異質(zhì)性

參數(shù)估算結(jié)果表明(見表2)，海州灣和萊州灣兩個(gè)海區(qū)中短吻紅舌鰨的VBGF參數(shù)分布存在一定差異。在2016—2017年，海州灣短吻紅舌鰨群體的極限體長Linf均值為(231.8±39.1)mm，生長參數(shù)K均值為(0.39±0.15)a-1；萊州灣短吻紅舌鰨群體的極限體長Linf均值為(206.3±39.7)mm，生長參數(shù)K均值為(0.44±0.16)a-1。萊州灣群體相比海州灣群體的生長速度更快，但體型較后者更小，即相比海州灣群體萊州灣群體更趨于小型化，呈現(xiàn)出一定程度的異質(zhì)性生長特征(見圖3a)。Z檢驗(yàn)結(jié)果表明，兩個(gè)群體的極限體長Linf和生長參數(shù)K均具有顯著差異(P<0.05)(見表2)。

2.3 海州灣短吻紅舌鰨群體生長參數(shù)的年際變化

生長參數(shù)估算結(jié)果表明(見表2)，海州灣短吻紅舌鰨極限體長Linf的均值范圍為(200.3～234.1)mm，其中最大值和最小值分別出現(xiàn)在2015和2019年，總體呈逐步下降趨勢。生長參數(shù)K的均值范圍為0.33～0.45 a-1，其中最大值和最小值分別出現(xiàn)在2019和2015年，與Linf相反，總體呈現(xiàn)逐步上升趨勢。對極限體長Linf和生長參數(shù)K的相關(guān)性檢驗(yàn)表明其具有顯著負(fù)相關(guān)關(guān)系(P<0.05)，可見海州灣短吻紅舌鰨整體呈體型變小和生長變快的趨勢(見圖3b)。短吻紅舌鰨的生長特征在2014—2017年變化趨勢較平穩(wěn)，而在2017—2019年變化幅度較大。其中2017—2018年生長參數(shù)的年際變化幅度最大，極限體長Linf由226.3 mm降低至209.2 mm，生長參數(shù)K由0.35 a-1上升至0.43 a-1。Z檢驗(yàn)結(jié)果表明，該群體的VBGF參數(shù)在2014—2019年際間(除2017年極限體長Linf)均具有顯著差異(P<0.05)。

表2 短吻紅舌鰨VBGF參數(shù)的異質(zhì)性

圖3 海州灣及萊州灣短吻紅舌鰨生長曲線

3 討論

魚類年齡與生長研究是漁業(yè)資源評估中的基礎(chǔ)性研究，研究結(jié)果對后續(xù)的評估與漁業(yè)管理具有重要意義，相關(guān)研究在國內(nèi)外十分普遍[25-27]。由于在數(shù)據(jù)有限的漁業(yè)中，年齡數(shù)據(jù)不易得到，相比之下體長數(shù)據(jù)更易獲得，因此基于體長的資源評估方法(Length-based method)得到廣泛應(yīng)用。ELEFAN方法作為一種基于體長的數(shù)據(jù)有限方法，可以根據(jù)目標(biāo)魚種的體長頻率分布來估算魚類的生物學(xué)參數(shù)，因此在國內(nèi)外受到廣泛應(yīng)用。本研究豐富了ELEFAN方法在漁業(yè)中的應(yīng)用，研究結(jié)果可為其他相關(guān)研究提供基礎(chǔ)數(shù)據(jù)。

海州灣群體和萊州灣群體生長參數(shù)的對比體現(xiàn)了不同海域魚類群體生長特征的空間異質(zhì)性。相較于海州灣群體萊州灣群體的體長組成整體更小，萊州灣群體具有更大的生長參數(shù)K和更小的極限體長Linf。一方面，萊州灣水深較海州灣淺[28]，屬半封閉海灣，水體交換能力差[29]，同時(shí)萊州灣沿岸工業(yè)的發(fā)展給該海域帶來嚴(yán)重的水體污染[30]。文獻(xiàn)[31]指出萊州灣的生態(tài)狀態(tài)較差，種類豐度低，生態(tài)環(huán)境受到較為嚴(yán)重的破壞，導(dǎo)致魚類棲息環(huán)境較差。另一方面，萊州灣的漁業(yè)資源受到長期的過度開發(fā)，魚類群體所受脅迫更嚴(yán)重，使得魚類群體中大個(gè)體數(shù)量減少，體長組成更偏向于小型化。且萊州灣生境異質(zhì)性程度較低[32]，導(dǎo)致其群落結(jié)構(gòu)異質(zhì)性較低，海州灣與萊州灣群落結(jié)構(gòu)的差異可能使短吻紅舌鰨的食餌組成和捕食難度有所不同，從而導(dǎo)致了兩個(gè)海域短吻紅舌鰨群體的異質(zhì)性生長特征。此外，竇碩增[33]的研究指出1982—1993年渤海灣雄性短吻紅舌鰨個(gè)體的極限體長Linf為250 mm，生長參數(shù)K為0.25 a-1；雌性個(gè)體的極限體長Linf為242 mm，生長參數(shù)K為0.29 a-1；這與本研究樣本的體長分布及生長參數(shù)估算結(jié)果相比具有明顯差異，可見近年來該海域短吻紅舌鰨資源狀況已有較大程度的變化。因此，環(huán)境因子、群落結(jié)構(gòu)及開發(fā)歷史等因素可能使海州灣和萊州灣海域短吻紅舌鰨的生長特征產(chǎn)生了較大的差異。有研究指出K值和自然死亡系數(shù)M的變化會(huì)對漁業(yè)生物的生物學(xué)參考點(diǎn)產(chǎn)生較大的波動(dòng)[34]，結(jié)合本研究結(jié)果，萊州灣群體和海州灣群體的生物學(xué)參考點(diǎn)應(yīng)存在較大差異。本研究建議應(yīng)當(dāng)基于此差異對兩海域制定不同的漁業(yè)管理策略，萊州灣海域應(yīng)設(shè)置更大的開捕體長，減少捕撈壓力，以促進(jìn)種群恢復(fù)。然而也有研究指出魚類的生長特征也會(huì)表現(xiàn)出緯度上的空間異質(zhì)性，南部的群體通常生長更快[35]，這與本研究的結(jié)果存在差異，其具體規(guī)律還有待進(jìn)一步深入探究。

ELEFAN方法首先對原始體長數(shù)據(jù)進(jìn)行重構(gòu)，使其出現(xiàn)一系列的“波峰”和“波谷”，再通過擬合生長曲線以估算生長參數(shù)，其特點(diǎn)在于還可以揭示因生長緩慢而具有體長重疊現(xiàn)象的高齡魚的真實(shí)年齡組成。本研究中海州灣短吻紅舌鰨樣本的體長均值與優(yōu)勢體長在2014—2019年呈無規(guī)律波動(dòng)變化，而ELEFAN方法的估算結(jié)果表明該群體具有小型化趨勢，可見VBGF參數(shù)相比體長分布能夠更好地反映漁業(yè)資源群體現(xiàn)狀。但值得注意的是，單次的ELEFAN方法容易受到采樣隨機(jī)性和算法隨機(jī)性的影響，使所估算的結(jié)果可能為局部最優(yōu)解而非全局最優(yōu)解。因此，本研究采用了Bootstrap的ELEFAN方法估算生長參數(shù)，相較于單次的ELEFAN，Bootstrap的ELEFAN方法可以得出生長參數(shù)估算結(jié)果的概率分布，以平均值或概率密度最高點(diǎn)為最終結(jié)果，可以在一定程度上避免干擾因素對估算的影響，使結(jié)果更具有“魯棒性”(Robustness)。此外，許多魚類在不同季節(jié)的生長速度具有一定程度的差異，即季節(jié)性生長。本研究在海州灣每年進(jìn)行兩個(gè)月份的調(diào)查，由于調(diào)查數(shù)據(jù)在時(shí)間上的局限，所使用的VBGF未考慮魚類的季節(jié)性生長現(xiàn)象，使用ELEFAN方法所測定的生長參數(shù)的不確定性尚待評估。今后可以結(jié)合耳石[43-44]等材料對魚類進(jìn)行年齡鑒定，以準(zhǔn)確推斷ELEFAN方法在漁業(yè)應(yīng)用中的不確定性。

在許多資源評估模型的應(yīng)用過程中，通常會(huì)選取相關(guān)魚種的其他研究的結(jié)果作為模型的輸入?yún)?shù)，而事實(shí)上在不同時(shí)空尺度下，魚類的生長參數(shù)是具有差異的，這種時(shí)空差異可能影響模型結(jié)果的真實(shí)性，從而為漁業(yè)資源評估與管理帶來一定的風(fēng)險(xiǎn)[2,45]。同時(shí)，對不同時(shí)空尺度的魚類生長特征的研究還可以揭示不同海域的生態(tài)環(huán)境狀況及其發(fā)展趨勢，從而為海洋環(huán)境及漁業(yè)資源保護(hù)提供依據(jù)。因此，對魚類生長特征時(shí)空異質(zhì)性的研究不僅可以掌握目標(biāo)魚種的資源動(dòng)態(tài)，更可以服務(wù)于漁業(yè)資源評估的其他方面。

4 結(jié)語與展望

本研究表明海州灣海域的短吻紅舌鰨群體呈小型化趨勢，同時(shí)海州灣群體和萊州灣群體的生長特征存在一定程度的空間異質(zhì)性，相關(guān)研究應(yīng)充分考慮該魚種生長特征的時(shí)空異質(zhì)性對漁業(yè)資源評估與管理的影響。本研究通過對短吻紅舌鰨生長特征的時(shí)空變化的研究，為該魚種的資源評估提供數(shù)據(jù)參考，并在模型參數(shù)不確定性方面為漁業(yè)管理策略的合理性提供重要保證。今后，可以通過對更多魚種的生長參數(shù)的研究，在充分考慮參數(shù)不確定性的條件下，結(jié)合海域理化條件、捕撈壓力、氣候變化、種間關(guān)系以及其他相關(guān)生物學(xué)特征等，通過構(gòu)建生態(tài)系統(tǒng)模型的方式對山東近海漁業(yè)資源群體進(jìn)行綜合分析，以更好地進(jìn)行漁業(yè)資源評估和漁業(yè)管理。