In emf-creaf/ecotraj: Ecological Trajectory Analysis

knitr::opts_chunk$set(echo = TRUE)

1. About this vignette

This vignette illustrates the ecological applications included in Sturbois et al. (2023) for the presentation of the Ecological Quality Assessment (EQA) framework. The vignette illustrates how to conduct an Ecological Quality Assessment through state-based and trajectory-based variants, and how to display the results graphically. The two data sets used in this vignette have been included in the ecotraj package for reproducibility.

2. Loading libraries

First of all, we load the libraries required for the examples, including ecotraj:

library(ecotraj)
library(vegan)
library(vegclust)
library(ape)
library(ggplot2)
library(reshape2)

3. Heathland conservation status in Landes et Marais de Glomel

3.1 Context

The nature reserve of Landes et Marais de Glomel (Brittany, France) is composed of temperate Atlantic wet heaths whose reference state is commonly considered dominated by plant communities associated to acid, nutrient poor soils that are at least seasonally water logged and dominated by Erica tetralix and E. ciliaris. This habitat is considered of community interest as part of the European directive on the conservation of natural habitats and of wild fauna and flora (97/62/CEE). For the definition of the reference envelope, surveys integrated the variability of reference states from youngest to early senescent stages, based on expert assessment with respect to the European Directive (97/62/CEE) requirements.

3.2 Loading data

The data set consists of 23 rows (stations) and 46 columns (species). The percent cover values derived from Braun-Blanquet ordinal scale were estimated for 43 species of vascular plants. The first five stations (rows) were used to define the reference envelope, and the next 18 stations (rows) where those for which the conservation status was to be assessed.

data(glomel)
#Showing a first lines /columns subset
head(glomel[,1:8])

3.3 Ecological Quality Assessment

We first create a compositional data matrix, by extracting species data from the initial data table:

glomel_comp <- as.matrix(glomel[,!(names(glomel) %in% c("ID", "Ref"))])
rownames(glomel_comp) <- glomel$ID 
dim(glomel_comp)

Using the compositional data, we can use function vegdist() from package vegan to calculate of Bray Curtis distances between ecosystem states:

glomel_bc <- vegan::vegdist(glomel_comp, method = "bray")

The distance matrix glomel_bc defines the multivariate space that represents the resemblance between ecological states.

To conduct EQA, we also need to define the conservation targets. Column Ref in glomel identifies the stations considered as reference. We define the state reference envelope using reference station IDs as follows:

glomel_env <- glomel$ID[glomel$Ref]

We are now ready to compare the assessed stations with respect to the state reference envelope using function compareToStateEnvelope():

glomel_assess <- compareToStateEnvelope(glomel_bc, glomel_env, m=1.7, distances_to_envelope = TRUE)
head(glomel_assess)

Function compareToStateEnvelope() returns SquaredDist, the squared distance to the envelope's centroid. The squared distances to the envelope averaged over the reference points defines the variability of the envelope:

mean(glomel_assess$SquaredDist[glomel_assess$Envelope])

This value should be equal to a call to function stateEnvelopeVariability() :

stateEnvelopeVariability(glomel_bc, glomel_env)

Function compareToStateEnvelope() also returns the Q statistic, which defines the quality of the tested stations according to their resemblance to the set of stations used as reference. Q is comprised between 0 (lowest quality) and 1 (highest quality). We can transform Q values into a qualitative assessment using:

glomel_assess$Status<-c(ifelse(glomel_assess$Q>=0.5,"Inside", "Outside"))

And the assessment for the 18 tested stations is:

glomel_assess[!glomel_assess$Envelope,]

3.4 Representing EQA results graphically

To represent the results in a graphical way, we start by conducting a Principal Coordinates Analysis (PCoA) using package ape:

pcoa_glom<-ape::pcoa(glomel_bc)

We can use function biplot() from package ape to show the ordination of stations with the species projected as arrows:

biplot(pcoa_glom, glomel_comp)

To generate a different representation, we first get coordinates of stations in the two first dimensions of the PCoA and copy them into the assessment data frame:

PCOA_DIM1_2 <- pcoa_glom[["vectors"]][,1:2]
glomel_assess$Dim1<-PCOA_DIM1_2[,1]
glomel_assess$Dim2<-PCOA_DIM1_2[,2]

