穩(wěn)定性同位素營養(yǎng)級估算的差異性與影響因素
——以鰱鳙為例

白子晗，王睿龍，王康

1.華中農(nóng)業(yè)大學(xué)生命科學(xué)技術(shù)學(xué)院， 武漢 430070； 2.中國科學(xué)院水生生物研究所/淡水生態(tài)與生物技術(shù)國家重點實驗室， 武漢 430072

穩(wěn)定性同位素技術(shù)因?qū)ι鷳B(tài)系統(tǒng)的影響小，能很好地反映生態(tài)系統(tǒng)中的能量流動和物質(zhì)循環(huán)，被廣泛應(yīng)用于記錄生態(tài)系統(tǒng)的特征和過程，如物種的飲食結(jié)構(gòu)、營養(yǎng)位置、資源獲取和分配模式及生態(tài)位特征等［1］。其中，營養(yǎng)級代表生物在食物網(wǎng)中所處的地位，可以用于分析食物網(wǎng)中食物鏈長度、能量流動方向和雜食程度等［2］，并逐漸成為了評估、檢測生態(tài)系統(tǒng)動態(tài)、生物多樣性變化的重要指標［3］。在2種常用的同位素標記（δ15N、δ13C）中，δ15N在捕食者與獵物之間一般體現(xiàn)為明顯且有規(guī)律的富集（通常為0.3%～0.4%），所以經(jīng)常用來估算消費者的營養(yǎng)級；而δ13C隨營養(yǎng)級升高變化不大，當(dāng)碳源的同位素特征不同時，可以用來評估不同營養(yǎng)來源的比例［4］。最簡單的營養(yǎng)級計算模型需要首先確定出基準生物的營養(yǎng)級λ以及它的δ15N值，再獲取目標生物的δ15N值和相鄰營養(yǎng)級之間的氮富集系數(shù)［5］。另外，考慮到消費者會從不只一種能量來源中獲取N元素，因而引入δ13C標記來估算不同營養(yǎng)來源的比例，從而構(gòu)建出雙基準的穩(wěn)定性同位素混合模型［4］。但是上述的簡單模型和混合模型的計算方法仍存在很多局限性，在模型的假設(shè)中，N的營養(yǎng)歧化因子（trophic dis?crimination factors， TDF）通常取定值0.34%，而C則被假定幾乎不發(fā)生營養(yǎng)歧化。但實際上，TDF值會受到如食物攝取率、發(fā)育階段、性別、體質(zhì)量等［6］多種因素的影響。因此，簡單和混合模型等傳統(tǒng)計算方法中對個體之間的變異的忽略可能會在估算營養(yǎng)關(guān)系時造成很大的誤判。

另外，傳統(tǒng)算法的估算結(jié)果往往會受到樣品采集質(zhì)量的影響。在試驗設(shè)計中，每個樣點的采樣個數(shù)一般為10～20。低于這個標準可能會導(dǎo)致數(shù)據(jù)的離散程度增大，結(jié)果的準確性降低，同時也無法反映出消費者營養(yǎng)級連續(xù)、線性的變化過程。除此之外，基準營養(yǎng)級在時間和空間上存在變異性，最終會通過食物鏈作用，影響生態(tài)系統(tǒng)中所有營養(yǎng)級生物的同位素特征［2］。

貝葉斯算法作為一種高效、便捷的計算方式，能在模型建立和參數(shù)后驗估計中引入個體變異性和抽樣誤差［7］，從而克服了傳統(tǒng)營養(yǎng)級計算方法的很多局限。但是對于其是否能改善低樣本數(shù)情況下傳統(tǒng)算法的計算結(jié)果，提高其準確性，彌補在試驗設(shè)計上的問題，還需要進一步的探索。

本研究以鰱（Hypophthalmichthys molitrix）、鳙（Hypophthalmichthys nobilis）2種濾食性魚類作為研究對象。鰱、鳙具有簡短且清晰的食物鏈，廣泛分布于我國各個流域，了解其不同區(qū)域的營養(yǎng)級狀況對整個生態(tài)系統(tǒng)有重要的意義。國內(nèi)外對鰱、鳙的食性、營養(yǎng)級都有一定程度的研究，目前主流觀點是鰱的食性偏向浮游植物，而鳙更偏向浮游動物。受食物來源營養(yǎng)等級的影響，鰱、鳙營養(yǎng)級一般在2～3之間波動，且鳙略高于鰱（圖1）。二者基于食物來源的豐度，改變攝食浮游生物的比例，以達到共存的狀態(tài)［8］。目前已有部分研究基于穩(wěn)定同位素技術(shù)對鰱、鳙的營養(yǎng)級進行了估算，但僅使用了上述的簡單和混合模型等傳統(tǒng)算法，尚無通過廣義貝葉斯算法對鰱、鳙營養(yǎng)級估計的研究。因此，本研究利用貝葉斯算法，對已發(fā)表文獻中鰱、鳙的營養(yǎng)級進行重新估算，擬探究傳統(tǒng)算法與貝葉斯算法結(jié)果的差異性，分析營養(yǎng)級估算的變異性與影響因素，以期為今后穩(wěn)定同位素估算營養(yǎng)級的研究提供理論參考。

