3D Geometric Morphometric Protocol - Quantifying the Morphology of Living and Extinct Vertebrates   

研究背景

自然界生物体形态上的差异可以反映出其功能、发育过程、对选择压力的反应等各个方面的区别。形态测量学 (morphometrics) 即是对于生物体的形态差异进行量化的方 法，主要包括传统形态测量 (traditional morphometrics) 和几何形态测量 (geometric morphometrics，以下简称GMM)。传统形态测量获取和分析的数据主要是对于距离的线性测量值，而几何形态测量获取和分析的对象则是一系列分布在二维或三维空间中的点的坐标 (Kendall, 1984; Adams et al., 2004; Zelditch et al., 2012)。与几何形态测量相比，传统形态测量具有以下缺陷：1)测量项间往往彼此不独立；2)缺少测量项间的空间关系信息；3)难以保证同源性；4)分析结果无法可视化因此难以进行直观的解释。相比之下，几何形态测量则能够提供对于形态更加准确、精确和全面的描述，同时结果的可视化也使其更易于被解读、展示和交流 (Zelditch et al., 2012)。
        形态与功能的关系一直是生物学者关注的重点之一。由于化石物种的生态习性无法直接观察，只能通过与现生物种进行形态-功能关系的类比，对其生态习性进行推测，因此功能形态学对于古生物学来说尤为重要。近年来，计算机断层 (CT) 扫描及表面扫描等技术的进步和推广，以及学界对于加强扫描数据开放获取的共识，使得研究人员获取大量现生及化石标本的三维数据进行分析成为可能 (Davies et al., 2017)。相比于过去的二维几何形态测量 (2D GMM)，基于三维模型的三维几何形态测量 (3D GMM)能够更加全面和准确的收集形态数据，更适宜于对多个结构进行综合分析，收集的数据也更具有重复使用性，因此近年来在包括鸟类、哺乳类和有鳞类等脊椎动物各个类群中均有所应用 (Fabre et al., 2014; Benson et al., 2017; Cooney et al., 2017; Felice and Goswami, 2018; Navalón et al., 2019; Watanabe et al., 2019)。本文将以探索现生鸟类下颌形态与功能之间的关系为例，演示使用三维几何形态测量方法对形态进行量化的基础操作流程。此实验流程同样适用于将化石物种的数据加入现生物种的数据集，用以推测化石物种的生态习性。
        注：本次分析仅用于演示三维几何形态测量方法的分析流程，未进行完整全面的样本选取，因此不能作为解释鸟类下颌形态与功能关系的最终研究结果。

实验步骤

1. 获取三维模型
   
