Characterizing plant assemblies

Identify plant assemblages highlighting dominant and indicator species of the pastoral mountains in Castril, Santiago and Pontones (CSP) pasturelands.

We hypothesise that pastoral commons shape plant diversity in CSP. Second, we hypothesise that livestock mobility, referred here to long- or short-transhumance, shape vegetation in CSP.

Species of the same genus having similar ecology have been grouped. This is the case of Helianthemum_perennial composed of Helianthemum_oelandicum, Helianthemum_cinereum and Helianthemum_appenninum. In this group, we did not include Helianthemum_salicifolium as it is an annual plant.

Also some species presenting difficulties of identification have been grouped within their genus. This is the case of Anthemis_arvensis and Anthemis_nobilis.

Before and after the species aggregation we have 305 and 241 taxa in the dataset.

Transect repartition

Assessing patterns in the ‘arrangement’ of plant transects by employing ordination methods in order to evaluate the effect of transhumance and commons on plant diversity.

(For the moment we just assess commons differences)

NMDS

library(vegan)
nmds_density_23 <- metaMDS(comm = community_data_density_agg, distance = "bray",  k=4, trymax=1000) 
print(nmds_density_23) # round(nmds_density_23$stress, digits = 3)
## 
## Call:
## metaMDS(comm = community_data_density_agg, distance = "bray",      k = 4, trymax = 1000) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     wisconsin(sqrt(community_data_density_agg)) 
## Distance: bray 
## 
## Dimensions: 4 
## Stress:     0.1275236 
## Stress type 1, weak ties
## Best solution was repeated 1 time in 20 tries
## The best solution was from try 13 (random start)
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'wisconsin(sqrt(community_data_density_agg))'
sites_density  <- as.data.frame(scores(nmds_density_23$points))
sites_density = sites_density  %>% mutate(site_23)%>% separate(site_23,into = c("common","code"),sep=1) 
sites_density <- within(sites_density, {common<-factor(common)})
colnames(sites_density) = c("NMDS1"  , "NMDS2" ,  "NMDS3"  , "NMDS4" ,  "common", "code")
thr_cat=c(1,1,1,2,3,1,1,1,1,2,3,1,1,1,3,3,1,2,3,3,3,3,2,1,1,1,3,3,1,2,3,3,3,3,1,3,1,1,2,2,3,3,3,1,1,2,2,3)
sites_density = sites_density %>% mutate(thr_cat = thr_cat)
sites_density <- within(sites_density, {thr_cat<-factor(thr_cat)})

Dimensions one and two of the NMDS

(mds_plot_1 = ggplot(sites_density, aes(x = NMDS1, y = NMDS2, label = rownames(sites_density),color=common)) +  geom_text_repel(alpha = 1, size = 2) +
  scale_colour_manual(values=c("blue","green3","red"))+
   scale_x_continuous(limits = c(-0.8,0.8), breaks = c(-0.8,-0.4,0.0,0.4,0.8)) + 
scale_y_continuous(limits = c(-0.8,0.8), breaks = c(-0.8,-0.4,0.0,0.4,0.8)) + theme_minimal() +  theme(legend.position="bottom"))
Figure 1: Site ordination

Figure 1: Site ordination

Species and sites ordination

envfitto identify species or environmental variables which are driving the pattern

csp.spp.fit <- envfit(nmds_density_23, community_data_density_agg, permutations = 999)
species_scores_csp <- as.data.frame(scores(csp.spp.fit, display = "vectors")) * ordiArrowMul(csp.spp.fit)
species_scores_csp <- cbind(species_scores_csp, species = rownames(species_scores_csp))
species_scores_csp <- cbind(species_scores_csp, pval = csp.spp.fit$vectors$pvals)
sig.species_scores_csp <- subset(species_scores_csp, pval<=0.05)
sig.species_scores_csp = sig.species_scores_csp%>% separate(species,into = c("genus","species","id_genus","id_sp"),remove = FALSE,sep ='_') %>%
  tidyr::unite(species,genus, species,sep = "_") %>% filter(id_genus == "ok")

Significant Species

