在Python中使用的Procrustean伪地标方法，用于处理大量叶片形状数据

《Applications in Plant Sciences》：Procrustean pseudo-landmark methods in Python to measure massive quantities of leaf shape data

　　摘要

前提
在研究彼此不同的叶片形状时，以生物学和统计学上有意义的方式同时比较叶片的整体形状和叶片边缘的各个点可能会很困难。

方法
为了解决这个问题，我们提出了一种简单易用的叶片形状分析方法，该方法基于Jupyter Notebook和Python，利用伪地标（pseudo-landmarks）和广义Procrustes分析（generalized Procrustes analysis）来测量和比较任何叶片的形状。为了展示我们的分析方法，我们创建了一个包含来自八个实验数据集的真实叶片的仓库。

结果
通过使用我们的叶片仓库，我们展示了如何利用伪地标来比较同一物种内和不同物种间的叶片形状，并通过主成分分析（principal component analysis）等降维技术来进行比较；同时，我们还展示了如何通过线性判别分析（linear discriminant analysis）利用伪地标预测叶片形状。我们的叶片形状分析还展示了叶片在莲座状生长过程中形状的变化，揭示了叶片发育过程中的形态变化（phyllotaxy）。最后，我们展示了如何通过结合形状数据和遗传数据来研究叶片形状变异与遗传多样性之间的关系。

讨论
通过使用广义Procrustes分析和伪地标，我们的叶片形状分析成为了一种强大的工具，可以用于在多个生物学、生态学、进化和发育尺度上研究叶片的形状。生物学研究中最基本但也是最重要的方面之一就是形状变异现象（Thompson, 1963）。在植物科学中，研究通常会测量同一物种内不同时间、空间以及其他生物学相关轴线上花朵、根、茎和叶片形状的差异（Niklas, 1994）。传统的形状变异分析方法包括直接测量描述形状和大小的特征，以及测量长度、宽度和面积之间的异速生长关系（allometric relationships）。尽管采样有限，但手动测量叶片形状可以为发育遗传学（Hammond, 1941a, 1941b）、古植物学对气候的响应（Bailey and Sinnott, 1915）以及出于分类目的区分密切相关的群体提供基础性见解（Galet, 1979），从而证明了定量分析形状变异的力量。在数字和计算技术出现之前，这些测量都是手工完成的。根据叶片的形状和大小，这些测量可能既耗时又复杂，从而可能限制了可测量样本的数量。在过去的几十年里，我们进入了大数据科学和前所未有的计算能力时代。如今，许多植物器官的形状和大小可以通过计算机视觉自动标注，并使用一系列形态测量方法进行测量（Bucksch et al., 2017; Amézquita et al., 2020）。几何形态测量技术包括基于地标的方法（Bookstein, 1991）和椭圆傅里叶描述符分析（EFA）（Kuhl and Giardina, 1982）。对于基于地标的方法，闭合轮廓通过平移、旋转、缩放和反射函数进行对齐，以最小化Procrustes距离（Gower, 1975）。真正的地标是生物学上的同源点，基于地标的分析常用于测量和描述重要的花朵特征（例如侧瓣和改良花瓣），以及测量萼片、花瓣和其他生殖器官的角度（Savriama, 2018）。对于叶片地标，形状相似的叶片组可能在裂片数量、小叶数量或末端点数量等方面具有相似性。叶片地标的一个例子是葡萄叶形状的定量研究——ampelography，其中地标方法利用了所有葡萄叶都有五条主要叶脉的事实，从而在叶片之间创建了一个同源坐标系统（Chitwood et al., 2014）。因此，当存在生物学上的同源点时，基于地标的分析非常强大，因为形状的生物学和统计测量结果是相互对应的。然而，对于包含非同源点的形状或分辨率不足以捕捉许多细微叶片边缘特征的叶片图像，基于地标的分析可能并不总是可行的。对于侧生器官，如裂片或小叶的数量、角度和位置等特征可能会在不同物种间或随时间和空间变化（Balant et al., 2024）。不幸的是，缺乏同源性可能会使叶片形状的比较在计算上和生物学上变得困难甚至不可能，或者限制了对肉眼容易观察到的细微特征的测量。当同源性不存在或难以描述时，可以使用等距离放置的伪地标（Dujardin et al., 2014）在叶片的基部和尖端点之间进行比较。在Procrustean分析（Gower, 1975; Schindelin et al., 2012）中，通过使用地标和伪地标，可以对形状进行平移、缩放、旋转和翻转，以最小化两个形状之间的总体距离。在广义Procrustes分析（GPA）中，我们计算Procrustes距离并将所有形状叠加到一个任意样本上，从叠加的形状中生成一个新的平均值，并重复此过程，直到两个连续的平均形状之间的Procrustes距离低于任意选定的低阈值。一旦使用GPA计算出平均形状，就可以将所有样本叠加起来进行形态测量分析。通过基于地标或伪地标映射形状之间的差异，我们可以创建大量形状的连续表示。与几何形态测量类似，EFA也使用地标/伪地标来描述曲线。在EFA中，正弦和余弦函数以谐波级数的形式交替使用来描述闭合曲线或轮廓。几何形态测量、GPA和EFA技术已被包含在多个软件程序中，用于分析特定类型的形状，从而能够分析和生成大型形状数据集。在R语言中，流行的形状分析包包括shapes（Dryden and Mardia, 2016）用于GPA，而MOMOCS（Bonhomme et al., 2014）用于使用图像和EFA比较和可视化开放和闭合轮廓。在Python中，通常使用OpenCV程序（Bradski, 2000）从图像中提取轮廓。此外，许多形状分析可以在ImageJ/Fiji中完全完成（Schindelin et al., 2012），在那里可以测量、勾勒和提取形状，以便进一步分析。相比之下，随着研究叶片形状变异的增多，也开发了许多针对叶片的图像识别和分割程序。广泛使用的程序如LAMINA（Bylesj? et al., 2008）、MorphoLeaf（Biot et al., 2016）、MowJoe（Failmezger et al., 2018）、LeafletAnalyzer（Liao et al., 2017）和MuLES（Bowman et al., 2023）利用图像识别、自动分割和分类来比较叶片形状变异。然而，由于这些程序的“即用型”特性，它们的定制程度有限，限制了可以分析的非模型植物物种的数量和类型。因此，随着分析和生物学问题的计算复杂性在图像类型、物种和时间点上的不断增加，对新且灵活的分析方法的需求也在不断增长（Bylesj? et al., 2008; Liao et al., 2017; Failmezger et al., 2018; Bowman et al., 2023）。在这里，我们介绍了使用Python编写的几何形态测量方法来比较和描述具有或没有同源性的叶片形状差异的方法，并在易于理解和可复制的Jupyter Notebook中进行了展示（Kluyver et al., 2016）。我们的方法利用了已经可用于研究叶片形状变异的构建块，并在Python中应用了伪地标和GPA。我们的叶片形状分析基于这样一个事实：所有叶片在其尖端和基部都包含可比较的同源点。在这项研究中，我们提出了一种快速、简单且成本效益高的方法，该方法具有高度定制性，几乎可以用于任何叶片形状。我们展示了如何使用伪地标和Procrustean叶片形状分析来快速测量大量叶片，创建叶片形状的理论表示，通过降维技术可视化关系，并使用遗传和发育模型展示遗传和发育因素对叶片形状的连续影响。我们使用了来自全球不同研究物种的多样化数据集，并结合了先前发表的数据，以证明伪地标和GPA在多种叶片形状上的有效性。