head(glomel_assess)

We then use package ggplot to achieve a nicer representation of the assessment. Symbol colors are used to illustrate stations that do achieve conservation objectives (in green) from those that do not achieve them (in red). The size of dots represents the squared distance to the centroid of the state reference envelope:

p<-ggplot(glomel_assess,
          mapping=aes(x=Dim1,y=Dim2,size=SquaredDist, color=Status, shape=Status))+
  geom_point()+
  scale_colour_manual(values=c("#00BFC4", "#F8766D"))+
  geom_text(glomel_assess, 
            mapping=aes(x=Dim1,y=Dim2,label=Observation),
            hjust=-0.3, vjust=-0.3,size=2.5, color="Black")+
  xlab(expression("PCoA1 (27.7%)"))+
  ylab(expression("PCoA2 (13.7%)"))

p<-p + theme_minimal()+
  theme(
    # Hide some graphical elements
    panel.border = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.line = element_line(colour = "black"))

p

4. Impact of fishing activities on marine habitats

4.1 Context

We used an experimental data set built by Tauran et al. (2020) to study the impact of fishing dredges and varying fishing pressures on maerl beds, in the bay of Brest (Brittany, France). The study follows a Before-After-Control-Impact design (Stewart-Oaten et al., 1986). Briefly, three control stations were surveyed and compared to five treatment stations characterized by different fishing dredges and pressure levels (Tauran et al., 2020). Three fishing dredges are compared: (1) a clam dredge (CD), 70 to 90 kg, 1.5 m wide, 40 teeth of 11 cm each; (2) a queen scallop dredge (QSD), 120 kg,1.8 m wide, with a blade; and (3) a king scallop dredge (KSD), 190 kg, 1.8 m wide, 18 teeth of 10 cm each every 9 cm. Pressure levels were measured as the number of dredge tows performed on the zone during the experimental dredging session: 0 (i.e. control), 10, or 30 dredge tows. Samples were collected from April 2016 to April2017. Session 1 was sampled just before the experimental dredging (t0); session 2 at t0 + 1 week; session 3 at t0 + 1 month and session 4 at t0 + 12 months. Nine replicates were sampled for all treatments and sessions with a Smith-McIntyre grab (0.1m2). Abundance data were pooled at the treatment levels (i.e. CTRL1, CTRL2, CTRL3, KSD_10, CD_10, CD_30, QSD_10 and QSD_30; Appendix_D).

4.2 Loading data

The data set has 32 observations (i.e. treatments*surveys) and 250 species. We begin by loading the data:

data(glenan)
#Showing a first lines /columns subset
head(glenan[,1:6])

4.3 Definition of the ressemblance space

We separate the columns that identify the treatment and surveys.

#vector Treatment
Treatment<-glenan$Treatment

#vector Surveys
Surveys<-glenan$Surveys

We also extract the compositional data and log-transform it as initially performed in Tauran et al. (2020).

#log transformation
grab_comp<-log1p(glenan[,-c(251:252)])

Using the compositional data, we can use function vegdist() from package vegan to calculate of Bray Curtis distances between ecosystem states:

grab_bc <- vegan::vegdist(grab_comp, method = "bray")

Like in the first example, the distance matrix grab_bc defines the multivariate space that represents the resemblance between ecological states of maerl beds.

4.4 State-based EQA

In addition to the resemblance space, to conduct the EQA we need to define the conservation targets. In this case, we define the states conforming our state reference envelope manually among the observation IDs:

#Definition of the state reference envelope
ID<-rownames(grab_comp)
grab_env <- ID[c(1,9:20,21,25,29)]

We are now ready to compare the assessed stations with respect to the state reference envelope using function compareToStateEnvelope():

#Comparing assessed stations with respect to the state reference envelope
grab_assess <- compareToStateEnvelope(grab_bc, grab_env, m=1.5, 
                                      distances_to_envelope = TRUE)

We can again check that the average of squared distances in the envelope is equal to the envelope variability:

mean(grab_assess$SquaredDist[grab_assess$Envelope])
stateEnvelopeVariability(grab_bc, grab_env)

