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

白子晗，王睿龍，王康

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

穩(wěn)定性同位素技術因對生態(tài)系統(tǒng)的影響小，能很好地反映生態(tài)系統(tǒng)中的能量流動和物質(zhì)循環(huán)，被廣泛應用于記錄生態(tài)系統(tǒng)的特征和過程，如物種的飲食結構、營養(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)級升高變化不大，當碳源的同位素特征不同時，可以用來評估不同營養(yǎng)來源的比例［4］。最簡單的營養(yǎng)級計算模型需要首先確定出基準生物的營養(yǎng)級λ以及它的δ15N值，再獲取目標生物的δ15N值和相鄰營養(yǎng)級之間的氮富集系數(shù)［5］。另外，考慮到消費者會從不只一種能量來源中獲取N元素，因而引入δ13C標記來估算不同營養(yǎng)來源的比例，從而構建出雙基準的穩(wěn)定性同位素混合模型［4］。但是上述的簡單模型和混合模型的計算方法仍存在很多局限性，在模型的假設中，N的營養(yǎng)歧化因子（trophic dis?crimination factors， TDF）通常取定值0.34%，而C則被假定幾乎不發(fā)生營養(yǎng)歧化。但實際上，TDF值會受到如食物攝取率、發(fā)育階段、性別、體質(zhì)量等［6］多種因素的影響。因此，簡單和混合模型等傳統(tǒng)計算方法中對個體之間的變異的忽略可能會在估算營養(yǎng)關系時造成很大的誤判。

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

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

本研究以鰱（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)定同位素技術對鰱、鳙的營養(yǎng)級進行了估算，但僅使用了上述的簡單和混合模型等傳統(tǒng)算法，尚無通過廣義貝葉斯算法對鰱、鳙營養(yǎng)級估計的研究。因此，本研究利用貝葉斯算法，對已發(fā)表文獻中鰱、鳙的營養(yǎng)級進行重新估算，擬探究傳統(tǒng)算法與貝葉斯算法結果的差異性，分析營養(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)，搜索模式設為“全文”，搜索關鍵詞“鰱”“鳙”“穩(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。當基準是初級生產(chǎn)者時λ值應為2；當基準是初級消費者時λ值應為1。檢查原文獻中的計算結果，對基準類型與λ值不對應的情況進行更正。通過軟件R（version 4.0.0）中tRophicPosition包對鰱、鳙以及基線生物的同位素數(shù)據(jù)進行處理，處理后生成營養(yǎng)級，所有數(shù)據(jù)總結于表1。

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

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

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

2 結果與分析

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

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

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

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

通過選用不同模型對貝葉斯算法下營養(yǎng)級后驗分布的置信區(qū)間范圍與采集樣本數(shù)之間的相關關系進行擬合，發(fā)現(xiàn)指數(shù)函數(shù)對觀測結果的解釋能力最強（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)級的準確范圍應該在2～3之間［32］，這也是許多研究者公認的，而部分搜集的研究中鰱、鳙營養(yǎng)級與這一結果有較大的偏差。如：本研究中淀山湖、荊州、宜昌樣點鳙的營養(yǎng)級超過3；黃石、雙江、養(yǎng)鹿樣點鰱、鳙的營養(yǎng)級超過3；太湖、湖口、萬州鰱、鳙的營養(yǎng)級甚至小于2（表1）。

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

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

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

另外，有10處樣點使用了POM進行標定。POM作為一種高周轉率的基準，短期或者某一時間點的樣品采集分析無法反映出基準作用［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)轉移過程中同位素分餾的水平（<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存在相關關系，這說明鰱、鳙所處的水體環(huán)境對其營養(yǎng)級具有一定程度的影響。Xu等［37-39］利用穩(wěn)定碳、氮同位素技術研究了我國東湖、巢湖、撫仙湖等淡水湖泊的水生食物網(wǎng)結構與營養(yǎng)級關系之后，提出外來營養(yǎng)物質(zhì)輸入是導致鰱、鳙、湖鱭等浮游食物網(wǎng)中次級消費者的穩(wěn)定同位素值升高的原因，這與本研究的結果相似。此外，生態(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)定同位素技術為研究生態(tài)系統(tǒng)各種群之間的能量流動、營養(yǎng)關系提供了新的技術手段，但是該方法的應用還有賴于科學的試驗設計和對同位素組成的時間、空間、種內(nèi)的全面理解，從其結果和采樣個數(shù)來看，原文獻的數(shù)據(jù)還不能很好地代表目標生物的營養(yǎng)特征，使鰱、鳙的營養(yǎng)級計算結果出現(xiàn)了較大的偏差。為展示研究對象連續(xù)、完整、準確的營養(yǎng)特征，今后可展開多月份連續(xù)、適當數(shù)量的采樣，同時還要考慮樣點的水體環(huán)境的影響。

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