##### Morphocode used to obtain the deformation grids by genetic group of figure 1 and to generate the graph of Figure 3 of the manuscript. 
Note: This script was developed and kindly provided by one of the reviewers, thank you so much by your support ######
setwd("")
getwd()
library(geomorph)
#Loading the data
ldmk<-readland.tps("Landmarks-concha.txt", specID="imageID")
#First a generalized Procrustes analysis
gpa.lands<-gpa.lands <- gpagen(ldmk)
#Obtain the Procrustes coordinates for each one of the individuals
procdata$coords
#Obtain the centroid size
procdata$Csize
#To obtain the PCA graph with the deformation grids 
pca.lands <- plotTangentSpace(gpa.lands$coords, label=TRUE)
#To assess the contribution of the components
pca.lands$pc.summary
# To plot the variation of the shape without the deformation grids 
xlab<-paste("PC1 ", "(", round(pca.lands$pc.summary$importance[2,1]*100,1),"% explained variation)", sep = "")
ylab<-paste("PC2 ", "(", round(pca.lands$pc.summary$importance[2,2]*100,1),"% explained variation)", sep = "")
plot(pca.lands$pc.scores[,1],pca.lands$pc.scores[,2],xlab = xlab, ylab = ylab)
#To obtain the eigenvalues of the eigenvectors
prinscores<-pca.lands$pc.scores
# To obtain the consensus configuration of the complete dataset (129 individuals) 
consensus <- gpa.lands$consensus
#Loading link matrix to designates the two landmarks to be connected by that link
links<-as.matrix(read.csv("snail links.csv"))
# To obtain the deformation plate of the consensus shape 
plotRefToTarget(consensus,gpa.lands$coords[,,1], links = links)
# To obtain the distribution of the variation of each landmark relative to its centroid 
plot(consensus,asp=1, type="n")
for(i in 1:dim(procdata$coords)[3]) points(procdata$coords[,,i])
plotAllSpecimens(procdata$coords)
#To obtain the previous graph but changing the color and size of the centroid of each landmark 

for(i in 1:length(gpa.lands$coords[,,3]))
  
  points(gpa.lands$coords[,,i], pch=20)

points(consensus, col="Red", cex=2, pch=20)
##############################################################
############## SECOND PART OF THE CODE#########################
##############################################################
#Load the data : 1) A file including the landmarks and 2) a file  (*.csv) including to which genetic group belongs each one of the shells. After it, performs a PGA and define the parameters of the graph 		
ldmk <-readland.tps("Landmarks-concha.txt", specID="imageID")
site<-read.csv("snail site.csv", header=T, row.names=1)
Y.gpa<-gpagen(ldmk) 
gp <- as.factor(paste(site$Group))
col.gp <- rainbow(length(levels(gp)))
names(col.gp) <- levels(gp)
col.gp <- col.gp[match(gp, names(col.gp))]
PCA <- plotTangentSpace(Y.gpa$coords, verbose = T)
xlab <- paste("Principal Component 1", "(", round(PCA$pc.summary$importance[2,1]*100, 1), "%)", sep="")
ylab <- paste("Principal Component 2", "(", round(PCA$pc.summary$importance[2,2]*100, 1), "%)", sep="")
mat <- matrix(c(4,5,0,1,1,2,1,1,3), 3)
layout(mat, widths=c(1,1,1), heights=c(1,1,0.6))
par(mar=c(4, 4, 1, 1))
plot(PCA$pc.scores[,1], PCA$pc.scores[,2], pch=21, cex=2, bg=col.gp, xlab=xlab, ylab=ylab, asp=T)
legend(-0.09, 0.07, legend= unique(gp), pch=19,  col=unique(col.gp))
ref <- mshape(Y.gpa$coords)
par(mar = c(0,0,0,0))
plotRefToTarget(ref,PCA$pc.shapes$PC1min)
plotRefToTarget(ref,PCA$pc.shapes$PC1max)
plotRefToTarget(ref,PCA$pc.shapes$PC2min)
plotRefToTarget(ref,PCA$pc.shapes$PC2max)
##############################################################################################
######################## THIRD PART OF THE CODE: PHYLOMORPHOSPACE ############################
##############################################################################################
setwd()
getwd()
library("geomorph")
library("rgl")
library("ape")
library("phytools")
#Load the tree obtained from te pairwise Fst matrix
tree<-read.tree("njfst.TREE")
d<-multi2di(tree)
#Load the consensus shape of each one of the genetic groups (*.csv)
base<-read.csv("shellmorpho.csv",row.names=1)
#To obtain the graph
plotGMPhyloMorphoSpace(d, base)