叶片形状分析所需和可选数据
在这项研究中，我们将数据分为给定数据集所需的和可选的信息。所需信息包括分配伪地标、旋转每片叶片以及在GPA过程中使用伪地标比较叶片所需的数据。所需数据包括叶片轮廓——可以是二值图像（例如.jpg、.png、.tif或.bmp）或包含XY笛卡尔坐标的文本文件（.txt）。每片叶片的尖端和基部的XY坐标也是必需的；这些地标用于旋转和定位每片叶片，使尖端位于顶部，基部位于底部。最后，还需要每厘米像素数（px_cm）来缩放数据集中的每片叶片。每厘米像素数通常是通过在叶片图像扫描中加入尺子并使用ImageJ测量每1厘米的像素数量来计算的。关于每个数据集如何成像、分割和标记叶片的描述，包括使用的硬件和软件，可以在附录S1中找到（详见本文的支持信息）。每个数据集的可选信息包括识别数据集中每片叶片的人口统计、位置和遗传数据所需的信息，但即使没有这些信息也可以进行叶片形状分析。对于每片叶片，可选信息可以包括物种名称；种群和基因型数据；遗传和转录组数据；位置数据如纬度、经度、海拔、温度和降水量；生长条件；以及测试处理。在进行这项分析时，我们创建了一个最终文档（.csv），其中包含了所需的尖端、基部和px_cm信息，以及可选的数据集和基因型信息（附录S2）。因此，这种分析具有高度的定制性，可以用来沿着任何额外的轴检查数据，以回答多种生物学问题。

我们收集了来自全球合作者的八个独特数据集来测试我们的叶片形状分析，最终得到一个包含总共n=5129片叶片和n=48个物种的叶片数据集（附录S1）。Malus Mill. spp.数据集包含了最多的叶片（n=1560片）（Migicovsky et al., 2017）。其次是Quercus L. spp.（n=1142片）、Arabidopsis thaliana (L.) Heynh.（n=771片）（Seren et al., 2017）、Lobelia L. spp.（n=637片）、Cannabis sativa L.（n=340片）（Balant et al., 2024）、Beta vulgaris L.（n=269片）、Erythroxylum coca Lam.（n=260片）和Capsella bursa-pastoris (L.) Medik.（n=150片）（附录S1）。有三个数据集包含多个物种：Malus spp.（n=35片）、Lobelia spp.（n=10片）和Quercus spp.（n=3片）。B. vulgaris、C. bursa-pastoris、A. thaliana、Lobelia spp.和C. sativa的数据集的叶片图像是从共同花园（生长室或温室）环境中收集的，而Malus spp.数据集的叶片图像是从美国农业部（USDA）的两个野外采集点收集的，Quercus spp.数据集的叶片则来自德克萨斯州和俄克拉荷马州的公共和私人土地以及同一地区多年收集的标本（附录S1）。E. coca数据集的叶片是由哥伦比亚国家警察在反毒品局建立的花园中种植和收集的（附录S1）。除了上述叶片外，我们还分析了来自Passiflora spp.的3299片叶片和来自Solanum L. spp.（番茄）的11,276个小叶（Chitwood et al., 2013; Chitwood and Otoni, 2017a, 2017b），以展示这些方法如何应用于叶片发育的连续模型和高分辨率的叶片特征遗传映射。

方法
叶片形状分析所需和可选数据
在这项研究中，我们将数据分为给定数据集所需的和可选的信息。所需信息包括分配伪地标、旋转每片叶片以及在GPA过程中使用伪地标比较叶片所需的数据。所需数据包括叶片轮廓——可以是二值图像（即.jpg、.png、.tif或.bmp）或包含XY笛卡尔坐标的文本文件（.txt）。每片叶片的尖端和基部的XY坐标也是必需的；这些地标用于旋转和定位每片叶片，使尖端位于顶部，基部位于底部。最后，还需要每厘米像素数来缩放数据集中的每片叶片。每厘米像素数，记为“px_cm”，通常是通过在叶片图像扫描中加入尺子并使用ImageJ测量每1厘米的像素数量来计算的。关于每个数据集如何成像、分割和标记叶片的描述，包括使用的硬件和软件，可以在附录S1中找到（详见本文的支持信息）。每个数据集的可选信息包括识别数据集中每片叶片的人口统计、位置和遗传数据所需的信息，但即使没有这些信息也可以进行叶片形状分析。对于每片叶片，可选信息可以包括物种名称；种群和基因型数据；遗传和转录组数据；位置数据如纬度、经度、海拔、温度和降水量；生长条件；以及测试处理。在进行这项分析时，我们创建了一个最终文档（.csv），其中包含了所需的尖端、基部和px_cm信息，以及可选的数据集和基因型信息（附录S2）。因此，这种分析具有高度的定制性，可以用来沿着任何额外的轴检查数据，以回答多种生物学问题。