mds_plot_2 = ggplot(sig.species_scores_csp, aes(NMDS1, NMDS2)) +
geom_segment(aes(x = 0, xend = NMDS1, y = 0,yend = NMDS2),
arrow = arrow(length = unit(0.1,"cm")), alpha = 0.5, colour = "grey40", lwd=0.2) +
geom_text_repel(data = sig.species_scores_csp, aes(x=NMDS1, y=NMDS2, label = species), cex = 2.5, direction = "both", segment.size = 0.25,alpha = 0.7) +
scale_x_continuous(limits = c(-0.8,0.8), breaks = c(-0.8,-0.4,0.0,0.4,0.8)) + 
scale_y_continuous(limits = c(-0.8,0.8), breaks = c(-0.8,-0.4,0.0,0.4,0.8)) +
theme_minimal()

mds_plot_2 + geom_point(data = sites_density, aes(x = NMDS1, y = NMDS2, color = common))+
  scale_colour_manual(values=c("blue","green3","red")) +
geom_text_repel(data = sites_density, aes(x=NMDS1, y=NMDS2, label = rownames(sites_density)), cex = 2.5, direction = "both", segment.size = 0.25) + theme_minimal() +  theme(legend.position="bottom") 
## Warning: ggrepel: 17 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
Figure 2: Species ordination

Figure 2: Species ordination

Dimensions one and three of the NMDS

(mds_plot_1_b = ggplot(sites_density, aes(x = NMDS1, y = NMDS3, label = rownames(sites_density),color=common)) +  geom_text_repel(alpha = 1, size = 2) +
  scale_colour_manual(values=c("blue","green3","red")) +
  scale_x_continuous(limits = c(-0.8,0.8), breaks = c(-0.8,-0.4,0.0,0.4,0.8)) + 
 scale_y_continuous(limits = c(-1,1), breaks = c(-0.9,-0.4,0.0,0.4,0.9)) + theme_minimal() +  theme(legend.position="bottom"))
Figure 3: Site ordination

Figure 3: Site ordination 

csp.spp.fit_b <- envfit(nmds_density_23, community_data_density_agg,choices=c(1,3), permutations = 999)
species_scores_csp_b <- as.data.frame(scores(csp.spp.fit_b, display = "vectors"))* ordiArrowMul(csp.spp.fit_b)
species_scores_csp_b <- cbind(species_scores_csp_b, species = rownames(species_scores_csp_b))
species_scores_csp_b <- cbind(species_scores_csp_b, pval = csp.spp.fit$vectors$pvals)
sig.species_scores_csp_b <- subset(species_scores_csp_b, pval<=0.05)
sig.species_scores_csp_b = sig.species_scores_csp_b%>% separate(species,into = c("genus","species","id_genus","id_sp"),remove = FALSE,sep ='_') %>%
  tidyr::unite(species,genus, species,sep = "_") %>% filter(id_genus == "ok")
mds_plot_2_b = ggplot(sig.species_scores_csp_b, aes(NMDS1, NMDS3)) +
geom_segment(aes(x = 0, xend = NMDS1, y = 0,yend = NMDS3),
arrow = arrow(length = unit(0.1,"cm")),alpha = 0.5, colour = "grey40", lwd=0.2) +
geom_text_repel(data = sig.species_scores_csp_b, aes(x=NMDS1, y=NMDS3, label = species), cex = 2.5, direction = "both", segment.size = 0.25,alpha = 0.7) +
scale_x_continuous(limits = c(-0.8,0.8), breaks = c(-0.8,-0.4,0.0,0.4,0.8)) + 
 scale_y_continuous(limits = c(-1,1), breaks = c(-0.9,-0.4,0.0,0.4,0.9)) +
theme_minimal()

mds_plot_2_b + geom_point(data = sites_density, aes(x = NMDS1, y = NMDS3, color = common))+ scale_colour_manual(values=c("blue","green3","red")) +
geom_text_repel(data = sites_density, aes(x=NMDS1, y=NMDS3, label = rownames(sites_density)), cex = 2.5, direction = "both", segment.size = 0.25) + theme_minimal() + theme(legend.position="bottom") 
## Warning: ggrepel: 4 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
Figure 4: Sites and Species ordination

Figure 4: Sites and Species ordination

Ellipses by Commons

1-2 axes

# function for ellipsess 
veganCovEllipse <- function (cov, center = c(0, 0), scale = 1, npoints = 100) 
{  theta <- (0:npoints) * 2 * pi/npoints
  Circle <- cbind(cos(theta), sin(theta))
  t(center + scale * t(Circle %*% chol(cov))) }