根据实验目的，选取实验所涉及的样本。首先建议在三维模型在线共享平台上查看目标标本是否已有开放获取的数据，如Morphosource (https://www.morphosource.org/)，Phenome10k (https://www.phenome10k.org/) 和Aves 3D (https://www.aves3d.org/) 等，从而节省自行扫描的时间和经济成本。读者在下载使用开放获取数据时，须注意相关数据的使用权限，如是否可用于商业用途等。对于暂时没有开放获取数据的样本，则须自行挑选标本进行CT或表面扫描等，并视实验需要进行分割建模。
        本次实验中涉及的三维模型均选自Bjarnason and Benson (2021)公开发表并存储于Morphosource平台上的数据集：https://www.morphosource.org/projects/00000C420。原数据集包括149件现生鸟类骨骼样本的扫描及三维重建数据 (Bjarnason and Benson, 2021)，涵盖绝大部分现生鸟类支系 (Prum et al., 2015) 。本次实验选取其中的45个样本进行分析，选取时尽可能涵盖更全面的系统发育支系，但具有不同喙部功能的类群。最终数据集中包括有现生鸟类45个种，45个属，36个科，23个目。本实验根据Navalón et al. (2019)将所选取样本的喙部功能分为5类：捡拾 (Grabbing/gleaning，15种)、撕扯 (Tearing，9种)、啄食 (Pecking/grazing，10种)、探食 (Probing，5种)、咬合破开 (Cracking/biting，6种)。以其中一个样本斑尾塍鹬Limosa lapponica为例进行标志点获取的演示：标本号 NHMUK S/1994.46.20；馆藏于Natural History Museum bird collection, Tring UK；数据集中名称为“Limosa_lapponica_NHMUK_S_1994.46”；三维模型见附件1，Checkpoint文件见附件2。


2. 获取标志点及半标志点
   三维几何形态测量中的目标样本形态由标志点 (Landmark) 和半标志点 (Semilandmark) 在三维空间中的坐标进行描述。标志点为离散的同源点 (标志点 的详细定义见Zelditch et al., 2012)，半标志点则是一系列构成同源曲线/平面的连 续的点，但曲线/平面上的单个半标志点可以不是同源点 (半标志点的详细定义见 Philipp and Philipp, 2013)。标志点和由半标志点构成的曲线/平面的选取应满足以 下原则：1) 同源；2) 能够提供足够多的形态信息；3) 在每个样本上都能够找到。标志点数目和曲线/平面的数目在各样本间应保持一致。
           常见的标志点一般有三类—Type I: 可以客观依据生物学意义进行定义的 点，如结构上的孔洞、结构的交点等；Type II: 具有一定生物学意义，但由几何学 进行定义的点，如结构的顶点、曲线的拐点等；Type III: 只有几何学意义，无法独 立进行定义的点，如连线的中点等。这三种类型的标志点可靠性依次降低，因此选取时尽可能选取第I类和第II类而避免第III类 (Bookstein, 1997; 史宇坤，2017)。
           本实验选取标志点和由半标志点构成的曲线来代表样本的形态，并使用滑动 (Sliding) 的方法对半标志点进行计算。由于鸟类下颌形态相对简单，本实验未使用由半标志点构成的平面，但相关信息可参阅其他参考文献 (Philipp and Philipp, 2013; Bardua et al., 2019)，相关方法在其他研究中也已经有所应用 (Felice and Goswami, 2018; Bardua et al., 2019; Goswami et al., 2019; Watanabe et al., 2019)。为了简化数据处理流程，本实验对各条曲线选取固定数目的半标志点。
           对于两侧对称的样本如本次实验的对象——鸟类下颌，为了最大化可取样的样 本量，且避免两侧下颌在标本制取过程中位置相对扭动造成的误差，可以并建议仅对单侧进行取样 (Cardini, 2016 and 2017)。本次实验即选取右侧下颌进行单侧取样，而用于演示的数据“Limosa_lapponica_NHMUK_S_1994.46”以及其他部分模型左侧保存较好，因此在取样前进行了镜像处理 (Bjarnason and Benson, 2021)。 实验共选取了6个标志点和由半标志点构成的4条曲线，标志点及半标志点的定 义和描述见图1及以下 (修改自Bjarnason and Benson, 2021)：
   
   
   图1. 下颌标志点和半标志点示意图 在一件华丽琴鸟 (Menura novaehollandiae， 标本号FMNH 336751，馆藏于Field Museum of Natural History, Chicago, USA) 的下颌上对标志点 (LM1-LM6) 和半标志点构成的曲线 (SL1-SL4) 进行标识，上方为背视，下方为腹视。绿色标志点为曲线起点，黄色半标志点构成曲线，红色标志点为曲线终点。原图来自于Bjarnason and Benson (2021)。
   
           标志点1：LM1_Mandibular symphysis_cranialmost point：两侧下颌联合部 的背前侧顶点，通常为下颌最前端的点；同时是下颌背缘曲线的起点
   (SL1) 和下 颌联合部曲线的起点 (SL4)。
           标志点2：LM2_Mandible_contact of dorsal margin of the mandible with the glenoid：下颌背缘与后端关节窝前缘的交点；同时是下颌背缘曲线的终点(SL1)。
           标志点3：LM3_Articular glenoid_craniolateral point：后端关节窝的前侧顶点；同时是后端关节窝后缘曲线的起点(SL2)。
           标志点4：LM4_Articular glenoid_medial process：后端关节窝的内侧顶 点；同时是后端关节窝后缘曲线的终点(SL2)。
           标志点5：LM5_Mandibular symphysis_caudoventral point：两侧下颌联合部 的腹后侧顶点；同时是下颌腹缘曲线的起点(SL3)和下颌联合部曲线的终点(SL4)。
           标志点6：LM6_Mandible_caudoventral point：下颌腹缘向后最远点；同时是 下颌腹缘曲线的终点(SL3)。
           曲线1：SL1_Dorsal crest of mandible：下颌背缘曲线，起于标志点1，终于 标志点2；包括12个半标志点。
           曲线2：SL2_Articular glenoid rim：下颌后端关节窝后缘曲线，起于标志点 3，终于标志点4；包括6个半标志点。
           曲线3：SL3_Ventral surface of mandible：下颌腹缘曲线，起于标志点5，终 于标志点6；包括6个半标志点。
           曲线4：SL4_Mandibular symphysis_rostroventral surface：两侧下颌联 合部腹侧曲线，起于标志点1，沿腹侧延伸，终于标志点5；包括3个标志点。
   
           获取二维标志点的免费软件较多，而获取三维标志点的免费软件则相对较少。部分用于处理分割CT扫描数据的综合性三维模型处理软件，同时也具有选取标志点功能，如本次实验中利用的原始数据集便是由Avizo 9.3 创建。本次实验则选用专门用于形态分析的软件Stratovan Checkpoint 为例，方便习惯于不同三维模型处理软件的读者使用。Checkpoint每年须付费更新证书，但可以在官网申请15天的免费试用。
           打开Checkpoint后，会自动进入Specimen界面 (也可在左侧边栏点击进入 其他界面)。点击Browse选定存储三维模型的文件夹，软件会显示相应文件夹内可用的模型列表。选中目标模型后，点击open specimen即可打开该模型进行操作。左侧边栏点击进入Surface界面即可开始进一步操作。和其他三维模型处理软件一样，可使用鼠标和键盘对模型进行缩放、旋转和移动(图2)。
   
   
   图2. Checkpoint选取标志点时的Surface工作界面 红色箭头所指为选取标志点按钮，和标志点放大缩小按钮。
   
   
   1)

   选取标志点：标志点在Checkpoint中称为Single points, 添加方式为点击 按钮后，按住shift键的同时在目标位置单击鼠标进行添加，标志点将显示为圆球。同时可以对圆球大小进行调整以减少遮挡，方便查看。选中已有标志点后可对其进行移动，当前选中的标志点和曲线会显示为黄色。在工作视图左上角选中相应的标志点或曲线后点击右键，可对其进行重命名和删除等操作。添加全部6个标志点并进行重命名 (图2)。
   2)

   选取曲线 (半标志点)：由半标志点构成的曲线在Checkpoint中称为Curve， 添加方式为点击 按钮后，按住shift键的同时在目标位置按顺序连续单击鼠 标，添加一系列半标志点。半标志点显示为彼此相连的方块，且上方可显示曲线的方向。本次实验中的曲线起终点均在后续R程序中由上文中已添加的标志点进行定义，因此在Checkpoint工作界面上的曲线不包括与起终点的连线 (图3)。读者也可以尝试手动选择曲线的起终点以在工作界面形成完整的连线，但须注意由此带来的半标志点数量变化，并在后续R程序中对滑动点进行定义时加以考虑，同时确保手动选取的曲线起终点与相应的标志点高度重合。
   
   
   图3. Checkpoint选取曲线时的Surface工作界面 红色箭头所指为选取曲线 按钮。
   
           Checkpoint会自动在手动选取的半标志点之间增加更多的半标志点，形成与模型更加贴合的曲线。因此，在添加完半标志点后建议左侧边栏点击进入Landmark界面，修改确定每条曲线所包括的半标志点的数量，然后再根据半标志点的连线修改各个半标志点的位置，使其在目标半标志点数量下，形成最贴合的曲线。在曲线形状简单时，则可以利用此自动添加功能减少手动选取半标志点的工作量，如图4中即手动为包括三个半标志点的曲线4选取了两个半标志点，中间的第三个半标志点 (不带箭头的蓝色小方块) 则由软件自动生成。
   
   
   图4. Checkpoint自动添加半标志点 红色箭头所指为自动添加的半标志点。
   
   
   3)

   输出前检查：输出结果前应在Landmark界面中检查确定：所有样本上的标志点和半标志点均严格按照同一顺序/方向进行添加，且命名和数量在各样本间保持一致。若曲线方向取反，可选中后点击flip按钮进行翻转 (图5)。对标志点和曲线命名时避免使用逗号“,”，否则在后续数据处理中会被读取为不同字段。
   
   
   图5. Checkpoint输出结果前检查 红色箭头所指为flip按钮。
   
   4)

   输出结果：左侧点击进入Export界面，将文件格式选为CSV，点击Export Landmarks输出文件 (图6)。此时所有的标志点和半标志点信息将被存储在 一个csv文件中。为简化后续流程，我们手动将文件中包含有冗余信息的第一行 (基本信息行)、第一列 (序号列) 、第三列 (type列) 、以及最后三列NX-NZ删除，形成一个仅包含所需标志点和半标志点信息的名为“Limosa_lapponica_NHMUK_S_1994.46.20.csv”的文件(图7，csv文件见附件3)，并将其添加进已有数据集的文件夹(附件中名为“Mandible landmark csv files”的文件夹)。本文附件中的标志点数据集中已包括有此数据，读者可以直接利用其进行后续步骤，也可以尝试替换成自己获取的数据。
   
   
   图6. Checkpoint结果输出界面
   
   
   图7.Limosa_lapponica_NHMUK_S_1994.46.20样本的csv文件截图
   
   注：标志点/半标志点的.csv文件中的格式在各个样本中必须严格保持一致 (即 图7中的格式)，包括行、列的数量，以及其中标志点/半标志点的数量、名称和顺序。
           所有的半标志点和标志点信息随后将被导入开源软件R (实验中所用为 3.6.0版本) 中进行进一步的分析。
   
   