我们收集了来自全球合作者的八个独特数据集来测试我们的叶片形状分析，最终得到一个包含总共n=5129片叶片和n=48个物种的叶片数据集（附录S1）。Malus Mill. spp.数据集包含了最多的叶片（n=1560片）（Migicovsky et al., 2017）。在Malus数据集之后，按叶片数量从多到少排序，我们的最终数据集还包括Quercus L. spp.（n=1142片）、Arabidopsis thaliana (L.) Heynh.（n=771片）（Seren et al., 2017）、Lobelia L. spp.（n=637片）、Cannabis sativa L.（n=340片）（Balant et al., 2024）、Beta vulgaris L.（n=269片）、Erythroxylum coca Lam.（n=260片）和Capsella bursa-pastoris (L.) Medik.（n=150片）（附录S1）。有三个数据集包含多个物种：Malus spp.（n=35片）、Lobelia spp.（n=10片）和Quercus spp.（n=3片）。B. vulgaris、C. bursa-pastoris、A. thaliana、Lobelia spp.和C. sativa的数据集的叶片图像是从共同花园（生长室或温室）环境中收集的，而Malus spp.数据集的叶片图像是从美国农业部（USDA）的两个野外采集点收集的，Quercus spp.数据集的叶片则来自德克萨斯州和俄克拉荷马州的公共和私人土地以及同一地区多年收集的标本（附录S1）。E. coca数据集的叶片是由哥伦比亚国家警察在反毒品局建立的花园中种植和收集的（附录S1）。除了上述叶片外，我们还分析了来自先前研究的Passiflora spp.的3299片叶片和Solanum L. spp.（番茄）的11,276个小叶（Chitwood et al., 2013; Chitwood and Otoni, 2017a, 2017b），以展示这些方法如何应用于叶片发育的连续模型和高分辨率的叶片特征遗传映射。

叶片形状分析的总体过程包括叶片成像和分割（取决于数据集）、创建二值叶片图像或XY坐标叶片轮廓、每片叶片的方向和缩放、使用GPA比较所有叶片，以及使用降维技术和发育或遗传建模进行进一步分析。在这项研究中，每个数据集都包含了初始分析所需的信息，以及数据集内进一步分析所需的信息。对于所有数据集，我们要求提供叶片图像，可以是二值图像和/或包含XY笛卡尔坐标的叶片轮廓。因此，每个研究小组根据各自的标准为每个物种选择了叶片图像（附录S1）。我们要求提供每张叶片图像的像素数与厘米数的比例信息。每个小组还提供了每片叶片的顶端和基部的XY坐标以便于定位。对于包含物种或基因型信息的数据集，这些信息也被包含在它们的元数据中。在我们的形状分析中，我们首先分别分析了每个数据集，如下所述（在“伪地标插值”部分）。我们在每片叶片的边缘等距离放置了100个伪地标，从而形成了叶片的轮廓。然后，我们为每个数据集保存了一个定向的、重新缩放的并带有地标的数组，并将所有数据集合并成一个数组（以厘米为单位进行缩放）。这个合并后的数组被用来进行广义Procrustes分析（GPA）。

提取单个叶片和地标
对于每个数据集，我们使用ImageJ/Fiji软件选择了叶片并进行了地标标记，除了Lobelia属的叶片，这些叶片是使用WinRhizo（Arsenault等人，1995年）软件进行成像和地标标记的，这是根据研究人员的偏好决定的。B. vulgaris、C. sativa和C. bursa-pastoris的数据集包括了整株植物的叶片，包括莲座叶和茎生叶（如果适用的话）（附录S1）。Malus属、A. thaliana、E. coca和Quercus属的数据集包含了从每株植物上收集的特定数量的叶片。我们使用点选择工具在ImageJ中手动标记了每片叶片的基部和顶端。每片叶片的基部和顶端都被保存为笛卡尔（XY）坐标，并保存在一个单独的.txt文件中。每片叶片的元数据被包含在一个包含所有叶片坐标的合并.csv文件中。

伪地标插值
为了比较不同物种、基因型、地点和时间的叶片形状，我们使用了伪地标来进行比较。在这个过程中，我们首先将叶片形状.csv文件和所有图像文件（作为.jpg文件）按数据集上传到Jupyter Notebook（Kluyver等人，2016年）中。为了从叶片顶端开始沿着叶片边缘放置伪地标，所有图像都被读取为灰度图像，并使用OpenCV软件（Bradski，2000年）选择了最大对象（即叶片）的轮廓。利用OpenCV中的findContours()函数，我们提取了每张叶片图像的XY坐标和面积。然后我们使用插值函数将叶片坐标转换为高分辨率的点数。在这个上下文中，插值函数利用从叶片轮廓生成的已知真实点集来估计每片叶片的1000个高分辨率点。接下来，我们找到了每片叶片的基部和顶端地标之间的欧几里得距离（即线段的长度），并对基部和顶端点进行了重新索引。然后我们将基点重置为零，并重新索引了伪地标，创建了一个包含所有叶片坐标的单一数组。接下来，我们使用重新索引的地标对每片叶片的左右两侧进行了插值，每侧估计了50个伪地标，总共100个伪地标。然后我们根据每张叶片图像提供的像素数与厘米数的比例，对每片叶片进行了定向和重新缩放。最后，我们将每个定向和重新缩放的叶片保存为一个数组。对于每个数据集，我们保存了一个单独的定向和重新缩放的叶片数组。我们创建了一个名为“cm_arr”的最终数组，其中包含了所有数据集所需和可选的数据。向每张图像添加伪地标的过程是这个流程中计算量最大的部分。然而，我们使用一台普通的计算机在大约10分钟内处理了超过5000片叶片，这表明该流程能够高效地分析5000片叶片。

广义Procrustes分析（GPA）
为了计算Procrustes均值，我们首先选择了伪地标的数量和维度。然后我们使用gpa_mean函数计算了平均叶片形状，该函数包括了我们的叶片形状数组、伪地标的数量和维度数量。这个平均叶片是通过将每片叶片与所有其他叶片进行比较，并使用伪地标来创建所有成像叶片的平均值而生成的。然后我们将所有叶片对齐到新生成的平均叶片上，从而创建了一个所有对齐叶片的数组。为了生成每个数据集的平均叶片，我们使用了相同的方法，但是将“cm_arr”数组替换为数据集数组。

叶片尺寸分析
对于每片叶片，我们计算了宽度、长度、面积（以平方厘米为单位）、 solidity（即面积与凸包面积的比率）、圆形度（即面积与周长的比率）以及纵横比（即宽度与长度的比率）。所有这些值都被保存在mdata.csv文件中。

