RRmorph - 3D Morphological Analyses with RRphylo

Marina Melchionna, Silvia castiglione, Pasquale Raia

Index

  1. RRmorph overview
  2. Preparing the data
  3. The analysis of shape
  4. rate.map
  5. conv.map
  6. interpol.mesh

1. RRmorph overview

RRmorph is a comprehensive toolkit designed to investigate the effects of evolutionary rates and morphological convergence on phenotypes. The package integrates methodologies from Phylogenetic Comparative Methods (specifically, RRphylo, Castiglione et al. 2018) and Geometric Morphometrics to analyse shape evolution.

The functions embedded in RRmorph are landmarks-based. They allow to quantify and visualize the impact of rate variation and morphological convergence on three-dimensional shapes. Earlier versions of the two core functions of this package, rate.map and conv.map, were developed to map morphological changes onto reconstructed 3D meshes, based on landmarks positions in space (Melchionna et al. 2021; Castiglione et al. 2022).

The latest updates to rate.map and conv.map offer the possibility to map morphological changes onto real 3D surfaces provided by the user, preserving the resolution and detail of the original meshes. As a 3D landmark-based toolkit, RRmorph requires 3D coordinate data obtained from 3D surfaces to register and analyze shape variation effectively.

2. Preparing the data

To illustrate the functionalities of RRmorph, we will use a subset of the supporting data published by Melchionna et al. (2024). The dataset consists of two different sets of three-dimensional landmarks collected on several specimens belonging to 51 living Primate species, including Catarrhini, Platyrrhini and Strepsirrhini. The first set includes 18 endocast landmarks, the second includes 41 landmarks sampled on cranial surfaces.

To run the examples, the user needs to load RRmorph data. RRmorphdata object includes four different elements:

library(RRmorph)
library(Morpho)
library(Rvcg)
library(rgl)
library(Arothron)

setwd("YOUR_DIRECTORY")
da<-"https://github.com/pasraia/RRmorph_example_data/raw/refs/heads/main/RRmorphdata.rda"
download.file(url=da,destfile = paste0(getwd(),"/RRmorphdata.rda"))
load(paste0(getwd(),"/RRmorphdata.rda"))

3. The analysis of shape

Landmarks coordinates need to be analysed via usual GM procedures. This involves applying a Procrustes registration, which standardizes the coordinates by removing differences in size, position, and orientation. To reduce the dimensionality and explore the shape variation, the coordinates are elaborated by using Principal Component Analysis (PCA), or Relative Warp Analysis (RWA) if one is interested in separating the effect of the affine and non-affine variation (Schlager, 2017). The Principal Components (PCs) or the Relative Components (RWs) and their related scores serve as input data for the RRmorph functions.

There are different packages to perform the PCA or the RWA with GM data (i.e. Morpho, geomorph, shapes). In this example, PCA is performed by using the function procSym in Morpho (Schlager, 2017).

pca_endo<-procSym(endo.set)
pca_cran<-procSym(crania.set)

plot(pca_endo$PCscores[,1:2],pch=16,col=as.factor(dat.prima$group),
     main = "PC1&2 plot - endocasts",asp = 1)
text(pca_endo$PCscores[,1:2],labels = dat.prima$species,pos = 3,cex=0.6)
plot(pca_cran$PCscores[,1:2],pch=16,col=as.factor(dat.prima$group),
     main = "PC1&2 plot - crania",asp=1)
text(pca_cran$PCscores[,1:2],labels = dat.prima$species,pos = 3,cex=0.6)

In the PC1/PC2 plot of the endocast shape, Strepsirrhini are clearly separated from Haplorrhini on PC1, except for Aotus trivirgatus, which occupies an intermediate position between the two groups. Along the PC2, Homo sapiens and Alouatta stand out placing at the two extremes of the axis. The clear partition between Haplorrini and Strepsirrhini is also evident along PC1 on the cranial shape. Cynocephaline primates (i.e., Mandrillus, Papio, Theropithecus) place opposite to Strepsirrhini along PC2. Homo sapiens is positioned towards one extreme of PC1.

4. rate.map

The rate.map function is based on the evolutionary rate computations performed via RRphylo (Castiglione et al. 2018). The RRphylo method provides estimates of phenotypic rates per branch and reconstructs phenotypic values at internal nodes. It can handle both continuous and ordinal traits, as well as univariate and multivariate data (Castiglione et al. 2018, 2020).

Using vectors of Principal Component (PC) or Relative Warps (RW) scores, rate.map identifies the axes associated with the highest and lowest evolutionary rates and reconstructs the morphology weighted accordingly. This allows visualization of where and how phenotypic changes occurred most prominently between any given taxa under an evolutionary framework.

The algorithm operates by comparing either a single species/node or a pair of species/nodes. When two species or nodes are selected, they are compared to their most recent common ancestor; if a single species or node is chosen, it is compared to its immediate parental node. The rate values (one per PC/RW axis), along the evolutionary path of the selected species/nodes, are then ranked from highest to lowest.

