利用幾何形態測量學方法分析外來近似物種形態差異與親疏關系

粟海軍

(1. 貴州大學 林學院, 貴陽 550025; 2.貴州大學 生物多樣性與自然保護研究中心, 貴陽 550025)

精確鑒別外來物種的地理分隔種群或姊妹種的分類學關系有利于推測判定其來源并為生物入侵防控管理服務[1]，物種親疏關系與分類學地位的鑒別，在實踐上主要利用傳統形態學進行判別，當形態特征不易分清時再輔以分子分析。但是，一些入侵多地形成許多地理隔離種群的外來物種，傳統形態學方法常難以鑒別[2]。某種麗松鼠屬松鼠(Callosciurusspp.)近幾十年來侵入歐洲，在法國、比利時和意大利對本地種和生態系統造成了嚴重危害[3]。這3個國家的外來松鼠種群在外形及毛色上極其相似，難以區分。幾何形態測量學(Geometric morphometrics，GMM)近年來已成為動植物學家研究物種形態細微變化及其他生物學問題的有力工具[4]，基于GMM方法可有效檢視形態上的細小差異并進行圖形化分析[5]。

嚙齒類下頜骨是一種在發育上具功能集成性的復合骨結構，其組件的生長與形態反應源自不同的胚胎起源和控制因素，如遺傳基因和環境壓力[6]，由于下頜骨在取食功能中的重要作用，在形態變化上較之其他骨結構更容易受到環境(食物)壓力影響[7]。因此，下頜骨最后的形態變化(大小與形狀)是發育過程中基因與環境影響交互作用的表現[8]。為此，本研究應用GMM，針對三地種群個體下頜骨的形態變化進行研究，以期通過對三者間細微的形態變化差異與變化趨勢分析，推斷三者間親疏關系。

圖1 下頜地標點設置圖示(示左下頜外側)

1 材料與方法

1.1 數據獲取與地標點設置

1.2 Procrustes擬合與形變來源分析

利用tpsSmall軟件采用廣義普魯克斯特分析(Generalized procrustes analysis，GPA)進行疊置(superimposition)擬合[10]，得到單純的形狀坐標點數據(Procrustes coordinates)與質心大小數據(Centroid size，CS)進行后續分析[11]。應用ANOVA(Procrustes ANOVA)方法將各個體形狀配置與統計形狀(consensus)的方差偏離按不同的主要影響因素進行分解，包括：個體(individual)為主因的形變(Ind)、左右對稱體各面(Side)為主因的形變(DA)、個體與左右邊交互影響(Individual×side)為主因的形變(FA)[12]；同時由于每個標本進行了兩次拍攝和數字化，大量的重復冗余產生的方差偏離也在其中，因此用于評估測量誤差(Measurement errors)[13]。此外，對于種群和性別因素影響的差異，以種群和性別作為影響因子，以形狀中各地標點為因變量(DVs)，采用MANOVA進行分析[14]，為滿足MANOVA分析前提，形狀變量的正態性通過Kolmogorov-Smirnov分析進行檢驗。

表1 地標點定義

在本研究的幾何形態測量運算中，16個地標點在二維坐標中總計32個自由度，在疊置過程中因縮放、兩坐標軸投射及旋轉逼近導致4個自由度喪失，因此形狀變量總自由度為28。

1.3 異速生長影響、形變和種群親疏關系判別

異速生長(Allometry)的影響緣自發育過程，因此無法在疊置時消除掉，因而需要在組間個體比較時去除異速生長造成的影響，即尺度校正(Size correction)[15]。利用shape與size數據進行多元回歸來判定異速生長是否存在[16]。使用完全的多因素協方差分析(MANCOVA)，以三地種群作為分組因素，log(CS)作為協變量，來檢測種群間回歸斜率的差異性[17]。當斜率相同時，利用合并組內回歸來控制異速生長的影響。利用主成分分析(Principal components analysis，PCA)在兩個水平上分析各被分析種群平均形狀與目標配置形狀間的相對形變[18]。利用1000次迭代的留一法交叉驗證判別函數分析(Leave-one-out cross-validated discriminant function analysis，DFA)來判別群間親疏關系[17]。在PCA分析中，依據通過扭曲外廓圖形(Warped outline drawings)產生的相應軸上的極端值來圖形化形變效果[19]，而在DFA中，扭曲外廓圖形展示的是被比較種群平均形狀間的形變差異。