主成分分析
为了对所有使用100个伪地标对齐的叶片进行主成分分析（PCA），我们将Procrustes对齐的数组从四维重新调整为二维，计算了我们所需的主成分（PCs）的数量，执行了PCA函数，并将PC1和PC2保存在原始.csv文件中。这个PCA的载荷是100个伪地标的组合。我们使用每片叶片的PC1和PC2值计算并可视化了数据集中所有叶片的形态空间，然后通过将所有叶片绘制在形态空间散点图上来可视化所有叶片在PC空间中的位置。

线性判别分析
为了对最终数据集中的所有叶片进行线性判别分析（LDA），我们首先根据数据集创建了一个新的数据框。我们从mdata.csv文件中导入了数据集标签，并将地标数据作为x值，数据集标签作为y值。然后我们拟合了LDA模型，检索了LDA的缩放和系数，并进行了基因型预测。接下来，我们根据原始.csv文件和预测的基因型值创建了一个混淆矩阵。然后我们使用scikit-learn版本1.8.0中的Stratified K-fold交叉验证包（Pedregosa等人，2011年）对我们的模型进行了验证，默认K=5，重复三次，随机状态为42。最后，我们将混淆矩阵保存在新的数据框中，并使用LD1和LD2在LD空间中可视化了混淆矩阵，比较了真实数据集和预测数据集以及所有叶片。

我们使用Passiflora属的数据集来展示叶片形状从植物基部到顶端的变化，该数据集包含了40个物种的先前发表的叶片形状数据（Chitwood和Otoni，2017a，2017b）。首先，我们定义了间隔的数量（即所有莲座叶中的最小叶片数量）；建议使用三个或更多节点的序列进行建模。然后我们找到了数据集中的伪地标总数和个体植物总数，并创建了一个包含间隔、植物数量和伪地标数量的数组。接下来，我们创建了空列表来存储物种和形态类型信息。叶片形状可以使用元数据中已经定义的任何类别的信息来进行建模。最后，我们首先找到了每株植物每片叶片的相对节点数，然后使用NumPy包中的polyfit函数（Harris等人，2020年）来建模叶片形状随发育过程的变化。该函数使用XY坐标拟合一个多项式。polyfit函数将相对节点作为x值，每个节点的100个伪地标作为y值，以及可能的最小多项式阶数（polyorder）。对于这个数据集，我们使用了二阶多项式。然后我们对每个节点的多项式值进行了建模，创建了每株植物40片建模叶片，每片叶片有100个伪地标，并将所有建模叶片保存为一个数组，所有物种和形态类型信息保存为一个数据框。我们通过按节点数量递增的顺序绘制每片建模叶片来可视化每个数据集和形态类型的40片建模叶片，并按形态类型进行了着色。然后我们通过将每个模型的间隔投影回形态空间PCA来可视化每株植物的发育模型。每个模型表示为一条曲线，其中单个叶片在曲线上连接起来，并根据它们在序列中的位置进行了着色。我们为每个序列中的建模叶片位置指定了一个通用的颜色映射，使用viridis颜色刻度（Garnier等人，2024年），其中黄色代表位于或接近植物基部的叶片，紫色代表位于或接近植物顶端的叶片。

我们使用了之前发表的S. lycopersicum数据集，其中包含了76条引入系（ILs）的叶片数据，每个引入系都包含来自S. pennellii的小规模引入。对于每个引入系，我们重新索引了每个IL，以表示每个引入系的存在与否，这些引入系分布在S. lycopersicum基因组中。类似于发育过程中的建模，我们使用伪地标对所有数据集中的右侧叶片、左侧叶片和顶端小叶进行了Procrustes对齐。在为右侧叶片、左侧叶片和顶端叶片生成了平均叶片后，我们计算了每个伪地标到每个小叶的平均叶片的欧几里得距离（记为位移）。然后我们计算了每个伪地标的位移与每个引入系的存在/缺失之间的Spearman等级相关系数，生成了一个P值。然后我们根据每个伪地标的相关引入系及其P值可视化了位移的程度。

GPA利用伪地标揭示了叶片形状之间的关系
使用100个伪地标和GPA进行的平均叶片比较揭示了一个缺乏边缘特征（如锯齿或裂片）的共识叶片，这是由于平均作用的结果（图1A）。当特征在样本之间无法与特定的伪地标对齐时，就会发生这种平均效应，这也可能归因于每个数据集内叶片之间的显著变异（图1B）。整体平均叶片与每个数据集的平均叶片之间的比较显示，每个数据集的平均叶片都是独特的（图1B）。C. bursa-pastoris、Lobelia属和Quercus属的数据集包含的叶片具有高水平的非同源内部变异，导致形状差异被平均掉，在这三个数据集的平均叶片中不可见。相比之下，Malus属、A. thaliana、E. coca和B. vulgaris的数据集的平均叶片在每个数据集中都类似于真实叶片，并包含锯齿等特征（图1B）。特别是，C. sativa的平均叶片是一个简单的掌状叶，其平均值有三个明显的裂片，这是由于这些特征在叶片上的位置相似；而B. vulgaris的平均叶片是一个简单的叶片，具有宽大的箭形叶片和明显的叶柄。

图1展示了所有数据集中叶片形状的多样性，包括每个数据集的平均叶片。（A）通过Procrustes对齐生成的整个数据集的平均叶片（左），以及叠加在每个单独叶片轮廓上的相同平均叶片（右）。（B）对于每个物种数据集，包括从每个物种数据集中提取的100个伪地标的10个随机叶片轮廓的对齐。每个叶片的顶端用橙色圆圈标记，基部用蓝色圆圈标记。每个叶片旁边的垂直黑线代表1厘米的刻度尺。对于每个物种数据集，整个数据集的平均叶片由黑色轮廓表示，该轮廓叠加在彩色物种平均叶片上。形状描述符揭示了基于总体形状的叶片之间的关系。