#data for ellipse, in this case using the common factor
df_ell.csp.common <- data.frame() #sets up a data frame before running the function.
for(g in levels(sites_density$common)){
  df_ell.csp.common <- rbind(df_ell.csp.common, cbind(as.data.frame(with(sites_density [sites_density$common==g,],
veganCovEllipse(cov.wt(cbind(NMDS1,NMDS2),
wt=rep(1/length(NMDS1),length(NMDS1)))$cov,center=c(mean(NMDS1),mean(NMDS2))))),common=g))}
# data for labelling the ellipse
NMDS.mean.csp=aggregate(sites_density[ ,c("NMDS1", "NMDS2")],list(common = sites_density$common), mean)
#plot 3 ellipses
ggplot(sites_density, aes(x = NMDS1, y = NMDS2,color))  + geom_point(aes(color = common)) + geom_path(data=df_ell.csp.common, aes(x=NMDS1, y=NMDS2,colour=common), linewidth=1, linetype=2) + annotate("text",x=NMDS.mean.csp$NMDS1,y=NMDS.mean.csp$NMDS2,label=NMDS.mean.csp$common) + scale_colour_manual(values=c("blue","green3","red"))+
  scale_x_continuous(limits = c(-0.8,0.8), breaks = c(-0.8,-0.4,0.0,0.4,0.8)) + 
scale_y_continuous(limits = c(-0.8,0.8), breaks = c(-0.8,-0.4,0.0,0.4,0.8)) + theme_minimal() +  theme(legend.position="bottom")
Figure 5: 1-2 axes

Figure 5: 1-2 axes 

df_ell.csp.common <- data.frame() 
for(g in levels(sites_density$common)){
  df_ell.csp.common <- rbind(df_ell.csp.common, cbind(as.data.frame(with(sites_density [sites_density$common==g,],
veganCovEllipse(cov.wt(cbind(NMDS1,NMDS3),
wt=rep(1/length(NMDS1),length(NMDS1)))$cov,center=c(mean(NMDS1),mean(NMDS3))))),common=g))}
# data for labelling the ellipse
NMDS.mean.csp=aggregate(sites_density[ ,c("NMDS1", "NMDS3")],list(common = sites_density$common), mean)
#plot 3 ellipses
ggplot(sites_density, aes(x = NMDS1, y = NMDS3,color))  + geom_point(aes(color = common)) + geom_path(data=df_ell.csp.common, aes(x=NMDS1, y=NMDS3,colour=common), size=1, linetype=2) + annotate("text",x=NMDS.mean.csp$NMDS1,y=NMDS.mean.csp$NMDS3,label=NMDS.mean.csp$common)+ scale_colour_manual(values=c("blue","green3","red"))+scale_x_continuous(limits = c(-0.8,0.8), breaks = c(-0.8,-0.4,0.0,0.4,0.8)) + 
 scale_y_continuous(limits = c(-1,1), breaks = c(-0.9,-0.4,0.0,0.4,0.9)) + theme_minimal() +  theme(legend.position="bottom") 
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Figure 6: 1-3 axes

Figure 6: 1-3 axes

Ellipses by Transhumance modality

df_ell.csp.thr <- data.frame() #sets up a data frame before running the function.
for(g in levels(sites_density$thr_cat)){
  df_ell.csp.thr <- rbind(df_ell.csp.thr, cbind(as.data.frame(with(sites_density [sites_density$thr_cat==g,], veganCovEllipse(cov.wt(cbind(NMDS1,NMDS2),wt=rep(1/length(NMDS1),length(NMDS1)))$cov,center=c(mean(NMDS1),mean(NMDS2))))),thr_cat=g))}
# data for labelling the ellipse
NMDS.mean.thr=aggregate(sites_density[ ,c("NMDS1", "NMDS2")],list(thr_cat = sites_density$thr_cat), mean)
#plot 
ggplot(sites_density, aes(x = NMDS1, y = NMDS2,color))  + geom_point(aes(color = thr_cat)) + geom_path(data=df_ell.csp.thr, aes(x=NMDS1, y=NMDS2,colour=thr_cat), linewidth=1, linetype=2) +   annotate("text",x=NMDS.mean.thr$NMDS1,y=NMDS.mean.thr$NMDS2,label=NMDS.mean.thr$thr_cat) + scale_colour_manual(values=c("blue","green3","red"))+
  scale_x_continuous(limits = c(-0.8,0.8), breaks = c(-0.8,-0.4,0.0,0.4,0.8)) + 
  scale_y_continuous(limits = c(-0.8,0.8), breaks = c(-0.8,-0.4,0.0,0.4,0.8)) + theme_minimal() +  theme(legend.position="bottom")