Procrustes擬合、Procrustes ANOVA, PCA, DFA和相關的圖形化均使用MorphoJ軟件(Verson.1.06c)[20]進行；Kolmogorov-Smirnov 檢驗、MANOVA 和MANCOVA均利用PASW(SPSS v18)軟件進行。用于檢驗共有斜率的MANCOVA則利用TpsRegr(version 1.40)[21]軟件進行分析。

2 結果與分析

2.1 形變來源

如表2所示，Procrustes ANOVA表明，樣本個體間(Ind)形變以極高的統計顯著性(P=0)解釋了包含大小和形狀的總形變的絕大部分(分別為99.2%和83.8%)；個體與左右邊交互影響(Individual×side)解釋總形變的0.5%(大小)和9.5%(形狀)；左右面(Side)為主因的形變(DA)解釋貢獻最小，而測量誤差的總形變貢獻僅為0.2%(CS)和4.3%(shape)，樣本個體均方超過測量誤差1283倍(CS)和39倍(shape)，因此，圖像拍攝與數字化誤差較小可忽略不計。

表2 三地外來松鼠種群下頜質心大小與形狀的形變模式

SS: 普氏平方和Procrustes sum of squares;MS: 普氏均方Procrustes mean square;F: Goodall′sF統計值Goodall′sFstatistic;P: Goodall′sFp-values

Kolmogorov-Smirnov檢驗表明，除對稱組件中16號地標點外(D=0.131,df=109,P<0.05)，其余所有變量均符合正態分布。MANOVA則顯示在兩水平上(即Ind和FA)各種群間均存在顯著差異，而種群間和種群內的樣本性別因素影響均無顯著性差異(表3)。

表3 種群和性別間對下頜骨對稱和非對稱組件差異檢驗

2.2 異速生長與形變

MANCOVA斜率檢驗表明3種群間異速生長模式是一致的(Wilks′ Lambda = 0.499,F(56,152)= 1.126;P= 0.283)。隨后組內回歸表明，異速生長效應解釋了形變總方差的2.99%(P= 0.0015)，根據圖形化結果(圖2)，可推斷當個體(大小)發育增大時，3個種群的個體都會呈現更長的冠狀突和角突、更前移的咬股嵴前緣和水平枝背部位置輕微向內收縮。

不管是在Ind還是在FA的形變上，異速生長都表現出顯著性影響(Ind: 2.99% predicted,P= 0.0009; FA: 3.72% predicted,P<0.0001)。去除異速生長影響進行尺度校正后，對于個體水平形變，前4個主成分解釋了71.4%的總方差(PC1=39.2%,PC2=17.1%,PC3=8.2%,PC4=6.9%)；對于FA形變，前4個主成分解釋了51.3%的總方差(PC1=24.1%,PC2=11.1%,PC3=8.4%,PC4=7.7%)。以主成分1和2的圖形化結果描述(圖3)，可知在個體水平上，冠狀突和角突部分對總形變貢獻最大，其次是臼齒槽、咬肌嵴位置和頦孔；在FA上，形變貢獻并不大，且主要發生在角突部分。從PCA散點圖看，以個體水平(Ind)差異可以清晰地將樣本歸為法國、比利時和意大利3組，而以FA上的變化則不能(圖3)。

圖為形狀數據對質心大小數據的總體回歸，下頜外廓示X軸兩端由開始形狀(淡色線)到目標形狀(深色線)的異速生長總形變模式(4倍夸大扭曲外廓)

A group-centered regression of shape (Procrustes coordinates) on Log(Centroid Size); a 4-fold exaggerated warped outline drawing shows the overall allometric pattern from starting shape (light line) to target shape (dark line) at the extremes of the X-axis

圖2包含異速生長影響的尺度依賴性形變

Figure 2 Size-dependent shape variation

下頜外廓示PC軸兩端2倍夸大扭曲的相對形變，淺色線示開始形狀，深色線示目標形狀

At the extremes of the PC-axes, the relative variations of the mandibles are shown as warped outline drawings with 2-fold amplification. Light lines indicate starting shapes and dark lines represent target shapes

圖3尺度校正后的下頜形變主成分分析

Figure 3 Size-corrected principal component analysis (PCA)

2.3 種群親疏關系判別

DFA結果(表4)表明3個種群間存在顯著差異，BEL與ITA種群之間較之其他對比種群(BEL-FRA; ITA-FRA)更接近(圖3和圖4)。從圖形化結果比較來看，法國種群樣本下頜有更粗壯的水平分枝(horizontal ramus)且上升枝(ascending ramus)部分的冠狀突和角突部分略有點向內收；比利時種群樣本的下頜冠狀突和角突則稍纖瘦且往腹背方向更為舒展；而意大利種群樣本的下頜形變特征則介于兩者間而更偏向比利時。