在我们的分析中，可以直接使用形状描述符（如圆形度）来测量形状，圆形度定义为...圆形度是衡量叶片形状波动的指标，常用于量化和描述叶片边缘的裂片程度。圆形度描述了叶片相对于完美圆形的圆度；因此，圆形度较低的叶片具有接近0的圆形度分数，是最不圆的，而圆形度较高的叶片具有接近1的分数，更加圆润。圆形度值会随着叶片的波浪状结构（如裂片和锯齿）、长而逐渐变细的叶柄或细长的叶片形状而降低。整个数据集的圆形度范围在0.1到0.85之间，平均圆形度为0.45，涵盖了所有可能的叶片形状（图2A、B）。在检查每个数据集内的圆形度分布时（图2B），C. bursa-pastoris（图2C）、Lobelia属植物、Quercus属植物、Malus属植物和A. thaliana的数据集显示出最大的圆形度变化（0.07–0.86），其次是B. vulgaris和E. coca（0.41–0.85），最后是C. sativa（0.07–0.33）。Cannabis sativa的叶片整体上具有最低的圆形度，并且裂片最为明显。在这里，圆形度特别受到E. coca叶片细长形状和逐渐变细的叶片基部的影响（图2D）。总体而言，当我们考虑数据集内部的形状变化时，C. bursa-pastoris、Lobelia属植物、Quercus属植物和Malus属植物的叶片彼此之间的相似度最低，而C. sativa、B. vulgaris和E. coca的叶片则最为相似。

叶片形状的多样性可以通过圆形度来描述。（A）每个物种数据集中叶片圆形度值的直方图。（B）按圆形度着色的形态空间主成分分析（PCA），由整个数据集的叶片组成。（C）按圆形度着色的Capsella bursa-pastoris数据集的形态空间PCA。（D）按圆形度着色的Erythroxylum coca数据集的形态空间PCA。在每个PCA中，叶片沿着连续的圆形度谱进行着色，低圆形度用浅粉色表示，中等圆形度用浅绿色到深绿色表示，高圆形度用深紫色表示。形态空间PCA揭示了基于叶片各部分的叶片之间的关系。

与直接使用形状描述符测量整体形状相比，我们还可以使用伪地标和降维技术来确定和预测叶片边缘的哪些部分对叶片之间的整体变化贡献最大。通过PCA，我们可以通过检查PC空间中形状的重叠来了解多个轴上的形状关系，从而可视化形态空间。在形态空间PCA中，伪地标组合的PC载荷可以揭示叶片边缘的哪些部分对每个PC轴上的叶片差异贡献最大。此外，通过进行逆PCA，我们可以重建对整体差异贡献最大的形状。在整个数据集的形态空间PCA中，叶片的宽度沿着PC1从左到右变化，左侧的叶片较薄，而右侧边缘的叶片较宽，尤其是那些具有多个裂片和小叶的叶片（附录S3）。这种叶片纵横比的变化通过PC1中的一个值来捕捉，解释了我们数据中48.6%的可观察变化。此外，叶片基部叶片层的数量和明显叶柄的存在沿着PC2从上到下变化。位于形态空间顶部的叶片没有明显的叶柄，并且包含最多的基部叶片层，而位于形态空间底部的叶片具有明显的叶柄，并且基部叶片层较少（附录S3）。叶柄的明显程度通过PC2解释了我们数据中11.3%的变化。PCA形态空间揭示了我们数据集之间的有趣关系（图3A）。Cannabis sativa的叶片与其他所有数据集分离。在形态空间中，C. sativa占据了由具有多个裂片的叶片代表的区域，而其他七个数据集占据了其余部分，这些数据集的叶片具有大量的裂片变化。相比之下，Quercus属植物和Malus属植物的叶片形成了具有部分重叠的明显组（图3B），而所有其他数据集的叶片在PC空间中显示出显著的重叠（图3B）。当我们检查每个数据集中单个叶片在形态空间中的分布时，出现了两个主要趋势。在数据集内部，单个叶片在PC空间中涵盖了较大的范围，表明数据集内的叶片形状变化大于数据集之间的变化。在各自的数据集中，Quercus属植物、E. coca和B. vulgaris的叶片彼此之间的相似度最低，而C. sativa、B. vulgaris和E. coca的叶片彼此之间的相似度最高。

所有叶片的形态空间PCA。投影的形态空间包括通过取PCA的逆得到的特征叶片。这些特征叶片是根据每个叶片图像的100个伪地标创建的。（A）整个数据集的形态空间PCA。（B）每个物种数据集的形态空间PCA与整个形态空间的对比图。点颜色对应于每个物种数据集。LDA预测了沿生物学相关轴的叶片之间的关系。

我们可以使用对叶片形状差异贡献最大的叶片部分，通过LDA来预测和解释沿生物学重要轴的叶片形状之间的关系。LDA与PCA类似，但LDA的轴不是基于解释数据整体变化最大的线性组合来定义的，而是基于最能区分定义组的特征线性组合来定义的。由于LD轴是由指定的组定义的，LDA可以用于根据多个连续变量（如Procrustean对齐的坐标值）来建模分类组身份。在使用LDA时，我们的目标是使用100个伪地标来预测叶片属于哪个数据集。值得注意的是，由于每个数据集中的叶片数量不等，我们使用了分层K折交叉验证来为训练组和测试组创建均匀的样本。分层采样有助于提高真正例预测的可能性，确保每个数据集在训练组和测试组中都有比例上的代表。如果使用随机采样而不是分层采样，较小的数据集（例如C. bursa-pastoris数据集[n=150]）可能不会得到比例上的代表。当我们进行LDA时，我们的数据集达到了87%的平均准确率，其中511片叶片被错误预测，而5129片叶片中有4618片被正确预测。我们的混淆矩阵显示，正确预测的叶片大多来自Malus属数据集（图4A）。C. bursa-pastoris的叶片预测正确率最低（n=67/150），并且经常被错误地预测为A. thaliana的叶片（n=60）（图4A）。在其他叶片数据集中，97%的E. coca叶片、96%的C. sativa叶片、92%的B. vulgaris叶片、91%的Quercus属叶片和80%的Lobelia属叶片在LDA中被正确预测（图4A）。与PCA中的PC类似，LD是我们预测变量的线性组合，即100个伪地标。当我们可视化预测结果时，我们看到数据集沿着LD1分离，C. sativa的叶片与其他所有数据集有明显的区分（图4B）。沿着LD2，虽然我们看到Malus属和Quercus属叶片之间有明显的分离，但我们也在簇的中心看到其余数据集之间的显著重叠（图4B、C）。虽然我们的结果表明大多数数据集的叶片可以使用伪地标可靠地预测，但我们的模型可能对C. bursa-pastoris叶片数量较少以及C. bursa-pastoris和A. thaliana叶片形状之间的相似性敏感。通过结合两种降维技术的结果，我们能够建模数据集之间的关系。当我们综合考虑PCA和LDA的结果时，我们看到C. sativa的叶片始终分离成一个明显的组，Malus属和Quercus属叶片在PC和LD空间中都有共同点，但可以区分成不同的组，而其余数据集在基于100个伪地标和圆形度的叶片形状上显示出显著的重叠。