圖1 鰱、鳙在浮游食物網(wǎng)中相互作用的示意圖Fig. 1 Schematic illustration of trophic interaction in a pelagic food web dominated by silver carps and bighead carps

1 材料與方法

1.1 數(shù)據(jù)搜集

2019年10月利用中國知網(wǎng)，搜索模式設(shè)為“全文”，搜索關(guān)鍵詞“鰱”“鳙”“穩(wěn)定同位素”“食物網(wǎng)”“營養(yǎng)級”，共搜集了15個樣點的鰱、鳙、基線（base?line）三類的δ13C、δ15N同位素數(shù)據(jù)及其輔助數(shù)據(jù)；基線分為初級生產(chǎn)者和初級消費者兩大類。初級生產(chǎn)者：顆粒有機物（particulate organic matter，POM）、藻類（algae）；初級消費者：河蜆（Corbicula flu?minea）、東北田螺（Viuiparus chuiYen）、銅銹環(huán)棱螺（Bellamya aeruginosa）。輔助數(shù)據(jù)包括樣本采集個數(shù)、樣點水體的總磷（TP）、總氮（TN）等環(huán)境數(shù)據(jù)，環(huán)境數(shù)據(jù)通過檢索與原文獻相同采樣時間和地點的文獻獲得。樣點主要分布在長江流域、太湖流域、錢塘江流域，包括萬州（30°48'N，108°24'E）、養(yǎng)鹿（31°5'N，108°34'E）、雙江（30°57'N，108°40'E）、高陽（31°5'N，108°41'E）、黃石（31°00'N，108°43'E）、宜昌（30°40'N，111°19'E）、荊州（29°47'N，112°55'E）、城陵磯（29°26'N，113°9'E）、鄂州（30°27'N，114°51'E）、湖口（29°45'N，116°14'E）、千島湖（29°33'N，118°57'E）、老虎潭（30°40'N，119°57'E）、太湖（31°16'N，120°7'E）、淀山湖（31°6'N，120°57'E），以及小興凱湖（45°20 'N，132°34'E）。

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

對搜集到的數(shù)據(jù)進行整理，表格數(shù)據(jù)直接錄入Excel，圖形數(shù)據(jù)經(jīng)過GetData Graph Digitizer軟件提取后錄Excel。當(dāng)基準是初級生產(chǎn)者時λ值應(yīng)為2；當(dāng)基準是初級消費者時λ值應(yīng)為1。檢查原文獻中的計算結(jié)果，對基準類型與λ值不對應(yīng)的情況進行更正。通過軟件R（version 4.0.0）中tRophicPosition包對鰱、鳙以及基線生物的同位素數(shù)據(jù)進行處理，處理后生成營養(yǎng)級，所有數(shù)據(jù)總結(jié)于表1。

表1 研究數(shù)據(jù)匯總與文獻來源Table 1 Summary of information about silver carp and bighead carp

1.3 數(shù)據(jù)整理分析

使用分位數(shù)圖（quantile-quantile plot）檢驗用于配對樣本t檢驗和相關(guān)分析的數(shù)據(jù)是否符合正態(tài)分布。采用配對樣品t檢驗對傳統(tǒng)算法和貝葉斯算法下的鰱、鳙營養(yǎng)級進行差異性分析；對貝葉斯算法下鰱、鳙營養(yǎng)級后驗估計的區(qū)間大小（95% CI）分別和魚樣本數(shù)、基準樣本數(shù)進行Pearson相關(guān)性分析，并計算相關(guān)系數(shù)；對貝葉斯算法下鰱、鳙的營養(yǎng)級后驗估計的眾數(shù)和水體環(huán)境中TP、TN數(shù)據(jù)進行Pearson相關(guān)性分析，并計算相關(guān)系數(shù)。選取一次函數(shù)、反比例函數(shù)、對數(shù)函數(shù)、指數(shù)函數(shù)及復(fù)合函數(shù)等作為候選模型，對貝葉斯算法下后驗分布的置信區(qū)間范圍與采集樣本數(shù)的相關(guān)性進行擬合，使用赤池信息準則（Akaike information criterion，AIC）篩選出最優(yōu)模型。

2 結(jié)果與分析

