Both year results : Question 1
SP
2023-10-12
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.
We merged datasets of both year. We have 336 species identified at genus level from which 243 are identified at species level.
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 aggregation we have 384 and 271 taxa respectively. Then we keep in the dataset those taxa that are identified at least at the genus level. So we remove 48 not identified plants that represent 0.4% of plant density considering all plants. Finally, we filtered out those taxa present in only one transect which was the case of 66 taxa.
After this dataset stabilisation we have 157 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.
NMDS
library(vegan)
bray_dist_log = vegdist(log1p(community_data_density_22_23), method="bray", binary=FALSE, diag=FALSE, upper=FALSE,na.rm = FALSE)
nmds_density_22_23 <- metaMDS(comm = bray_dist_log, autotransform = FALSE, k=3,weakties = TRUE, model = "global", maxit = 300, try = 40,trymax = 100)The fit of a NMDS ordination can be assessed by plotting the original dissimilarities against the (Euclidean) ordination distances (Clarke 1993).
## [1] 0.1625238
sites_density <- as.data.frame(scores(nmds_density_22_23$points))
sites_density = sites_density %>% mutate(transect = rownames(nmds_density_22_23$points)) %>% separate(transect,into = c("common","code","replicate"),sep=c(1,3)) %>% mutate(year = case_when(replicate == 'C' | replicate == 'B' ~ 23, replicate == "A" ~ 22))
sites_density <- within(sites_density, {common<-factor(common)}) ; sites_density <- within(sites_density, {year<-factor(year)})
colnames(sites_density) = c("NMDS1" , "NMDS2" , "NMDS3", "common", "code","replicate","year")
thr_cat=c(1,1,1,1,1,1,1,1,1,2,2,2,3,3,3,1,1,1,3,3,3,1,1,1,1,1,1,2,2,2,2,2,2,3,3,3,3,3,1,1,1,1,1,1,3,3,3,3,3,3,1,1,1,2,2,2,3,3,3,3,3,3,3,3,3,3,3,3,2,2,1,1,1)
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
Figure 1: Site ordination
Species and sites ordination
envfitto identify species or environmental variables
which are driving the pattern
Significant Species
Figure 2: Species ordination
Dimensions one and three of the NMDS
Figure 3: Site ordination
Figure 4: Sites and Species ordination
Ellipses by year
##
## 22 23 Sum
## 25 48 73
Figure 5: 1-2 axes
Ellipses by Commons
##
## C P S Sum
## 18 20 35 731-2 axes
Figure 6a: 1-2 axes
Figure 6b: 1-3 axes
Ellipses by Transhumance modality
Three categories :
Pastoral units used by Short-distance Transhumant Shepherds (STS)
Pastoral units used by Long-distance Transhumant Shepherds (LTS) that are nearby STS
Isolated and and distant pastoral units used by STS that are surrounded by LTS
Pastoral units in which there were a transhumant practice change in recent years (i.e. STS –>LTS or LTS –>STS )
Pastoral units used by LTS
##
## C P S Sum
## 1 12 6 12 30
## 2 3 6 5 14
## 3 3 8 18 29
## Sum 18 20 35 73
Figure 7a: 1-2 axes
Figure 7b: 1-3 axes
Alpha diversity and NMDS
S = specnumber(community_data_density_22_23_bef_agg) # Species richness
D <- diversity(community_data_density_22_23_bef_agg, "simpson")
H <- diversity(community_data_density_22_23_bef_agg,"shannon")
J <- H/log(S) # Pielou's evenness
invsimp <- diversity(community_data_density_22_23_bef_agg, "inv") # inverse Simpson
diversity_csp = data.frame(S,H,J,D,invsimp)
csp.envfit <- envfit(nmds_density_22_23, diversity_csp, permutations = 999)
head(csp.envfit)## $vectors
## NMDS1 NMDS2 r2 Pr(>r)
## S 0.71962 0.69436 0.1978 0.002 **
## H 0.97832 -0.20712 0.1525 0.005 **
## J 0.71500 -0.69913 0.1215 0.008 **
## D 0.92483 -0.38039 0.1492 0.004 **
## invsimp 0.83778 -0.54601 0.1506 0.003 **
## ---
## 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: 0x000001d82e935728>
## <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
Figure 8a: 1-2 axes
Figure 8b: 1-3 axes
Hierarchical cluster analysis
Determining Optimal Clusters
Average Silhouette Method
Figure 9: Silhouette Method
UPGMA cluster analysis
hc_bray_upgma=hclust(bray_dist_log, method = "average") # UPGMA
plot(hc_bray_upgma, ylab= "Bray-Curtis coefficient", xlab="Transects") # hang = -1
rect.hclust(hc_bray_upgma, k = 3, border = 2:5)
Figure 9: UPGMA
Indicator Species by cluster
\[a_{ij} = N individuals_{ij} /N individuals_{i}\] \[f_{ij} = N sites_{ij} /N sites_{j}\]
\[IndVal_{ij} = a_{ij} * f_{ij}\]
(Dufrêne and Legendre, 1997)
For each species i in each site group j:
\(a_{ij}\) is the relative average abundance of species i in cluster j compared to all cluster in the study. a also called specificity
\(f_{ij}\). is the relative frequency of occurrence of species i in cluster j. f is a measure of fidelity
library(labdsv)
indval_3 = indval(community_data_density_22_23, cutree(hc_bray_upgma, k = 3), numitr = 1000)
indval_3_df = indval_3$maxcls[indval_3$pval <= 0.05] %>% as.data.frame(.) ; colnames(indval_3_df) = c("cluster")
indval_3_df = indval_3_df %>% mutate(species = rownames(indval_3_df))%>% group_by(cluster) %>% dplyr::summarise_all(list(~toString(unique(.)))) ; kable(indval_3_df)| cluster | species |
|---|---|
| 1 | Merendera_montana_ok_ok, Poa_ligulata_ok_ok |
| 2 | Pilosella_pseudopilosella_ok_ok, Androsace_maxima_ok_ok, Ononis_spinosa_ok_ok, Bromus_squarrosus_ok_ok, Orlaya_daucoides_ok_ok |
| 3 | Teucrium_webbianum_ok_ok, Convolvulus_lineatus_ok_ok, Brachypodium_distachyon_ok_ok, Phlomis_lychnitis_ok_ok, Aphyllanthes_monspeliensis_ok_ok, Asphodelus_cerasiferus_ok_ok, Santolina_chamaecyparissus_ok_ok, Artemisia_campestris_ok_ok, Filago_fuscescens_ok_ok, Stipa_apertifolia_ok_ok |
Ward’s Minimum Variance Clustering
Ward algorithm requires distances to be metric (i.e. possible to display in metric Euclidean space) and Bray-Curtis is not metric, instead of raw Bray-Curtis distances we use their square-root transformation.
hc_bray_ward <- hclust(sqrt(bray_dist_log), method = "ward.D2")
plot(hc_bray_ward, ylab= "Height", xlab="Transects")
rect.hclust(hc_bray_ward, k = 3, border = 2:5)
Figure 10: Ward’s Cluster
Indicator species by Ward’s Cluster
library(labdsv)
indval_3 = indval(community_data_density_22_23, cutree(hc_bray_ward, k = 3), numitr = 1000)
indval_3_df = indval_3$maxcls[indval_3$pval <= 0.05] %>% as.data.frame(.) ; colnames(indval_3_df) = c("cluster")
indval_3_df = indval_3_df %>% mutate(species = rownames(indval_3_df))%>% group_by(cluster) %>% dplyr::summarise_all(list(~toString(unique(.)))) ; kable(indval_3_df)| cluster | species |
|---|---|
| 1 | Arenaria_tetraquetra_ok_ok, Asperula_aristata_ok_ok, Coronilla_minima_ok_ok, Draba_hispanica_ok_ok, Helianthemum_perennial_ok_ok, Ononis_pusilla_ok_ok, Seseli_montanum-subsp-granatense_ok_ok, Teucrium_aureum_ok_ok, Thymus_serpylloides_ok_ok, Cuscuta_sp_ok_n.i, Anthyllis_vulneraria_ok_ok, Armeria_spp_ok_ok, Leucanthemopsis_pallida_ok_ok, Scabiosa_andryaefolia_ok_ok |
| 2 | Eryngium_campestre_ok_ok, Koeleria_vallesiana_ok_ok, Pilosella_pseudopilosella_ok_ok, Aegilops_geniculata_ok_ok, Androsace_maxima_ok_ok, Bombycilaena_erecta_ok_ok, Microthlaspi_perfoliatum_ok_ok, Ononis_spinosa_ok_ok, Bromus_squarrosus_ok_ok, Plantago_ruderal_ok_n.i, Echinaria_capitata_ok_ok, Taeniatherum_caput-medusae_ok_ok, Santolina_chamaecyparissus_ok_ok |
| 3 | Xeranthemum_inapertum_ok_ok, Alyssum_annual_ok_n.i, Bromus_ruderal_ok_ok, Erodium_cicutarium_ok_ok, Poa_bulbosa_ok_ok, Trifolium_annual_ok_n.i, Vulpia_spp_ok_ok, Medicago_annual_ok_ok, Leontodon_spp_ok_ok, Linaria_annual_ok_n.i, Anthemis_sp_ok_n.i, Crepis_foetida_ok_ok, Helianthemum_annual_ok_n.i, Pilosella_castellana_ok_ok, Rumex_bucephalophorus_ok_ok, Dactylis_glomerata_ok_ok, Valerianella_coronata_ok_ok, Arrhenatherum_sp_ok_n.i |
Species Specificity and Fidelity Analysis
The best partition would be the one that maximizes both the proportion of clusters with significant indicator species (i) and the sum of indicator values (Borcard et al. 2018)
## indval_2_cluster
## 1 2 Sum
## 18 25 43## [1] 20.6459## indval_3_cluster
## 1 2 3 Sum
## 14 13 17 44## [1] 18.66162## indval_4_cluster
## 1 2 3 4 Sum
## 8 22 4 10 44## [1] 18.52678## indval_5_cluster
## 1 2 3 4 5 Sum
## 8 21 3 12 5 49## [1] 19.38576## indval_6_cluster
## 1 2 3 4 5 6 Sum
## 7 20 2 11 7 4 51## [1] 19.19567## indval_7_cluster
## 1 2 3 4 5 6 7 Sum
## 5 18 2 12 7 8 5 57## [1] 20.53568## indval_8_cluster
## 1 2 3 4 5 6 7 8 Sum
## 5 15 2 4 7 11 6 4 54## [1] 19.282637 clusters following this method
Figure 11: Ward’s method with 7 clusters
indval_7_df = indval_7$maxcls[indval_7$pval <= 0.05] %>% as.data.frame(.) ; colnames(indval_7_df) = c("cluster")
indval_7_df = indval_7_df %>% mutate(species = rownames(indval_7_df))%>% group_by(cluster) %>% dplyr::summarise_all(list(~toString(unique(.)))) ; kable(indval_7_df)| cluster | species |
|---|---|
| 1 | Arenaria_tetraquetra_ok_ok, Asperula_aristata_ok_ok, Helianthemum_perennial_ok_ok, Jurinea_humilis_ok_ok, Anthyllis_vulneraria_ok_ok |
| 2 | Pilosella_pseudopilosella_ok_ok, Xeranthemum_inapertum_ok_ok, Achillea_odorata_ok_ok, Androsace_maxima_ok_ok, Bombycilaena_erecta_ok_ok, Hornungia_petraea_ok_ok, Microthlaspi_perfoliatum_ok_ok, Ononis_spinosa_ok_ok, Bromus_squarrosus_ok_ok, Festuca_gr.rubra_ok_n.i, Orlaya_daucoides_ok_ok, Plantago_ruderal_ok_n.i, Echinaria_capitata_ok_ok, Sanguisorba_sp_ok_n.i, Crataegus_sp_ok_n.i, Taeniatherum_caput-medusae_ok_ok, Rosa_sp_ok_n.i, Centaurea_castellanoides_ok_ok |
| 3 | Koeleria_vallesiana_ok_ok, Santolina_chamaecyparissus_ok_ok |
| 4 | Poa_ligulata_ok_ok, Aegilops_geniculata_ok_ok, Erodium_cicutarium_ok_ok, Medicago_annual_ok_ok, Leontodon_spp_ok_ok, Anthemis_sp_ok_n.i, Crepis_foetida_ok_ok, Geranium_molle_ok_ok, Helianthemum_annual_ok_n.i, Herniaria_cinerea_ok_ok, Pilosella_castellana_ok_ok, Rumex_bucephalophorus_ok_ok |
| 5 | Draba_hispanica_ok_ok, Ononis_pusilla_ok_ok, Teucrium_aureum_ok_ok, Thymus_serpylloides_ok_ok, Arenaria_leptoclados.serpyllifolia_ok_n.i, Cuscuta_sp_ok_n.i, Armeria_spp_ok_ok |
| 6 | Centaurea_jaennensis_ok_ok, Coronilla_minima_ok_ok, Merendera_montana_ok_ok, Silene_conica.colorata_ok_ok, Valeriana_tuberosa_ok_ok, Leucanthemopsis_pallida_ok_ok, Scabiosa_andryaefolia_ok_ok, Spergula_pentandra_ok_ok |
| 7 | Bromus_ruderal_ok_ok, Poa_bulbosa_ok_ok, Carduus_platypus_ok_ok, Valerianella_coronata_ok_ok, Scandix_spp_ok_ok |