The highest and lowest rates, along with their corresponding PC/RW axes, contribute most significantly to morphological variation from the ancestral to the descendant shapes. These axes are selected using the Extremum Distance Estimator (EDE) method (Christopoulos, 2022) implemented via the ede function (inflection package). By using the selected PCs/RWs, a three-dimensional mesh is reconstructed, starting from the consensus shape. Since RRphylo can estimate the phenotype at nodes (aces), the surface of the ancestral shape is reconstructed as well, applying the Ball Pivoting algorithm (Bernardini et al. 1999; vcgBallPivoting, Rvcg package). If no prior surface is provided, rate.map automatically performs this operation. The surface, however, can be also computed externally, using different software. In such cases, it is preferable to use the consensus shape configuration as a reference for reconstruction.

With vcgBallPivoting, the surface can be realized and visualized as follow:

mshapeE<-vcgBallPivoting(pca_endo$mshape)
mshapeC<-vcgBallPivoting(pca_cran$mshape)

mfrow3d(nr=1,nc=2,sharedMouse=TRUE)
shade3d(mshapeE,col="lightgray",specular="black")
wire3d(mshapeE,col="black",specular="black")
next3d()
shade3d(mshapeC,col="lightgray",specular="black")
wire3d(mshapeC,col="black",specular="black")

The comparison between reconstructed shapes is based on a triangle-by-triangle area difference between the triangular meshes of each species/node and its most recent common ancestor. Since positive RRphylo rates means larger coordinate values, mesh triangles mapping to morphologically expanding areas will appear larger than their respective counterparts and vice versa.

Note: the reconstructed surfaces of the compared species/nodes are intentionally altered, as the function utilizes only a subset of the morphological variation carried by the selected axes. This approach emphasizes which aspects of shape variation evolved most rapidly rather than using evolutionary rates to quantify the overall magnitude of phenotypic change.

In the following example, we first compute the evolutionary rates with RRphylo function (see RRphylo vignette for more details). This analysis relies on the phenotypic vectors retrieved earlier (i.e., pca_endo and pca_cran objects). Next, we run rate.map function to compare the endocasts of Macaca mulatta and Homo sapiens. Each species will be compared to its most recent common ancestor, generating two vectors of triangle by triangle area differences as the output. The length of the two vectors will be equal to the number of triangles composing the surfaces.

The reconstructed surfaces will be colored according to the area difference. The default color palette goes from dark red, which means area reduction, to dark blue, which conversely means area expansion. Zero values are always white.

PCscore_endo<-RRphylo::treedataMatch(tree.prima,pca_endo$PCscores)$y
RRendo<-RRphylo::RRphylo(tree.prima,PCscore_endo)

PCscore_cran<-RRphylo::treedataMatch(tree.prima,pca_cran$PCscores)$y
RRcran<-RRphylo::RRphylo(tree.prima,PCscore_cran)

rm_endo<-rate.map(x = c("Macaca_fuscata","Homo_sapiens"),
                  RR = RRendo,
                  scores = PCscore_endo,
                  pcs = pca_endo$PCs, 
                  mshape = pca_endo$mshape)

Whenever real surfaces are provided by the user, rate.map transfers the triangle by triangle area differences from the reconstructed surfaces to the three-dimensional meshes. To this aim, the landmarks/semilandmarks sets must be further provided to the function. The transfer is accomplished via interpolation. The interpolation algorithm is performed with interpolMesh, embedded in RRmorph package (see below for an extended explanation). The real surfaces, and the related sets, must be provided as named lists.

endo.list<-arraytolist(endo.set[,,c("Macaca_fuscata","Homo_sapiens")])

rm_endoS<-rate.map(x = c("Macaca_fuscata","Homo_sapiens"),
                   RR = RRendo,
                   scores = PCscore_endo,
                   pcs = pca_endo$PCs, 
                   mshape = pca_endo$mshape,
                   refsur = endo.sur,
                   refmat = endo.list)

5. conv.map

Morphological convergence can be investigated using the function conv.map. A crucial prerequisite is that the user must have prior knowledge of which species or clades exhibit convergence to avoid potentially misleading inference. We strongly recommend running the function search.conv from the RRphylo package to identify convergent taxa before proceeding. The functioning and the interpretation of results of search.conv are described in Castiglione et al. 2019 and in the function vignette.

SC<-search.conv(RR = RRcran, y = PCscore_cran) 
# Please note this result is not significant because we cut the data to reduce their size