在所有的15個樣點中，鰱的δ15N取值范圍在0.759%～1.593%，最大值出現(xiàn)在太湖（1.593%），最小值出現(xiàn)在湖口（0.759%）；鳙的δ15N取值范圍在0.563%～1.68%，最大值出現(xiàn)在小江流域雙江（1.68%），最小值出現(xiàn)在鄂州（0.563%）。

傳統(tǒng)算法下鰱的營養(yǎng)級范圍在1.57～3.94、鳙為1.17～4.92，最小值和最大值均分別出現(xiàn)在萬州和雙江；貝葉斯算法下鰱的營養(yǎng)級范圍在1.64～4.10，最小值、最大值出現(xiàn)地點與傳統(tǒng)算法下一致，鳙的營養(yǎng)級范圍在1.30～5.03，最小值出現(xiàn)在鄂州，不同于傳統(tǒng)算法，最大值出現(xiàn)地（雙江）與傳統(tǒng)算法一致。

經(jīng)檢驗用于配對樣本t檢驗和相關(guān)分析的數(shù)據(jù)符合正態(tài)分布，傳統(tǒng)算法下與貝葉斯算法下的鰱、鳙營養(yǎng)級沒有顯著性差異（分別為P鰱=0.501 7，P鳙=0.897 6，表2）。貝葉斯算法下營養(yǎng)級的95%置信區(qū)間大小與魚采集個數(shù)具有顯著性負相關(guān)（P<0.001***，R=-0.758，表3），與基準采集個數(shù)也同樣表現(xiàn)出顯著性負相關(guān)（P<0.01**，R=-0.535，表3）。在水體環(huán)境方面，貝葉斯算法下鰱和鳙的營養(yǎng)級與TP、TN等環(huán)境指標未表現(xiàn)出顯著的相關(guān)性（表4）。

表2 貝葉斯算法和傳統(tǒng)算法結(jié)果的配對樣本t檢驗Table 2 Paired sample t test of bayesian algorithm and traditional algorithm results

表3 營養(yǎng)級范圍與樣本數(shù)的相關(guān)分析Table 3 Correlations between range of trophic positions and sampling number

表4 營養(yǎng)級與環(huán)境變量的相關(guān)分析Table 4 Correlations between trophic position and envionment variables

通過選用不同模型對貝葉斯算法下營養(yǎng)級后驗分布的置信區(qū)間范圍與采集樣本數(shù)之間的相關(guān)關(guān)系進行擬合，發(fā)現(xiàn)指數(shù)函數(shù)對觀測結(jié)果的解釋能力最強（AIC=63.676，表5）。根據(jù)模擬曲線（圖2），與樣本數(shù)=1時相比，將樣本數(shù)增加至6，可以使貝葉斯算法下營養(yǎng)級的后驗估計的區(qū)間大小收斂90%以上。

圖2 貝葉斯算法中營養(yǎng)級95%置信區(qū)間的大小隨魚樣本數(shù)變化的模擬曲線Fig. 2 Simulation curve of correlation between the size of 95% confidence interval of trophic position in Bayesian algorithm and sampling number of fish

表5 模型擬合信息Table 5 Summary of model fitting information

3 討論

3.1 貝葉斯算法與傳統(tǒng)算法的差異

在2016年武漢東湖的一項試驗中顯示，主要攝食浮游生物的鰱、鳙營養(yǎng)級的準確范圍應(yīng)該在2～3之間［32］，這也是許多研究者公認的，而部分搜集的研究中鰱、鳙營養(yǎng)級與這一結(jié)果有較大的偏差。如：本研究中淀山湖、荊州、宜昌樣點鳙的營養(yǎng)級超過3；黃石、雙江、養(yǎng)鹿樣點鰱、鳙的營養(yǎng)級超過3；太湖、湖口、萬州鰱、鳙的營養(yǎng)級甚至小于2（表1）。

傳統(tǒng)方法計算的鰱、鳙營養(yǎng)級誤差較大，而使用貝葉斯算法計算的鰱、鳙營養(yǎng)級與傳統(tǒng)方法比較前后無明顯差異，說明貝葉斯算法無法提升傳統(tǒng)計算方法的結(jié)果，或者說即使是利用貝葉斯算法也無法沖淡試驗設(shè)計對研究結(jié)果的影響。

3.2 樣品采集對營養(yǎng)級估算的影響