所有叶片的线性判别分析。（A）基于100个伪地标的每个物种数据集中实际叶片和预测叶片的混淆矩阵。（B）每个物种数据集中预测叶片的散点图。（C）从LD1到LD5的每个物种数据集中预测叶片的散点图。通过建模叶片的发育过程和遗传数据，可以更深入地了解形状变化的机制。

除了PCA和LDA这样的降维技术外，还可以使用叶片沿茎的生长位置来建模连续生长系列的叶片。对连续生长系列的叶片进行建模可以了解叶片形状在发育过程中的变化，这与基因表达的变化直接相关。这种连续建模的一个例子可以在Passiflora数据集的叶片中看到。为了展示从茎基部（最早出现的叶片）到顶端（最晚出现的叶片）的叶片形状变化，我们使用了之前发表的数据（Chitwood和Otoni, 2017a, 2017b）来根据相对节点位置对叶片进行建模。总共分析了3299片叶片。每株植物的叶片/节点总数通常不同，这可能导致相同节点数的叶片处于不同的发育阶段。因此，为了比较生长顺序和叶片发育，必须使用相对节点数（即从茎基部开始的节点数除以总节点数）来上下文化每株植物的叶片生长顺序。我们将每个伪地标的x和y坐标值作为每个茎的相对节点数的函数来建模叶片形状的发育过程，使用二次多项式。使用这些连续模型，我们然后按照从茎基部到顶端的顺序可视化了每个Passiflora物种的40片叶片，应用了为Passiflora叶片形态类型开发的颜色方案（图5A）（Chitwood和Otoni, 2017b）。在Chitwood和Otoni（2017b）中，形态类型是根据沿维管系统的Procrustes调整后的地标定义的，使用窦和裂片以及椭圆形傅里叶描述符来定义叶片轮廓。这些分析的综合结果揭示了七种Passiflora叶片形态类型（图5A）。形态类型A、E、F和G在发育过程中的叶片形状变化很小。相比之下，形态类型B、C和D在发育过程中的变化较大，靠近基部的叶片比靠近顶端的叶片更均匀。当我们使用PCA形态空间比较单个Passiflora叶片时，我们看到形态类型之间有显著的分离（图5B）。形态类型E、F和G有重叠，代表更完整的叶片；形态类型A和B在PC空间中有重叠，其特征叶片（或理论叶片）宽度大于长度，并且没有明显的叶片尖端；形态类型C单独分离，由三裂叶片形态表示；而形态类型D没有与其他形态类型明显分离，可以在所有三个簇中找到（图5B）。当我们使用连续曲线在PC空间中比较Passiflora叶片的发育模型时，我们看到形态类型分为三个不同的组（图5C）。一般来说，每种形态类型和物种的最初叶片裂片较少，更完整，起源于形态空间的中心（图5C）。后期发育的叶片表现出更大的形状多样性，叶片边缘结构有变化（图5C）。通过建模发育过程中的形状变化，我们可以了解叶片形状在发育过程中的变化时间和位置。

除了PCA和LDA这样的降维技术外，还可以使用叶片沿茎的生长位置来建模连续生长系列的叶片。对连续生长系列的叶片进行建模可以了解叶片形状在发育过程中的变化，这与基因表达的变化直接相关。这种连续建模的一个例子可以在Passiflora数据集的叶片中看到。为了展示从茎基部（最早出现的叶片）到顶端（最晚出现的叶片）的叶片形状变化，我们使用了之前发表的数据（Chitwood和Otoni, 2017a, 2017b）根据相对节点位置对叶片进行建模。总共分析了3299片叶片。每株植物的叶片/节点总数通常不同，这可能导致相同节点数的叶片处于不同的发育阶段。因此，为了比较生长顺序和叶片发育，必须使用相对节点数来上下文化每株植物的叶片生长顺序（即从茎基部开始的节点数除以总节点数）。我们将每个伪地标的x和y坐标值作为每个茎的相对节点数的函数来建模叶片形状的发育过程，使用二次多项式。使用这些连续模型，我们然后按照从茎基部到顶端的顺序可视化了每个Passiflora物种的40片叶片，应用了为Passiflora叶片形态类型开发的颜色方案（图5A）（Chitwood和Otoni, 2017b）。在Chitwood和Otoni（2017b）中，形态类型是根据沿维管系统的Procrustes调整后的地标定义的，使用窦和裂片，并使用椭圆形傅里叶描述符来定义叶片轮廓。这些分析的综合结果揭示了七种Passiflora叶片形态类型（图5A）。形态类型A、E、F和G在发育过程中的叶片形状变化很小。相比之下，形态类型B、C和D在发育过程中的变化较大，靠近基部的叶片比靠近顶端的叶片更均匀。当我们使用PCA形态空间比较单个Passiflora叶片时，我们看到形态类型之间有显著的分离（图5B）。形态类型E、F和G有重叠，代表更完整的叶片；形态类型A和B在PC空间中有重叠，其特征叶片（或理论叶片）宽度大于长度，并且没有明显的叶片尖端；形态类型C单独分离，由三裂叶片形态表示；而形态类型D没有与其他形态类型明显分离，可以在所有三个簇中找到（图5B）。当我们使用连续曲线在PC空间中比较Passiflora叶片的发育模型时，我们看到形态类型分为三个不同的组（图5C）。一般来说，每种形态类型和物种的最初叶片裂片较少，更完整，起源于形态空间的中心（图5C）。后期发育的叶片表现出更大的形状多样性，叶片边缘结构有变化（图5C）。通过建模发育过程中的形状变化，我们可以了解叶片形状在发育过程中的变化时间和位置。除了我们合作数据集和之前发表的西番莲数据中包含的其他叶物种外，还使用了番茄（Solanum lycopersicum）的叶子来研究叶形状的遗传和发育基础（Chitwood等人，2013年）。Chitwood等人（2013年）使用了76个S. lycopersicum导入系（ILs）的叶形状数据，这些导入系中包含了S. pennellii基因组的小片段，这些片段以序列方式整合到了S. lycopersicum基因组中，揭示了末端叶、近端左侧叶和远端右侧叶在形状上的遗传差异。总共分析了11,276个小叶。由于每个导入系都包含一个明确的S. pennellii基因组导入，因此导入系与S. lycopersicum背景之间的显著统计差异表明了其遗传基础（即导入）。Chitwood等人（2013年）表明，可以将数量性状位点（QTLs）映射到整体形状描述符或小叶形状的PC总结上。由于每个导入系只包含一次栽培番茄和野生番茄之间的基因组导入，因此这些导入系在此分析中起到了QTL的作用。在这里，我们分别连续映射每个小叶上的每个点；这使我们能够识别出与小叶形状连续变化相关的QTL。为了将伪地标转换为连续的性状值，我们计算了样本中的伪地标与其在整体平均值上的对应位置的欧几里得距离，其中超出平均值的距离为正，低于平均值的距离为负。我们使用Spearman等级相关性来关联导入的存在与否与每个伪地标从末端叶、远端右侧叶或远端左侧叶的整体平均值上的位移。我们可以将每个伪地标从其相应平均小叶上的位移程度可视化为热图，沿着小叶（叶区域）绘制伪地标，并根据遗传图谱（染色体导入）和交互作用的统计强度及其方向进行着色（图6B）。我们连续映射的结果揭示了与之前描述的QTLs相对应的末端叶、右侧叶和左侧叶的区域（Chitwood等人，2013年）。最重要的是，我们现在可以识别出之前描述的QTLs（由导入系表示）作用于哪些小叶的区域。IL2-1在近端小叶基部表现出更窄的叶子，超出了S. lycopersicum的特征。IL9-1-2在末端叶和侧叶的远端表现出更窄的叶子，同样超出了S. lycopersicum的特征。IL4-3是最极端的番茄导入系之一，它在侧叶的远端区域表现出更宽的叶子，但在末端叶上也有一定程度的宽叶现象，类似于S. pennellii。IL5-4在末端叶上表现出更宽的叶子，在基部也表现出宽叶现象。因此，通过使用简单的相关性、伪地标和明确的遗传资源，我们能够模拟与番茄复叶特定部分相关的遗传结构。

