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 270 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 156 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.
First NMDS
We recorded the touch of each plant in each point of transect. So we record not just the presence of each plant in each point of the transect proposed by Daget & Poissonet (1974). By considering the touch of each plant, we take into account rather the plant density by species. Which may be important to evaluate grazing effects (Castañeda et al. 2023)
We carried out Bray-Curtis similarity index as distance measure between plant transects.
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.1626598
Dimensions one and two of the NMDS
Figure 1: Site ordination
To remember the coding, the initial letter C,S and P represent Castril, Santiago and Pontones, respectively. The following number (01 up to 12) is the transect number in each common. The letters A, B and C are the replicates in each pastoral unit. Replicates A were carried out during the spring 2022. Replicates B and C were carried out in 2023. The last number(1,2 or 3) represents the Tranhumance modality explained below.
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
Axis 1 discriminates plant associations characterized mainly by perennial plants as Arenaria tetraqueta, Helianthemum spp, Thymus serpylloides and Teucrium aureum from plant associations characterized by annual plants as for example Vulpia spp, Xeranthemum inapertum and Bombycilaena erecta.
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
Our results reveal that axis 3 may be related to commons. As Castril differentiates from Santiago and Pontones.
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
This result suggests that axis 2 may be related to transhumance modalities. As modality 1 differentiates from modality 3.
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.72438 0.68940 0.1913 0.002 **
## H 0.97481 -0.22305 0.1523 0.006 **
## J 0.71156 -0.70263 0.1232 0.016 *
## D 0.92059 -0.39054 0.1497 0.005 **
## invsimp 0.83090 -0.55642 0.1523 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: 0x0000013cfcc10498>
## <environment: namespace:stats>Significant alpha diversity Variables
Figure 8a: 1-2 axes
Figure 8b: 1-3 axes
Second NMDS
Jaccard similarity index for binary species data i.e. presence/absence of a plant in a transect
jacc_dist <- vegdist(presence_22_23_agg_ok_2, "jaccard", binary = TRUE)
nmds_presence <- metaMDS(comm = jacc_dist, autotransform = FALSE, k=3,weakties = TRUE, model = "global", maxit = 300, try = 40,trymax = 100) # any dissimilarity index in vegdist can be used for metaMDS## [1] 0.1785058
Figure 1_bis: Site ordination
Species and sites ordination
envfitto identify species or environmental variables
which are driving the pattern
Warning in cor(Hperm, take): l’écart type est nul
Significant Species
## Warning: ggrepel: 11 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
Figure 2_bis: Species ordination
Dimensions one and three of the NMDS
Figure 3_bis: Site ordination
Figure 4_bis: Sites and Species ordination
Ellipses by Commons
1-2 axes
Figure 6a_bis: 1-2 axes
Figure 6b_bis: 1-3 axes
Ellipses by Transhumance modality
Figure 7a_bis: 1-2 axes
Figure 7b_bis: 1-3 axes
Hierarchical cluster analysis (HCA)
HCA using Bray-Curtis distance as we consider here plant density which take into account total touch of each plant in a transect.
Determining Optimal Clusters
Average Silhouette Method
Figure 9: Silhouette Method with quantitative data i.e. plant density
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
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, so instead of using 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 = 2, border = 2:5)
Figure 10: Ward’s Cluster
Regarding the coding, the initial letter C,S and P represent Castril, Santiago and Pontones, respectively. The following number (01 up to 12) is the transect number in each common. The letters A, B and C are the replicates in each pastoral unit. Replicates A were carried out during the spring 2022. Replicates B and C were carried out in 2023. The last number(1,2 or 3) represents the Tranhumance modality.
Indicator species by Ward’s Cluster
library(labdsv)
indval_3 = indval(community_data_density_22_23, cutree(hc_bray_ward, k = 2), 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, Avenula_bromoides_ok_ok, Cuscuta_sp_ok_n.i, Anthyllis_vulneraria_ok_ok, Dianthus_brachyanthus_ok_ok, Armeria_spp_ok_ok, Linaria_perennial_ok_n.i, Leucanthemopsis_pallida_ok_ok, Scabiosa_andryaefolia_ok_ok |
| 2 | Astragalus_spp_ok_ok, Eryngium_campestre_ok_ok, Koeleria_vallesiana_ok_ok, Xeranthemum_inapertum_ok_ok, Aegilops_geniculata_ok_ok, Androsace_maxima_ok_ok, Bombycilaena_erecta_ok_ok, Bromus_ruderal_ok_ok, Erodium_cicutarium_ok_ok, Ononis_spinosa_ok_ok, Poa_bulbosa_ok_ok, Taraxacum_sp_ok_n.i, Trifolium_annual_ok_n.i, Vulpia_spp_ok_ok, Medicago_annual_ok_ok, Bromus_squarrosus_ok_ok, Bufonia_tenuifolia_ok_ok, Leontodon_spp_ok_ok, Anthemis_sp_ok_n.i, Crepis_foetida_ok_ok, Helianthemum_annual_ok_n.i, Convolvulus_lineatus_ok_ok, Convolvulus_arvensis_ok_ok, Carduus_platypus_ok_ok, Paronychia_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 and the sum of indicator values (Borcard et al. 2018).
# 2 clusters
indval_2 = indval(community_data_density_22_23, cutree(hc_bray_ward, k = 2), numitr = 1000)
indval_2_cluster = factor(indval_2$maxcls[indval_2$pval <= 0.05]);addmargins(table(indval_2_cluster))## indval_2_cluster
## 1 2 Sum
## 18 24 42## [1] 20.50126The sum of indicator values is 20.5 and 42 species are highlighted as indicator species: 18 and 24 for cluster 1 and 2, respectively.
## indval_3_cluster
## 1 2 3 Sum
## 14 12 17 43## [1] 18.52722## indval_4_cluster
## 1 2 3 4 Sum
## 8 22 5 11 46## [1] 18.65864## indval_5_cluster
## 1 2 3 4 5 Sum
## 8 17 3 12 6 46## [1] 17.99044## indval_6_cluster
## 1 2 3 4 5 6 Sum
## 7 19 3 11 5 4 49## [1] 18.22482## indval_7_cluster
## 1 2 3 4 5 6 7 Sum
## 5 15 2 11 7 11 4 55## [1] 19.71672## indval_8_cluster
## 1 2 3 4 5 6 7 8 Sum
## 5 13 1 4 7 11 6 4 51## [1] 18.126437 clusters maximizes the number of clusters with significant indicator species (55) and the sum of indicator values is the second highest (19.7).
IndVal within 7 clusters
Figure 11: Ward’s method with 7 clusters
Indicator species within the seven clusters
| 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, Aegilops_geniculata_ok_ok, Androsace_maxima_ok_ok, Bombycilaena_erecta_ok_ok, Microthlaspi_perfoliatum_ok_ok, Ononis_spinosa_ok_ok, Bromus_squarrosus_ok_ok, Acinos_annual_ok_n.i, Plantago_ruderal_ok_n.i, 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, 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, Ranunculus_paludosus_ok_ok, Valeriana_tuberosa_ok_ok, Crucianella_angustifolia_ok_ok, Cerastium_gibraltericum_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 |
HCA with Jaccard similarity index
Jaccard similarity index for binary species data i.e. presence/absence of a plant in a transect.
Figure 12: Silhouette Method with presence/absence data
Figure 13: Ward’s method using Jaccard distances highlighting 2 clusters
Indicator species within the two clusters
| cluster | species |
|---|---|
| 1 | Arenaria_tetraquetra_ok_ok, Asperula_aristata_ok_ok, Centaurea_jaennensis_ok_ok, Coronilla_minima_ok_ok, Draba_hispanica_ok_ok, Ononis_pusilla_ok_ok, Poa_ligulata_ok_ok, Seseli_montanum-subsp-granatense_ok_ok, Silene_conica.colorata_ok_ok, Teucrium_aureum_ok_ok, Thymus_serpylloides_ok_ok, Silene_legionensis_ok_ok, Avenula_bromoides_ok_ok, Cuscuta_sp_ok_n.i, Anthyllis_vulneraria_ok_ok, Dianthus_brachyanthus_ok_ok, Armeria_spp_ok_ok, Linaria_sp_ok_n.i, Saxifraga_carpetana_ok_ok, Leucanthemopsis_pallida_ok_ok, Scabiosa_andryaefolia_ok_ok, Parentucellia_latifolia_ok_ok |
| 2 | Eryngium_campestre_ok_ok, Xeranthemum_inapertum_ok_ok, Aegilops_geniculata_ok_ok, Androsace_maxima_ok_ok, Bombycilaena_erecta_ok_ok, Erodium_cicutarium_ok_ok, Taraxacum_sp_ok_n.i, Vulpia_spp_ok_ok, Medicago_annual_ok_ok, Bufonia_tenuifolia_ok_ok, Crucianella_angustifolia_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, Herniaria_cinerea_ok_ok, Convolvulus_arvensis_ok_ok, Orlaya_daucoides_ok_ok, Carduus_platypus_ok_ok, Brachypodium_distachyon_ok_ok, Dactylis_glomerata_ok_ok, Valerianella_coronata_ok_ok |