本研究中，貝葉斯算法的置信區(qū)間與魚類和基線的采集樣本數(shù)目均具有顯著相關(guān)性，說明采樣的強度確實會影響后續(xù)分析的結(jié)果。在本研究所涉及的15處樣點中，除老虎潭、千島湖、淀山湖外，其他樣點采樣數(shù)量總體偏少，多數(shù)為單次采樣、無周期性采樣，個別樣點采樣數(shù)量為1，如雙江、養(yǎng)鹿、黃石，這可能導(dǎo)致了在鰱、鳙營養(yǎng)級估計上的偏差。另外在少數(shù)位點，如鄂州、養(yǎng)鹿等，2種算法獲得的結(jié)果有較大差異。這可能是因為較低的樣本量使得傳統(tǒng)算法中所依賴的均值和貝葉斯算法所擬合的樣本分布之間出現(xiàn)了偏差。而曲線擬合的結(jié)果表明，維持采樣數(shù)在6個以上有助于提高貝葉斯算法的計算結(jié)果的準確性。

另外，有10處樣點使用了POM進行標定。POM作為一種高周轉(zhuǎn)率的基準，短期或者某一時間點的樣品采集分析無法反映出基準作用［33］，通常需要分時段、高頻率地長期檢測，其δ15N值也通常會體現(xiàn)出明顯的季節(jié)和空間變化。例如，Gu［34］將42個無季節(jié)性采樣的湖泊和36個進行季節(jié)性采樣的湖泊進行了對比，發(fā)現(xiàn)與進行季節(jié)性采樣的湖泊相比，無季節(jié)性采樣的湖泊中POM的δ15N存在持續(xù)消耗的現(xiàn)象，而且平均消耗量大于0.2%，甚至可以掩蓋浮游動物與POM之間的營養(yǎng)聯(lián)系。Gu等［35］還發(fā)現(xiàn)，在32個研究對象中，POM的δ13C的季節(jié)性變化幅度在1%到20%之間，大大超過每次營養(yǎng)轉(zhuǎn)移過程中同位素分餾的水平（<0.1%）。因此在采樣期間獲取能夠反映POM季節(jié)性周期變化的樣本非常重要，但是本文引用的很多研究僅采集了單次、少量的POM樣品，無法滿足其作為基準的標準；食物網(wǎng)富集指數(shù)也僅采用了定值0.34%［36］，這可能是造成其營養(yǎng)級計算誤差的一個重要原因。

3.3 水體環(huán)境對營養(yǎng)級估算的影響

貝葉斯算法下鰱、鳙營養(yǎng)級與環(huán)境因子TN、TP存在相關(guān)關(guān)系，這說明鰱、鳙所處的水體環(huán)境對其營養(yǎng)級具有一定程度的影響。Xu等［37-39］利用穩(wěn)定碳、氮同位素技術(shù)研究了我國東湖、巢湖、撫仙湖等淡水湖泊的水生食物網(wǎng)結(jié)構(gòu)與營養(yǎng)級關(guān)系之后，提出外來營養(yǎng)物質(zhì)輸入是導(dǎo)致鰱、鳙、湖鱭等浮游食物網(wǎng)中次級消費者的穩(wěn)定同位素值升高的原因，這與本研究的結(jié)果相似。此外，生態(tài)系統(tǒng)的其他組分也受到富營養(yǎng)化的影響，如Gu［34］的研究指出，由于水體營養(yǎng)水平不同，富營養(yǎng)湖泊POM的δ15N的變化幅度（0.26%～2.5%，平均值為1.03%）遠大于貧營養(yǎng)湖泊（0.15%～0.85%，平均值為0.42%），因此，水體營養(yǎng)環(huán)境對營養(yǎng)級估算的影響不可忽略。

綜上所述，穩(wěn)定同位素技術(shù)為研究生態(tài)系統(tǒng)各種群之間的能量流動、營養(yǎng)關(guān)系提供了新的技術(shù)手段，但是該方法的應(yīng)用還有賴于科學(xué)的試驗設(shè)計和對同位素組成的時間、空間、種內(nèi)的全面理解，從其結(jié)果和采樣個數(shù)來看，原文獻的數(shù)據(jù)還不能很好地代表目標生物的營養(yǎng)特征，使鰱、鳙的營養(yǎng)級計算結(jié)果出現(xiàn)了較大的偏差。為展示研究對象連續(xù)、完整、準確的營養(yǎng)特征，今后可展開多月份連續(xù)、適當(dāng)數(shù)量的采樣，同時還要考慮樣點的水體環(huán)境的影響。

采樣之前，要明確研究整體所代表的時間尺度，基于不同的時間尺度，分析系統(tǒng)中、短期和長期的結(jié)構(gòu)與功能變化也要依據(jù)不同檢測手段，進行相互補充、相互優(yōu)化［36］。在研究長期結(jié)構(gòu)時，同位素樣品在采樣期間需要進行密集采樣，每個樣點采集10～20個，每月或每2個月采樣1次；考慮到高周轉(zhuǎn)率樣品同位素的不穩(wěn)定性，中、短期樣品采集個數(shù)與采樣頻率應(yīng)適當(dāng)增加，每個樣點采集20～30個，每2周采樣1次。