图6：在Solanum物种中根据平均小叶形状与数量性状位点（QTL）之间的关系映射伪地标位置。（A）从左到右：远端左侧叶、末端叶和远端右侧叶。每个小叶的广义Procrustes分析（GPA）平均值以洋红色显示，叠加的小叶轮廓以黑色表示。（B）小叶伪地标的遗传基础，显示了远端左侧叶（顶部）、末端叶（中间）和远端右侧叶（底部）的结果。在左侧，每个小叶都标有分配到小叶轮廓上10个区域的伪地标。对于Benjamini–Hochberg程序假发现率（FDR）调整后的P值小于0.05的QTL，显示了一条垂直于相关平均伪地标位置的线段，并根据P值和方向进行着色，显示了特定小叶位置的遗传效应。在每种小叶类型旁边，显示了一个热图（根据-log base 10 Benjamini–Hochberg FDR调整后的P值进行着色），表示每个小叶区域（由y轴上的阿拉伯数字表示）和染色体上的遗传区间（由x轴上的罗马数字表示）。

**讨论**

**基于伪地标和Procrustes分析的优势**

如上所述，使用伪地标和基于Procrustes的叶形状分析的一个关键优势是能够包括高度变化和不规则的叶形状。包含地标的几何形态测量分析依赖于物种和形状之间的生物学同源性。此外，传统的地标通常数量有限，且往往是手工选择的。对于许多物种来说，实现真正地标所需的生物学同源性水平使得传统的几何形态测量分析变得极其困难。例如，C. bursa-pastoris的叶子在裂片数量、裂片角度和裂片深度上表现出极大的变化，使得在叶子之间找到同源性几乎不可能。相比之下，通过利用所有叶子共有的两个地标——基点和尖点——并在这些锚点之间使用等距放置的伪地标，我们不仅可以增加可比较的数据点数量，还可以对齐数据点以推断物种和形状之间的新同源点。对齐和比较形状的能力的重要性不容忽视。在模拟形状变化时，确定平均形状以及个体与平均值的位移对于注释和测量形状变化至关重要。使用线性建模，Procrustes距离和伪地标的组合可以用来确定形状之间的差异。此外，伪地标的组合显示了沿叶缘最负责叶子之间变化的位置，揭示了潜在的组织选择和表型可塑性位置，并创建了捕捉额外叶形状复杂性的特征，如果只考虑维度或其他摘要特征则会丢失这些信息。总体而言，使用基于伪地标和Procrustes的方法是一种简单、快速且信息丰富的叶形状变化测量方法。

**使用伪地标描述和测量叶形状的另一个优势是，伪地标可以用作PCA和LDA的输入，分别创建形态空间的可视化和预测模型。通过创建PCA形态空间，我们不仅可以重现连续定义叶子之间变化的理论叶形状，还可以测量和可视化我们数据集中叶子之间的关系，这些关系跨越多个正交轴。如前所述，PCA形态空间存在边界（Chitwood和Mullins，2022年），这些边界揭示了叶形状在发育方面的可能限制。这可以在远离数据集中心的形态空间边缘观察到，其中特征叶会以实际叶子从未出现的方式伸展或折叠。这种现象也出现在使用伪地标方法难以捕捉的叶形状上，例如非常薄的叶子或具有许多不规则裂片的叶子。尽管如此，这些困难的形状仍然代表了远离常见祖先叶型的极端叶形态。此外，应使用尽可能少的伪地标数量来定义叶形状，因为更多的比较点并不会提高形状分辨率。同样，对于LDA，通过使用伪地标的线性组合，可以使用那些最有助于形状分离的叶部分进行预测。LDA模型得出的分类也有助于理解单个叶形状之间的关系，彼此更相似的叶形状以及与整体平均叶形状更相似的叶形状更容易被混淆，而更独特和极端的叶形状则更容易被预测。**

**地标的生物学相关性和形状**