下頜外廓示各比較種群間平均形狀相對差異，扭曲形變2倍夸大

The relative differences between the mean shapes of each population are displayed by warped outline drawings corresponding to the color of the legend, exaggerated 2-fold in scale factor for better visibility

圖4 不同外來種群間的判別函數分析

3 討論

3.1 外來種麗松鼠的形變來源

定向不對稱(DA)、波動不對稱(FA)和反對稱是生物形態中的非對稱類型[13]，前兩者能夠被Procrustes ANOVA檢測出來。DA在GMM中指的是下頜骨兩面的系統結構差異性[13]，盡管DA具統計顯著性，但由于它對總形變貢獻非常小，因而常常被忽略[22]。本研究中同樣如此，DA僅解釋總形變的0.1%(CS)和2.4%(形狀)。FA是僅次于個體間差異(Ind)的對總形變貢獻最大的因素(貢獻總形變的9.5%)，這也與其他嚙齒類的研究相似，如黃頸鼠(Apodemusflavicolli)達10.1%[23]，狐松鼠(Sciurusniger)達7.9%[22]。另外，樣本年齡可能也會對總形變有一定影響，但本研究均采用的成體樣本，且個體發育的異速生長的影響已考慮進來，因此成體后的年齡差異影響應微乎其微[12]。

3.2 異速生長與形變

異速生長效應對總形變的貢獻雖小(2.99%)，卻具有統計極顯著性，因此不能夠在總形變分析中被忽略，因為異速生長效應都會模糊掉形變差異的來源[23]，因此，進行尺度校正是十分必要的。去除異速生長效應后，總體形變模式能夠通過PCA分析反映出來。Atchley等曾指出，下頜骨表形特征發育過程極大地受到外部環境影響[24]。本研究中，3處外來種群個體間形變主要發生在包含冠狀突和角突的上升枝，而FA層次的形變主要發生在角突部位，角突部位是翼內肌(medial pterygoid)和部分咬肌(masseter)著生部位，這一區域具發育不穩定性，對環境壓力是十分敏感的，因其參與構成的肌肉咬合功能極大地受到環境因子(如食譜條件、取食壓力等)影響[25]。但是，從PCA結果來看，在個體間(Ind)層次，70%以上的總形變被前4個主成分解釋，PC1和PC2即能清晰地將種群樣本分離歸組(圖3)，而FA層次卻不能有效分組，由于FA形變主要來源于外部環境壓力[26]，因此，可推斷3個種群個體雖然受不同的環境壓力(不同的生境、食物條件、氣候等)影響產生了形狀變化，但變化尚不十分明顯。

3.3 種群親疏關系判別

由于至今這種麗松鼠(Callosciurussp.)從何時何地引入歐洲三國的仍未明確[3]，因此確定其分類地位和親疏關系具有重要意義。Mazzamuto等利用Barcoding及頭骨傳統形態測量比較了這3種種群樣本[27]，分子結果表明，法國種群與赤腹松鼠(Callosciuruserythraeus)高度相似，而意大利和比利時種群則具有不同于法國種群的相同的單倍型，但兩者頭骨形態測量特征卻又與已知的赤腹松鼠樣本相似，因此其研究結果仍未得到明確結論。本文進一步拓展了Mazzamuto等的形態學判別研究。DFA結果(圖4)表明，盡管意大利與比利時種群在馬氏距離上存在重疊，但兩者仍具有顯著差異，而法國種群則與這兩個種群明顯不同，應屬不同的分類學地位，這支撐了Mazzamuto等的分子結果[27]。因此，我們推斷意大利與比利時種群可能來源于同一外來物種或同一來源，但兩者已產生部分分化，至少在下頜形態上因地理隔離和生態適應而產生了分化。

利用GMM方法能夠細致的分析和比較出外來種群各自形態上的差異，從而推斷其親疏關系，GMM的準確性與精確性，Cardini等曾做過深入討論與展示[28]，GMM的優勢也使其在近年中迅速成為生物形態研究的熱門方法。受限于缺乏已知分類地位的樣本比較，本研究尚無法確定這些外來種群的分類學地位，但通過GMM手段有效克服了分子方法研究表型變化的不足，為進一步厘清入侵種群來源提供了有益參考。