3.  普罗克鲁斯叠加
   针对获取的标志点和半标志点信息所做的第一步处理是进行普罗克鲁斯叠加分析 (Generalized Procrustes Analysis, GPA)，即对各样本的坐标点进行平移 (translate，统一位置)、旋转 (rotate，统一方向)、缩放 (scale，统一大小) 等操 作——这些操作可以消除样本上与形状 (shape) 无关的因素 (Gower, 1975; Rohlf and Slice, 1990; Zelditch et al., 2012)。Zelditch et al. (2012)书中第52页图3.1和 第76页图4.1对GPA中的这一系列操作有详细的图解。另外，在使用本实验中所采用的滑动半标志点方法时，各半标志点在所构成曲线上的滑动也在GPA中完成，从而减少主观对曲线进行定义的影响 (Philipp and Philipp, 2013)。
           GPA将输出各个样本经过叠加操作后的“普罗克鲁斯坐标” (Procrustes coordinates)，以及各个样本的原始大小信息。GMM中使用中心大小即centroid size代表体型大小信息，中心大小是样本上所有标志点到中心点距离的平方差求和后的平方根，是在几何学意义上独立于形状 (shape) 的指标 (Zelditch et al., 2012)。关于中心大小的详细图解可见Zelditch et al. (2012)中第59页图3.7。中心大小信息可在 后续统计学分析中参与运算，但本实验中选择剔除大小信息，仅对形状即shape信 息进行分析讨论。详细的几何形态学理论知识参阅文末“推荐阅读”部分，R的基础操作见其在线指南“An Introduction to R”及其他资料。本实验中使用R包“geomorph”中的gpagen()函数完成GPA (Adams and Otárola‐Castillo, 2013)。R 程序及注释见附件4，同时解释如下：
   
           #下载安装R包“geomorph” (Adams and Otárola‐Castillo, 2013)。 “geomorph”是专门用于分 析几何形态学数据的R包，本实验中使用的是version 3.1.2，但所用的函数也可能同时包含在其他版本中。
                require( geomorph )
   
           #创建函数用于在后续可视化过程中，在三维形态上对半标志点进行连线以完整显示出曲线。
               plot.semilandmark.lines <- function( landmark.matrix , curves , col = "grey" ) {
               lines.plot.temp <- rbind( curves[ , 1:2 ] , curves[ , 2:3 ] )
               for( j in 1:nrow( lines.plot.temp ) ) {lines3d( landmark.matrix[ lines.plot.temp[j,] , ] , lwd = 1 , col = col ) }
               }
   
           #将工作路径设置为存放标志点文件夹的文件夹
          setwd( "~/desktop/GMM_protocol" )
   
           #使用define.sliders()函数对各条曲线上进行滑动处理的半标志点进行定义
           sliders = rbind(define.sliders(c(1,7:18,2), write.file = F),
                                   define.sliders(c(3,19:24,4), write.file = F),
                                   define.sliders(c(5,25:30,6), write.file = F),
                                   define.sliders(c(1,31:33,5), write.file = F))
   
           #获取所有样本的名称
           specimen.names <- gsub( ".csv" , "" , list.files( "Mandible landmark csv files" , full.names =   ) )
           landmark.files <- list.files( "Mandible landmark csv files" , full.names = T )
   
           #读取所有csv文件中的标志点和半标志点信息
           landmark.list <- lapply( landmark.files , read.csv , row.names = 1 )
           names( landmark.list ) <- specimen.names
           landmark.array <- array( unlist( landmark.list ) , dim = c( nrow( landmark.list[[1]] ) ,  ncol( landmark.list[[1]] ) , length( landmark.list ) ) )
   
           dimnames( landmark.array )[[ 1 ]] <-rownames( landmark.list[[1]] )
           dimnames( landmark.array )[[ 2 ]] <-colnames( landmark.list[[1]] )
           dimnames( landmark.array )[[ 3 ]] <- names( landmark.list )
   
           #使用gpagen()函数对读取的信息进行GPA，使用上面定义的滑动点sliders作为参数
           GPA.fit <- gpagen( landmark.array, curves = sliders, ProcD = F )
   
           #绘制GPA结果中某个样本和平均形态的对比图，此处选取样本Accipiter_nisus_NHMUK_S1982.149.1作为示例。
   
           #获取GPA结果中的平均形态
                mean<-GPA.fit$consensus
   
           #选取并绘制GPA结果中的目标样本
               choose.taxon <- "Accipiter_nisus_NHMUK_S1982.149.1"
               plot3d( GPA.fit$coords[,, choose.taxon], size = 8 , col = "red", box = F, axes = F, xlab = "", ylab = "", zlab = "")
               aspect3d( "iso" )
               plot.semilandmark.lines (GPA.fit$coords[,, choose.taxon], curves = sliders, col = "red")
   
            #绘制GPA结果中的平均形态并保存
               points3d( GPA.fit$consensus , size = 8 , col = "grey" )
               plot.semilandmark.lines (GPA.fit$consensus, curves = sliders, col = "grey")
               rgl.postscript("shape
               comparison_Accipiter_nisus_NHMUK_S1982.149.1.pdf",fmt="pdf",drawText = TRUE)
   
   
4. 数据分析
   普罗克鲁斯叠加后形成的样本坐标即可进行多变量分析，如主成分分析 (Principal Component Analysis)，典型变量分析 (Canonical Variate Analysis)，聚类分析
   (Cluster Analysis) 等。分析方法的深入学习可参阅Zelditch et al. (2012), Johnson and Wichern (2001) 以及向东进等 (2005)。本实验中使用“geomorph”中的
   plotTangentSpace()函数进行主成分分析的演示——通过坐标变化，用少数变量即 主成分Principal Components (PCs) 集中表达出样本间的主要差异，便于研究者进行判断 (Zelditch et al. 2012)。本实验使用数据集中鸟类的喙部功能进行分组，分组信息见附件中名为“taxa_list_GMM protocol_selected.csv”的文件。
   
           #使用plotTangentSpace()函数对GPA的结果进行主成分分析
           PCA.fit <- plotTangentSpace (GPA.fit$coords)
   
           #以下尽量使用R中自带函数如plot()用于图表绘制，读者可选用其他R包中的绘图函数进行自行调整。
   
           #对PCA结果图中的X轴和Y轴分别进行命名和标注：X轴为PC1，Y轴分别为PC2和PC3，并对各个主成分所代表的差异占比进行标注。
           X <- "PC1"
           XLAB <- paste( X , " (" , round( PCA.fit$pc.summary$importance[ "Proportion of Variance" , X ] *100, 2) , "%)" , sep = "" )
           Y1 <- "PC2"
           Y1LAB <- paste( Y1 , " (" , round( PCA.fit$pc.summary$importance[ "Proportion of Variance" , Y1 ] *100, 2) , "%)" , sep = "" )
           Y2 <- "PC3"
           Y2LAB <- paste( Y2 , " (" , round( PCA.fit$pc.summary$importance[ "Proportion of Variance" , Y2 ]*100 , 2) , "%)" , sep = "" )
   
           #读取PCA结果和喙部功能分组信息
           write.csv(file="PCA_mandible.csv",PCA.fit$pc.scores)
           tt<-read.csv(file="PCA_mandible.csv",header=T)
   
           beakcat <- read.csv(file="taxa_list_GMM
           protocol_selected.csv", header=T)
           beak<-beakcat[,2]
   
           PCA.beak<-data.frame(tt,beak)
   
           #在同一视图中绘制PC1-PC2和PC1-PC3结果图，添加图例并分别导出未标注属名及标注属名的结果图(图8，9)
   
   
   图8. PCA结果图 显示各样本在由PC1、PC2和PC3构成的形态空间中的位置 (未标注属名) 。括号中标注的数值为各PC所代表的差异占比。
   
   
   图9. PCA结果图 显示各样本在由PC1，PC2和PC3构成的形态空间中的位置 (标注属名) 。若用于正式发表则属名应改为斜体。
   
               dev.new( width = 4.5 , height = 9 )
               split.screen( c(2,1) )
   
               screen( 1 )
               par( mar = c(4,4,1,1) )
               plot( PCA.beak$PC1, PCA.beak$PC3, bty = "l" , xlab = XLAB , ylab = Y2LAB, col = beak)
   
            #可对图中所有样本以属名进行标注，以便未来查看
           genus<-beakcat[,3]
           text(PCA.beak$PC1, PCA.beak$PC3, labels = genus, cex=0.7, pos=4, offset=0.2)
   
               screen( 2 )
               par( mar = c(4,4,1,1) )
               plot( PCA.beak$PC1, PCA.beak$PC2, bty = "l" , xlab = XLAB , ylab = Y1LAB, col = beak)
   
               legend("bottomleft", legend = c("Cracking/biting", "Grabbing/gleaning", "Pecking/grazing", "Probing", "Tearing"),
                   pt.cex = 1.5, cex = 0.8, bty = "n", pch = 1,
                   col=1:length(beak))
   
   
5. 结果可视化
   GMM的一大优势是可以对分析结果进行可视化，从而直观地显示出样本间形态的差异，通常使用薄板样条法 (Thin-plate Spline) 实现 (Bookstein, 1989; Zelditch et al., 2012)。本实验使用“geomorph”中的plotRefToTarget()函数进行对以上PCA 结果可视化的演示 (图10)，绘制出的形态为各个PC轴上最大值和最小值反向映 射出的形态。下列以绘制PC1最大值对应的下颌形态为例列出代码，读者可根据需要将代码中的PC1最大值 (PC1max) 更换为其他值如PC1最小值 (PC1min)、 PC2最大/小值 (PC2max/min)、 PC3最大/小值 (PC3max/min) 等。包含图10中 所有结果的完整代码见附件4。
   
   
   图10. 分析结果可视化演示图 A. PC1最小值对应的下颌侧面形态；B. PC1最大值对应的下颌侧面形态；C. PC2最小值对应的下颌侧面形态；D. PC2最大值对应的下颌侧面形态；E. PC3最小值对应的下颌侧面形态；F. PC3最大值对应的下颌侧面形态；G.样本Accipiter_nisus_NHMUK_S1982.149.1与平均形态的对比。红色轮廓为目标形态，灰色轮廓为平均形态。
   
           #使用plotRefToTarget()函数，以平均形态为参考，绘制出PC1最大值所对应的三维形态
         plotRefToTarget(mean,PCA.fit$pc.shapes$PC1max,gridPars=gridPar(pt.bg = "grey", pt.size = 0.5,link.lwd = 0.5,tar.pt.bg = "red", tar.pt.size = 0.8, tar.link.col = "red", tar.link.lwd = 0.8), method="points")
   
           #将平均形态和PC1最大值形态中的半标志点分别连接成不用颜色的曲线
               plot.semilandmark.lines (GPA.fit$consensus, curves = sliders, col = "grey")
               plot.semilandmark.lines (PCA.fit$pc.shapes$PC1max, curves = sliders, col = "red")
   
           #将模型调整至合适的视角后，保存视图
              rgl.postscript("pc1max-lateral.pdf",fmt="pdf",drawText = TRUE)
   
           #也可将模型导出为可交互的HTML文件，方便读者观察
           browseURL(
               paste("file://", writeWebGL(dir=file.path(tempdir(), "webGL"),
                                                   width=1500), sep="")
           )
   
   
6. 结果讨论
   PCA结果显示第一、二、三主成分反映了近 90%的形态差异，因此在它们构成的 坐标系里可以观察到样本主要的形态分异状况。其中，PC1反映了63.33%的形态差异，因此在对结果进行诠释时尤为重要。根据可视化的结果，PC1正值指示下颌背腹向纤细，而负值则指示下颌背腹向加深 (图10)。由图可见，喙部用于探食 (probing) 的鸟类在形态空间中明显位于PC1值为正值而PC3值为负值的位置 (图8和图9)，显示其下颌形态整体尤为纤细，且前端普遍上翘 (图10)；喙部用于 咬合破开 (Cracking/biting) 和撕扯 (Tearing) 的鸟类绝大部分位于PC1值为负值 的位置(图8和图9)，显示二者的下颌形态整体都背腹向更深，更加粗壮(图10)。此外，咬合破开型鸟类位于PC2值基本为正的位置，显示其下颌前端较尖，而撕扯型鸟类基本位于PC2值为负而PC3值为正的位置 (图8和图9)，显示其下颌前端具有突然的缩狭，且稍向下弯 (图10)。与此相对，喙部用于捡拾 (Grabbing/gleaning) 和啄食 (Pecking/grazing) 的鸟类占据的形态空间则有着大量重合，同时在PC1、PC2、PC3三个轴上都位于PC值为零左右的位置 (图8和图9)，显示其形态非常接近平均形态。
   如果此时将一个未知喙部功能的化石样本加入数据集进行计算，则可以根据其在形态空间中所处的位置对其喙部功能及相关生态习性进行初步的推测，但更为严格的检验则需要进一步的分析：PCA作为一种坐标转换工具，只能对生态习性和形态之间的关系提供一定的参考，但无法对相关假说进行严格的统计学检验。如果需要对假说进行进一步严格的统计学检验，则需要使用其他方法如phylogenetic multivariate ANOVA (Adams, 2014) 和phylogenetic flexible discriminant function analysis (Schmitz and Motani, 2011)。同时建议在可能的情况下，尽量将系统发育关系加入统计学分析中 (Felsenstein, 1985)，相关方法的应用参见(Schmitz and Motani, 2011; Close and Rayfield, 2012; Benson et al., 2017; Navalón et al., 2019; Choiniere et al., 2021)。
   

推荐阅读

1. 史宇坤. (2017). 形态测量学 (Morphometrics) 常用方法及其在微体古生物学中的应用. 微体古生物学报 34 (2): 179-191
   
2. Adams, D. C., Rohlf, F. J. and Slice, D. E. (2004). Geometric morphometrics: ten years of progress following the ‘revolution.’ Ital J Zool 71: 5-16.
   
3. Bjarnason, A. and Benson, R. (2021). A 3D geometric morphometric dataset quantifying skeletal variation in birds. MorphoMuseum 7: e125
   
4. Davies, T. G., Rahman, I. A., Lautenschlager, S., Cunningham, J. A., Asher, R. J., Barrett, P. M., Bates, K. T., Bengtson, S., Benson, R. B. and Boyer, D. M. (2017). Open data and digital morphology. P Roy Soc B-Biol Sci 284: 20170194.
   
5. Zelditch, M. L., Swiderski, D. L. and Sheets, H. D. (2012). Geometric Morphometrics for Biologists: a primer (Second Edition). Academic press. London.

致谢

This research is part of a project that has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 101024572.

参考文献

1. 史宇坤. (2017). 形态测量学 (Morphometrics) 常用方法及其在微体古生物学中的应用. 微体古生物学报 34 (2): 179-191
2. 向东进，李宏伟，刘小雅. 实用多元统计分析. (2014). 中国地质大学出版社.
3. Johnson, R. A. and Wichern, D. W. (2001). 实用多元统计分析. 清华大学出版社有限公司. 北京.
4. Adams, D. C. (2014). A method for assessing phylogenetic least squares models for shape and other high‐dimensional multivariate data. Evolution 68(9): 2675-2688.
5. Adams, D. C. and Otárola‐Castillo, E. (2013). Geomorph: an R package for the collection and analysis of geometric morphometric shape data. Methods Ecol Evol 4: 393-399.
6. Adams, D. C., Rohlf, F. J. and Slice, D. E. (2004). Geometric morphometrics: ten years of progress following the ‘revolution.’ Ital J Zool 71: 5-16.
7. Bardua, C., Felice, R. N., Watanabe, A., Fabre, A. C. and Goswami, A. (2019). A practical guide to sliding and surface semilandmarks in morphometric analyses. Integr Organism Biol 1: obz016.
8. Benson, R. B. J., Starmer‐Jones, E., Close, R. A. and Walsh, S. A. (2017). Comparative analysis of vestibular ecomorphology in birds. J Anat 231: 990-1018.
9. Bjarnason, A. and Benson, R. (2021). A 3D geometric morphometric dataset quantifying skeletal variation in birds. MorphoMuseum 7: e125.
10. Bookstein, F. L. (1989). Principal warps: thin-plate splines and the decomposition of deformations. IEEE T Pattern Anal 11: 567-585.
11. Gower, J. C. (1975). Generalized procrustes analysis. Psychometrika 40(1): 33-51.
12. Cardini, A. (2016). Lost in the other half: improving accuracy in geometric morphometric analyses of one side of bilaterally symmetric structures. Syst Biol 65: 1096-1106.
13. Cardini, A. (2017). Left, right or both? Estimating and improving accuracy of one‐side‐only geometric morphometric analyses of cranial variation. J Zool Syst Evol Res 55: 1-10.
14. Choiniere, J. N., Neenan, J. M., Schmitz, L., Ford, D. P., Chapelle, K. E. J., Balanoff, A. M., Sipla, J. S., Georgi, J. A., Walsh, S. A., Norell, M. A., Xu, X., Clark, J. M. and Benson, R. B. J. (2021). Evolution of vision and hearing modalities in theropod dinosaurs. Science 372: 610-613.
15. Close, R. A. and Rayfield, E. J. (2012). Functional morphometric analysis of the furcula in mesozoic birds. PLoS One 7: e36664.
16. Cooney, C. R., Bright, J. A., Capp, E. J. R., Chira, A. M., Hughes, E. C., Moody, C. J. A., Nouri, L. O., Varley, Z. K. and Thomas, G. H. (2017). Mega-evolutionary dynamics of the adaptive radiation of birds. Nature 542: 344-347.
17. Davies, T. G., Rahman, I. A., Lautenschlager, S., Cunningham, J. A., Asher, R. J., Barrett, P. M., Bates, K. T. et al. (2017). Open data and digital morphology. P Roy Soc B-Bioll Sci 284: 20170194.
18. Fabre, A., Goswami, A., Peigné, S. and Cornette, R. (2014). Morphological integration in the forelimb of musteloid carnivorans. J Anat 225: 19-30.
19. Felice, R. N. and Goswami, A. (2018). Developmental origins of mosaic evolution in the avian cranium. Proc Natl Acad Sci 115: 555-560.
20. Felsenstein, J. (1985). Phylogenies and the comparative method. Am Nat 125: 1-15.
21. Goswami, A., Watanabe, A., Felice, R. N., Bardua, C., Fabre, A. C. and Polly, P. D. (2019). High-density morphometric analysis of shape and integration: the good, the bad, and the not-really-a-problem.Integr Comp Biol 59: 669-683.
22. Kendall, D. G. (1984). Shape manifolds, procrustean metrics, and complex projective spaces. B Lond Math Soc 16: 81-121.
23. Bookstein, F. L. (1997). Morphometric tools for landmark data: geometry and biology. Cambridge University Press.
24. Navalón, G., Bright, J. A., Marugán‐Lobón, J. and Rayfield, E. J. (2019). The evolutionary relationship among beak shape, mechanical advantage, and feeding ecology in modern birds. Evolution 73: 422-435.
25. Gunz, P. and Mitteroecker, P. (2013). Semilandmarks: a method for quantifying curves and surfaces. Hystrix 24(1): 103-109.
26. Prum, R. O., Berv, J. S., Dornburg, A., Field, D. J., Townsend, J. P., Lemmon, E. M. and Lemmon, A. R. (2015). A comprehensive phylogeny of birds (Aves) using targeted next-generation DNA sequencing. Nature 526: 569-573.
27. Rohlf, F. J. and Slice, D. (1990). Extensions of the procrustes method for the optimal superimposition of landmarks. Systc Zool 39: 40-59.
28. Schmitz, L. and Motani, R. (2011). Nocturnality in dinosaurs inferred from scleral ring and orbit morphology. Science 332: 705-708.
29. Sutton, M., Rahman, I. and Garwood, R. (2014). Techniques for virtual palaeontology. John Wiley & Sons, Ltd. ISBN: 9781118591130.
30. Watanabe, A., Fabre, A. C., Felice, R. N., Maisano, J. A., Müller, J., Herrel, A. and Goswami, A. (2019). Ecomorphological diversification in squamates from conserved pattern of cranial integration. P Natl Acad Sci USA 116: 14688-14697.
31. Zelditch, M. L., Swiderski, D. L. and Sheets, H. D. (2012). Geometric Morphometrics for Biologists: a primer (Second Edition). Academic press. London.