As before we derive a qualitative assessment from Q values:

grab_assess$Status<-c(ifelse(grab_assess$Q>=0.5,"Inside", "Outside"))

grab_assess$Status[c(1,9:20,21,25,29)]<-"Reference"

grab_assess$Status<- factor(grab_assess$Status, levels = c("Reference","Inside","Outside"))

head(grab_assess)

We then represent the results of the state-based EQA on in the two first dimensions of a Principal Coordinates Analysis:

#Get coordinates of stations in the two first dimensions of the PCoA
pcoa<-ape::pcoa(grab_bc)
PCOA_DIM1_2<-pcoa[["vectors"]][,1:2]
grab_assess$Dim1<-PCOA_DIM1_2[,1]
grab_assess$Dim2<-PCOA_DIM1_2[,2]

#Information about Treatment and Surveys 
grab_assess$Treatment<-Treatment
grab_assess$Surveys<-Surveys


p<-ggplot(grab_assess, aes(x=Dim1, y=Dim2, group=Status)) +
  geom_point(aes(shape=Status, color=Status, size=SquaredDist))+
  scale_shape_manual(values=c(1, 16, 17))+
  scale_color_manual(values=c("#00BFC4", "#00BFC4", "#F8766D"))+
  geom_path(aes(x=Dim1,y=Dim2,group=Treatment),color="grey", arrow = arrow(length = unit(0.03, "cm")))+
  geom_text(grab_assess,mapping=aes(x=Dim1,y=Dim2,label=ID),hjust=-0.3, vjust=-0.3,size=2.5, color="Black")+
  xlab(expression("PCoA1 (13.3%)"))+
  ylab(expression("PCoA2 (9.9%)"))

p<-p+ theme_minimal()+
  theme(
    # Hide some graphical elements
    panel.border = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.line = element_line(colour = "black"))

p

Since we are dealing with temporal data, we can create a bar plot showing the evolution of the squared distance to the centroid of the state reference envelope during the study period for each treatment:

grab_assess$Treatment<-factor(grab_assess$Treatment,c("CTRL1","CTRL2","CTRL3",
                                                      "CD_10","CD_30", 
                                                      "QSD_10","QSD_30",
                                                      "KSD_10"))
grab_assess$Status<-c(ifelse(grab_assess$Q>=0.5,"Inside", "Outside"))
grab_assess$Status<- factor(grab_assess$Status, levels = c("Inside","Outside"))

ggplot(grab_assess, aes(x=Surveys, y=SquaredDist, color=Status, fill=Status)) + 
  geom_bar(stat = "identity",width=0.5)+
  scale_colour_manual(values=c("#00BFC4", "#F8766D"))+
  scale_fill_manual(values=c("#00BFC4", "#F8766D"))+
  scale_y_continuous("Squared distance to the reference envelope")+
  theme_minimal()+
  facet_wrap(~Treatment)+
  theme(axis.text.x = element_text(angle=90))

In addition, we can create the heatmap showing the evolution of the distance between control and dredged stations during the experimentation:

grabbc_sel <- melt(as.matrix(grab_bc),rnames = c("row", "col"))

#selection of the observation of interest (i.e. control stations for heat map row, and dredged stations for heat map column)
grabbc_sel<-grabbc_sel[!grabbc_sel$Var2 %in% c("CTRL1-S1","CTRL2-S1","CTRL3-S1","CTRL4-S1",
                                                "CTRL1-S2","CTRL2-S2","CTRL3-S2","CTRL4-S2",
                                                "CTRL1-S3","CTRL2-S3","CTRL3-S3","CTRL4-S3",
                                                "CTRL1-S4","CTRL2-S4","CTRL3-S4","CTRL4-S4"),]   

grabbc_sel<-grabbc_sel[grabbc_sel$Var1 %in% c("CTRL1-S1","CTRL2-S1","CTRL3-S1","CTRL4-S1",
                                               "CTRL1-S2","CTRL2-S2","CTRL3-S2","CTRL4-S2",
                                               "CTRL1-S3","CTRL2-S3","CTRL3-S3","CTRL4-S3",
                                               "CTRL1-S4","CTRL2-S4","CTRL3-S4","CTRL4-S4"),]