正如伪地标能够将叶缘上的特定位置精确转换为可以在所有叶子之间比较的特征，从而确定哪些特征最受遗传或发育因素的影响一样，它们也可以用来识别叶形态的功能方面。古植物学和共同花园研究已经充分证明，叶形状的变化影响了许多对植物适应性至关重要的过程。生理学研究表明，深裂的叶子具有较低的中柱组织与大型导脉的比例，这降低了水分在叶子中的传输阻力（Xu等人，2009年；Ding等人，2020年；Sakurai和Miklavcic，2021年）。在裂片较少的叶子中，较小的叶脉增加了水分阻力，因为水分必须在细胞之间横向传输（Ding等人，2020年；Sakurai和Miklavcic，2021年）。此外，小叶中的热对流比大叶更快，边界层的厚度随着叶长的增加而增加，导致热量散失的差异更大（Bednarz和van Iersel，2001年；Peppe等人，2011年）。在深裂的叶子中，热量必须穿越叶片的距离减小，导致裂片叶子的热量传递大于整片叶子（Zwieniecki等人，2004年）。因此，叶结构与功能之间的生物学权衡已有充分记录（de Boer等人，2016年；Gutschick，2016年）。因此，通过了解叶子在不同尺度上的变化部分，我们可以推断出叶形状变化可能导致的进化和生理功能。

**伪地标和Procrustes分析的未来应用**

随着伪地标和基于Procrustes的分析在叶形状研究中的广泛应用，这些分析的未来用途包括创建合成叶子（Balant等人，2024年；Chitwood等人，2025年）或理论特征叶，这些特征叶在越来越多的PC上采样，可以用来推断叶形状的未测量连续变化。随着使用自动化分割方法从共同花园和标本馆的图像叶子中测量更多叶子，我们将能够预测形状边界并在生物学相关的轴上可视化理论叶形状。通过创建发育模型，我们同样可以捕捉到连续的发育变化，先前的研究表明这种变化是加性的，并且与物种间的叶形状差异正交（Chitwood等人，2016a，2016b，2025年）。例如，通过模拟节点上的叶子，无论是茎还是莲座状叶，我们可以生成叶形状的连续模型，作为叶发育时间指数的函数，该指数测量两个连续叶子形成之间的时间间隔。叶发育时间指数根据叶形状等形态特征描述植物的年龄，也可以用来预测异形性和个体发育的轨迹。通过伪地标对叶形状的形态测量也可以用于高分辨率的遗传学研究。伪地标和维度减少输出（如PCs）可以直接用于全基因组关联研究或RNA-Seq分析。然而，正如我们在这里展示的，通过测量伪地标与其相应平均值之间的欧几里得距离位移，可以独立测试各个点的遗传基础。此外，Procrustes距离本身是两种形状之间整体相似性的度量，可以用来测量组内和组间的变化，相对于各自的GPA平均值。最后，伪地标以及这里讨论的任何其他测量组合都可以用于深度学习，其中伪地标可以用来创建人工神经网络，以更深入地理解叶形状模式识别。通过样本之间的距离定义的形态空间，结合定义分子谱空间的类似空间，为使用图学习从基因组信息预测表型提供了可能性。

**行动号召：构建所有叶形状的形态空间**

在这项研究中，我们展示了检查叶形状变化可以是简单、廉价且快速的。重要的是，通过使用所有叶子共有的基点和尖点，并在叶子的每一侧等距放置伪地标，任何叶子都可以与其他叶子进行整体相似性比较。因此，随着更广泛的植物科学界继续研究叶形状，我们鼓励对叶形状感兴趣的研究人员不仅将他们的分析和数据与这里介绍的分析技术结合使用，还将他们的数据与这里共享的叶形状结合起来，创建一个不断增长的叶形状形态空间。将分割算法应用于数字化的植物标本，有望使来自不同物种和环境的数百万片叶子数据得以共享，从而可以将特定研究中的叶片形状与所有叶片进行比较。深入研究遗传、发育、生理或进化过程的研究将能够将这些宝贵信息与所有叶片形状的总体情况相结合。现在是一个前所未有的时机，可以实现对植物叶片形态在遗传、发育、进化和功能方面的综合分析。有关如何提交您的叶片形状数据的信息，请参阅我们在GitHub上的说明（https://doi.org/10.5281/zenodo.17635690；请参阅数据可用性声明）。

**作者贡献**

A.H. 收集了每个数据集的实验数据，使用C. bursa-pastoris进行了实验数据采集，创建了分析结果，准备了包括Jupyter Notebook和GitHub仓库在内的叶片形状分析材料，并起草和编辑了手稿。D.H.C. 进行了初步的叶片形状分析，参与了后续分析的讨论，并审阅和编辑了手稿。A.C.、C.B.、S.H. 和 S.W. 生成了Lobelia spp.数据集中的叶片，并对手稿提供了反馈。Cl.H.、Ar.P. 和 Ch.C. 生成了Quercus spp.数据集中的叶片，并对手稿提供了反馈。J.H.、L.S. 和 Z.M. 生成了Malus spp.数据集中的叶片，并对手稿提供了反馈。A.Y.H. 和 R.U.C. 生成了A. thaliana数据集中的叶片，并对手稿提供了反馈。M.B. 和 Al.P. 生成了C. sativa数据集中的叶片，并对手稿提供了反馈。A.D. 生成了B. vulgaris数据集中的叶片，并对手稿提供了反馈。E.B.J.、Cl.H.、Cl.C.、E.A. 和 G.D. 生成了C. bursa-pastoris数据集中的叶片，并对手稿提供了反馈。Y.A.S.R. 和 E.P. 生成了E. coca数据集中的叶片，并对手稿提供了反馈。所有作者均批准了手稿的最终版本。

**致谢**

作者感谢美国国家科学基金会（项目编号2039489和2310356，资助A.Y.H.）和美国国立卫生研究院（项目编号R35 GM142829，资助E.B.J.）的资助。A.D. 参加了由美国国家科学基金会资助的“Plant Genomics @ MSU Research Experiences for Undergraduates (REU)”项目。我们感谢美国农业部农业研究服务局的Rachel Naegele为B. vulgaris数据集生成叶片数据并对手稿提供反馈。

**数据可用性声明**

本研究中使用的叶片轮廓、各个数据集的元数据，以及Jupyter Notebook和配套代码，均可在GitHub仓库（https://doi.org/10.5281/zenodo.17635690）中找到。