conv.map requires phenotypic vectors (i.e., PC or RW scores) for the species under study. These vectors are computed relative to the origin of the PC axes (the consensus shape). The angle they form is a correlation coefficient. The angle ranges from 0° to 180°, according with the following: - an angle close to 0° indicates strong morphological convergence - angles around 90° suggests morphological dissimilarity - angles close to 180° indicate phenotypes evolving in opposite directions. conv.map computes angles between paired PC/RW vector by comparing two species at time. For each species pair, the function iteratively excludes one PC/RW at a time and recalculates the angle. If the angle increases after exclusion, it implies the excluded PC/RW contributes significantly to convergence, and is therefore retained. Conversely, if the angle decreases, the PC/RW is considered irrelevant to convergence and discarded.

Using the selected PCs/RWs, conv.map automatically reconstructs the three-dimensional shapes of the species pair under examination. These reconstructed surfaces are then compared to each other by calculating the triangle-by-triangle area differences, which are color-coded based on the amount of area difference. Regions without significant differences are assigned a by default a deep blue colour, shading towards white as the differences increase.

As explained in the rate.map section, if the user provides real surfaces along with their corresponding landmark sets, the algorithm will interpolate and transfer the computed area differences onto the real meshes and colour it accordingly. Below, we provide an example using real surfaces. The arguments x1 and x2 specify the clades to be examined. If convergence is being assessed within a single clade, the user must specify the species to be plotted in x1. In this example, we investigate morphological convergence between howler monkeys (genus Alouatta) and Lemuroidea (see Melchionna et al. 2024). The species names entered in x1 and x2 correspond to those whose surfaces we wish to visualize for the best convergent areas.

cran.list<-arraytolist(crania.set[,,c("Alouatta_guariba",
                                    "Alouatta_pigra",
                                    "Hapalemur_griseus",
                                    "Propithecus_verreauxi")])

cm_crann<-conv.map(x1=c("Alouatta_guariba","Alouatta_pigra"),
                   x2=c("Hapalemur_griseus","Propithecus_verreauxi"),
                   scores = PCscore_cran,pcs = pca_cran$PCs,
                   mshape = pca_cran$mshape,refmat = cran.list,refsur = crania.sur)

The resulting 3D plot will display a grid comparing two species at time, one from x1 and one from x2, displayed in the upper and the lower triangle respectively.

Besides the 3D visualization, conv.map provides additional outputs:

#>                                        real.angle selected   others  ang.diff
#> Alouatta_guariba-Propithecus_verreauxi   73.36717 10.36342 100.7186 -90.35513
#> Alouatta_pigra-Propithecus_verreauxi     67.94427 10.41381  97.0251 -86.61129
#> Alouatta_pigra-Hapalemur_griseus         80.60542 16.13478 114.4362 -98.30137
#> Alouatta_guariba-Hapalemur_griseus       84.81633 17.78638 115.5649 -97.77855
#>                                        p.value
#> Alouatta_guariba-Propithecus_verreauxi   0.002
#> Alouatta_pigra-Propithecus_verreauxi     0.003
#> Alouatta_pigra-Hapalemur_griseus         0.015
#> Alouatta_guariba-Hapalemur_griseus       0.010
#>                       Alouatta_guariba Alouatta_pigra Hapalemur_griseus
#> Alouatta_guariba            0.00000000             NA        0.06112345
#> Alouatta_pigra                      NA     0.00000000        0.05346897
#> Hapalemur_griseus           0.06112345     0.05346897        0.00000000
#> Propithecus_verreauxi       0.04398796     0.04564388                NA
#>                       Propithecus_verreauxi
#> Alouatta_guariba                 0.04398796
#> Alouatta_pigra                   0.04564388
#> Hapalemur_griseus                        NA
#> Propithecus_verreauxi            0.00000000

6. interpol.mesh

The function interpol.mesh is embedded within both rate.map and conv.map. As a stand-alone tool interpol.mesh can be used to interpolate any vector of values from a reconstructed mesh (sur argument) to a real one (refsur argument), by using a set of landmarks as target (refmat argument). The values can be linked to either the triangles or the vertices of sur, as specified by the argument element. The algorithm starts by locating a set of points (NNps) on the real surface (refsur) provided by the user. Each NNps corresponds to a single nearest neighbor for each vertex of the reconstructed surface (or barycenter if element="triangles"). Interpolation is then performed by selecting the nearest-neighbor vertices for each point on the real surface and computing the weighted mean of their values based on distance.

In the example below the values are associated to the triangles.

The first step of the interpolation is computing the barycenter for each triangle of sur (blue dots in panel A) and linking the area difference value to its coordinates. The position of each barycenter is then approximated on refsur by selecting the nearest-neighbor vertex (green dot in panel B). This operation is possible because the real surfaces are superimposed to the reconstructed meshes. Then, for each point of refsur (red dot in panel C), the algorithm selects the closest projected barycenters (green dots). The mean value of the k adjoining points, weighted by the Euclidean distances is computed as to represent the interpolated value at the focal point. Note that by default interpol.mesh select k=4 nearest neighbor points, but it can be changed by the user. interpol.mesh returns a vector of interpolated values as long as the number of triangles of the given real surfaces. Then the real surface can be colored by the function col2mesh (panel D).

References