Figure 7a: 1-2 axes

Figure 7a: 1-2 axes 

df_ell.csp.thr_2 <- data.frame() #sets up a data frame before running the function.
for(g in levels(sites_density$thr_cat)){
  df_ell.csp.thr_2 <- rbind(df_ell.csp.thr_2, cbind(as.data.frame(with(sites_density [sites_density$thr_cat==g,], veganCovEllipse(cov.wt(cbind(NMDS1,NMDS3),wt=rep(1/length(NMDS1),length(NMDS1)))$cov,center=c(mean(NMDS1),mean(NMDS3))))),thr_cat=g))}
# data for labelling the ellipse
NMDS.mean.thr=aggregate(sites_density[ ,c("NMDS1", "NMDS3")],list(thr_cat = sites_density$thr_cat), mean)
#plot 
ggplot(sites_density, aes(x = NMDS1, y = NMDS3,color))  + geom_point(aes(color = thr_cat)) + geom_path(data=df_ell.csp.thr_2, aes(x=NMDS1, y=NMDS3,colour=thr_cat), linewidth=1, linetype=2) +   annotate("text",x=NMDS.mean.thr$NMDS1,y=NMDS.mean.thr$NMDS3,label=NMDS.mean.thr$thr_cat) + scale_colour_manual(values=c("blue","green3","red"))+
   theme_minimal() +  theme(legend.position="bottom")
Figure 7b: 1-3 axes

Figure 7b: 1-3 axes

Alpha diversity and NMDS

S = specnumber(community_data_density) # Species richness
D <- diversity(community_data_density, "simpson")
H <- diversity(community_data_density,"shannon")
J <- H/log(S) # Pielou's evenness
invsimp <- diversity(community_data_density, "inv") # inverse Simpson
diversity_csp = data.frame(S,H,J,D,invsimp)

csp.envfit <- envfit(nmds_density_23, diversity_csp, permutations = 999)
head(csp.envfit)
## $vectors
##            NMDS1    NMDS2     r2 Pr(>r)    
## S       -0.36814 -0.92977 0.4300  0.001 ***
## H       -0.82317 -0.56779 0.2237  0.005 ** 
## J       -0.96009  0.27969 0.1169  0.058 .  
## D       -0.94063 -0.33944 0.1598  0.018 *  
## invsimp -0.95742 -0.28869 0.1211  0.047 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
## 
## $factors
## NULL
## 
## $na.action
## function (object, ...) 
## UseMethod("na.action")
## <bytecode: 0x0000020b3658f838>
## <environment: namespace:stats>
csp.envfit.score <- as.data.frame(scores(csp.envfit, display = "vectors"))# * ordiArrowMul(csp.envfit) 
csp.envfit.score <- cbind(csp.envfit.score, env.variables = rownames(csp.envfit.score))
csp.envfit.score <- cbind(csp.envfit.score, pval = csp.envfit$vectors$pvals) 
sig.env.scrs <- subset(csp.envfit.score, pval<=0.05)

Significant alpha diversity Variables

mds_plot_4 = ggplot(sig.env.scrs, aes(NMDS1, NMDS2)) +
geom_segment(aes(x = 0, xend = NMDS1, y = 0,yend = NMDS2), arrow = arrow(length = unit(0.2,"cm")),alpha = 0.5, colour = "grey20", lwd=0.3) +
geom_text_repel(data = sig.env.scrs, aes(x=NMDS1, y=NMDS2, label = env.variables), cex = 5, direction = "both", segment.size = 0.4,alpha = 0.7) +
#  scale_x_continuous(limits = c(-0.9,0.9), breaks = c(-0.9,-0.6,-0.3,0.0,0.3,0.6,0.9)) + #  scale_y_continuous(limits = c(-0.9,0.9), breaks = c(-0.9,-0.6,-0.3,0.0,0.3,0.6,0.9)) +
  theme_minimal()  

mds_plot_4 + geom_point(data = sites_density, aes(x = NMDS1, y = NMDS2, color = common)) + scale_colour_manual(values=c("blue","green3","red"))+ geom_text_repel(data = sites_density, aes(x=NMDS1, y=NMDS2, label = rownames(sites_density)), cex = 2.5, direction = "both", segment.size = 0.25) + theme_minimal() + theme(legend.position="bottom")
Figure 8a: 1-2 axes