grabbc_sel$Var1<-factor(grabbc_sel$Var1,
                        c("CTRL1-S1","CTRL2-S1","CTRL3-S1","CTRL4-S1",
                          "CTRL1-S2","CTRL2-S2","CTRL3-S2","CTRL4-S2",
                          "CTRL1-S3","CTRL2-S3","CTRL3-S3","CTRL4-S3",
                          "CTRL1-S4","CTRL2-S4","CTRL3-S4","CTRL4-S4"))

grabbc_sel$Var2<-factor(grabbc_sel$Var2,
                        c("KSD_10-S1","CD_10-S1","CD_30-S1","QSD_10-S1","QSD_30-S1",
                          "KSD_10-S2","CD_10-S2","CD_30-S2","QSD_10-S2","QSD_30-S2",
                          "KSD_10-S3","CD_10-S3","CD_30-S3","QSD_10-S3","QSD_30-S3",
                          "KSD_10-S4","CD_10-S4","CD_30-S4","QSD_10-S4","QSD_30-S4"))

ggplot(grabbc_sel, aes(x=Var2, y=Var1, fill=value)) + 
  geom_raster() + 
  scale_fill_viridis_c()+
  theme(axis.text.x=element_text(angle = 90, hjust = 0))

4.5 Trajectory-based EQA

One of the nice features of EQA is the possibility to conduct the quality assessment taking into account ecosystem dynamics. In this case, the conservation target is defined by a set of trajectories conforming a trajectory reference envelope, which in our case are the three control stations:

#Definition of the trajectory reference envelope 
grab_dynenv<-c("CTRL1","CTRL2","CTRL3")

We are now ready to compare the assessed trajectories (those subject to treatment) with respect to the trajectory reference envelope using function compareToTrajectoryEnvelope():

##Comparing assessed trajectories with respect to the trajectory reference envelope
grab_dynassess <- compareToTrajectoryEnvelope(grab_bc,Treatment, grab_dynenv, 
                                              m=1.5, distances_to_envelope = TRUE)

We can also check that the average of squared distances in the envelope is equal to the envelope variability:

mean(grab_dynassess$SquaredDist[grab_dynassess$Envelope])
trajectoryEnvelopeVariability(grab_bc, sites = Treatment, envelope = grab_dynenv)

As before, we can create a bar plot showing the squared distance to the centroid of the trajectory reference envelope during the study period for each treatment:

grab_dynassess$Treatment<-factor(grab_dynassess$Site,c("CTRL1","CTRL2","CTRL3",
                                                      "CD_10","CD_30", 
                                                      "QSD_10","QSD_30",
                                                      "KSD_10"))
grab_dynassess$Status<-c(ifelse(grab_dynassess$Q>=0.5,"Inside", "Outside"))
grab_dynassess$Status<- factor(grab_dynassess$Status, levels = c("Inside","Outside"))

ggplot(grab_dynassess, aes(x=reorder(Treatment, SquaredDist), 
                           y=SquaredDist, color=Status, fill=Status)) + 
  geom_bar(stat = "identity",width=0.5)+
  scale_colour_manual(values=c("#00BFC4", "#F8766D"))+
  scale_fill_manual(values=c("#00BFC4", "#F8766D"))+
  scale_y_continuous("Squared distance to the reference envelope")+
  theme_minimal()+
  theme(axis.text.x = element_text(angle=90))

5. References

• Sturbois, A., De Cáceres, M., Bifolchi, A., Bioret, F., Boyé, A., Gauthier, O., Grall, J., Grémare, A., Labrune, C., Robert, A., Schaal, G., Desroy, N. (2023). Ecological Quality Assessment: a general multivariate framework to report the quality of ecosystems and their dynamics with respect to reference conditions. Ecosphere 14, e4726. https://doi.org/10.1002/ecs2.4726

• Tauran, A., Dubreuil, J., Guyonnet, B., Grall, J., 2020. Impact of fishing gears and fishing intensities on maerl beds: An experimental approach. Journal of Experimental Marine Biology and Ecology 533, 151472. https://doi.org/10.1016/j.jembe.2020.151472

emf-creaf/ecotraj documentation built on April 17, 2025, 5:42 a.m.