Figure 8a: 1-2 axes 

csp.envfit_b <- envfit(nmds_density_23, diversity_csp,choices=c(1,3), permutations = 999)
csp.envfit.score_b <- as.data.frame(scores(csp.envfit_b, display = "vectors")) #* ordiArrowMul(csp.envfit_b) 
csp.envfit.score_b <- cbind(csp.envfit.score_b, env.variables = rownames(csp.envfit.score_b))
csp.envfit.score_b <- cbind(csp.envfit.score_b, pval = csp.envfit$vectors$pvals) 
sig.env.scrs_b <- subset(csp.envfit.score_b, pval<=0.05)

mds_plot_4_b = ggplot(sig.env.scrs_b, aes(NMDS1, NMDS3)) +
geom_segment(aes(x = 0, xend = NMDS1, y = 0,yend = NMDS3), arrow = arrow(length = unit(0.2,"cm")), alpha =0.5,colour = "grey20", lwd=0.3) +
geom_text_repel(data = sig.env.scrs_b, aes(x=NMDS1, y=NMDS3, label = env.variables), cex = 5, direction = "both", segment.size = 0.4,alpha = 0.7) + theme_minimal()  

mds_plot_4_b + geom_point(data = sites_density, aes(x = NMDS1, y = NMDS3, color = common)) + scale_colour_manual(values=c("blue","green3","red"))+ geom_text_repel(data = sites_density, aes(x=NMDS1, y=NMDS3, label = rownames(sites_density)), cex = 2.5, direction = "both", segment.size = 0.25) + theme_minimal() + theme(legend.position="bottom")
Figure 8b: 1-3 axes

Figure 8b: 1-3 axes

Transect aggregated by Pastoral Unit

df_ell.csp.thr <- data.frame() #sets up a data frame before running the function.
for(g in levels(sites_density_b$thr_cat)){
  df_ell.csp.thr <- rbind(df_ell.csp.thr, cbind(as.data.frame(with(sites_density_b [sites_density_b$thr_cat==g,], veganCovEllipse(cov.wt(cbind(NMDS1,NMDS2),wt=rep(1/length(NMDS1),length(NMDS1)))$cov,center=c(mean(NMDS1),mean(NMDS2))))),thr_cat=g))}
# data for labelling the ellipse
NMDS.mean.thr=aggregate(sites_density_b[ ,c("NMDS1", "NMDS2")],list(thr_cat = sites_density_b$thr_cat), mean)
#plot 3 ellipses
ggplot(sites_density_b, aes(x = NMDS1, y = NMDS2,color))  + geom_point(aes(color = thr_cat)) + geom_path(data=df_ell.csp.thr, aes(x=NMDS1, y=NMDS2,colour=thr_cat), linewidth=1, linetype=2) +   annotate("text",x=NMDS.mean.thr$NMDS1,y=NMDS.mean.thr$NMDS2,label=NMDS.mean.thr$thr_cat) + scale_colour_manual(values=c("blue","green3","red"))+
  theme_minimal() +  theme(legend.position="bottom")
Figure 9: 1-2 axes

Figure 9: 1-2 axes

Species driving the pattern

mds_plot_2_c = ggplot(sig.species_scores_csp, aes(NMDS1, NMDS2)) +
geom_segment(aes(x = 0, xend = NMDS1, y = 0,yend = NMDS2),
arrow = arrow(length = unit(0.1,"cm")), alpha = 0.5, colour = "grey40", lwd=0.2) +
geom_text_repel(data = sig.species_scores_csp, aes(x=NMDS1, y=NMDS2, label = species), cex = 2.5, direction = "both", segment.size = 0.25,alpha = 0.7) + theme_minimal()

mds_plot_2_c + geom_point(data = sites_density_b, aes(x = NMDS1, y = NMDS2, color = thr_cat))+
  scale_colour_manual(values=c("blue","green3","red")) +
geom_text_repel(data = sites_density_b, aes(x=NMDS1, y=NMDS2, label = rownames(sites_density_b)), cex = 2.5, direction = "both", segment.size = 0.25) + theme_minimal() +  theme(legend.position="bottom")
Figure 10: 1-2 axes

Figure 10: 1-2 axes