Using neonPlantEcology

library(dplyr)
#> Warning: package 'dplyr' was built under R version 4.3.2
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(tidyr)
#> Warning: package 'tidyr' was built under R version 4.3.2
library(stringr)
#> Warning: package 'stringr' was built under R version 4.3.2
library(ggplot2)
#> Warning: package 'ggplot2' was built under R version 4.3.3
library(neonPlantEcology)
library(vegan)
#> Loading required package: permute
#> Loading required package: lattice
#> Warning: package 'lattice' was built under R version 4.3.3
#> This is vegan 2.6-4
library(sf)
#> Warning: package 'sf' was built under R version 4.3.3
#> Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
library(ggpubr)

Doing a community ecology analysis using npe_community_matrix and {vegan}

Setup

Load the tidyverse library along with neonPlantEcology. There is a pre-installed dataset which contains the neonUtilities-generated list object for the Jornada Experimental Range (“JORN”) and the Santa Rita Experimental Range (“SRER”) from domain 14. This can be downloaded directly by uncommenting the following code, or just call data("D14").

# D14 <- npe_site_ids(domain = "D14") |>
#   npe_download_plant_div(sites = .)

data("D14")

Wrangle the data

Use npe_community_matrix to get the data into the needed format. By default, it aggregates annually to the (20m x 20m) plot scale. Use npe_cm_metadata to get all of the information from the rownames into a more interpretable format. Plot centroids, along with additional plot-level metadata, can be loaded with data("plot_centroids").

comm <- npe_community_matrix(D14)
#> Loading required namespace: lubridate
comm[1:4, 1:4]
#>                       2PLANT      ALLIO   AMAC       AMFI
#> JORN_001_plot_2017 1.2500000 0.00000000 0.0000 0.00000000
#> JORN_001_plot_2021 0.5833333 0.00000000 0.0000 2.00000000
#> JORN_001_plot_2022 0.0000000 0.08333333 0.0000 0.08333333
#> JORN_001_plot_2016 0.0000000 0.00000000 2.3125 0.00000000

data("plot_centroids")
plot_centroids <- sf::st_set_geometry(plot_centroids, NULL)
metadata <- npe_cm_metadata(comm) |>
  dplyr::left_join(plot_centroids)
#> Joining with `by = join_by(plotID)`
#> Warning in dplyr::left_join(npe_cm_metadata(comm), plot_centroids): Detected an unexpected many-to-many relationship between `x` and `y`.
#> ℹ Row 1 of `x` matches multiple rows in `y`.
#> ℹ Row 1001 of `y` matches multiple rows in `x`.
#> ℹ If a many-to-many relationship is expected, set `relationship =
#>   "many-to-many"` to silence this warning.

vegan analysis

first we’ll do a species accumulation curve for each site.


# checking the metadata file to see where JORN stops and SRER begins
which(metadata$site == "JORN") |> range()
#> [1]   1 354
which(metadata$site == "SRER") |> range()
#> [1] 355 729

# performing and plotting the species accumulation curve analysis
sp_jorn <- vegan::specaccum(comm[1:152,])
sp_srer <- vegan::specaccum(comm[152:347,])
plot(sp_srer, col = "blue")
plot(sp_jorn, col = "gold", add=T)

NMDS

The following code reproduces figure 3 in Mahood et al 2024 (reference ).

nmds <- metaMDS(comm,trace = F)

nmds_sites <- nmds$points |>
  as_tibble(rownames = "rowname") |>
  left_join(metadata)
#> Joining with `by = join_by(rowname)`

ggplot(nmds_sites, aes(x=MDS1, y=MDS2, color = eventID, size = elevation, shape = site)) +
  geom_point() +
  theme_classic() +
  scale_color_brewer(palette = "Set1") +
  theme(panel.background = element_rect(fill=NA, color= "black"))

getting summary info by site using npe_summary

Species richness by biogeographical origin

di <- npe_summary(D14, scale = "site", timescale = "all")

dj<- di |>
  dplyr::select(site, eventID, starts_with("nspp")) |>
  tidyr::pivot_longer(cols = starts_with("nspp")) |>
  dplyr::filter(!name %in% c("nspp_notexotic", "nspp_total")) |>
  dplyr::transmute(origin = str_remove_all(name, "nspp_"),
            nspp = value,
            site=site)

p1<-ggplot(dj, aes(x=nspp, y=site, fill = origin)) +
  geom_bar(stat = "identity", color = "black", lwd = .2) +
  xlab("Species Richness") +
  ylab("Site")

Relative cover by biogeographical origin

 dk<- di |>
  dplyr::select(site, eventID, starts_with("rel_cover")) |>
  pivot_longer(cols = starts_with("rel_cover")) |>
  filter(!name %in% c("rel_cover_notexotic", "rel_cover_total")) |>
  transmute(origin = str_remove_all(name, "rel_cover_"),
            relative_cover = value,
            site=site)

p2<- ggplot(dk, aes(x=relative_cover, y=site, fill = origin)) +
  geom_bar(stat = "identity", color = "black", lwd = .2) +
  xlab("Relative Cover")

Alpha biogeographical diversity

dl <- di|>
  dplyr::select(site, eventID, starts_with("shannon")) |>
  pivot_longer(cols = starts_with("shannon")) |>
  filter(!name %in% c("shannon_notexotic", "shannon_total", "shannon_family")) |>
  transmute(origin = str_remove_all(name, "shannon_"),
            shannon = value,
            site=site)

p3<- ggplot(dl, aes(x = shannon, y = site, fill = origin)) +
  geom_bar(stat = 'identity', color = "black", lwd = .2) +
  xlab("Shannon-Weiner Diveristy")

Make a multipanel figure

ggpubr::ggarrange(p1, 
                  p2 + theme(axis.title.y = element_blank(), 
                                 axis.text.y = element_blank(), 
                                 axis.ticks.y = element_blank()), 
                  p3 + theme(axis.title.y = element_blank(), 
                                 axis.text.y = element_blank(), 
                                 axis.ticks.y = element_blank()), 
                  nrow=1, ncol=3, 
                  common.legend = TRUE, 
                  widths = c(1.35,1,1))

Family cover through time using npe_longform

data("D14")
lf <- npe_longform(D14, scale = "plot", timescale = "annual")

lf_f <- lf |>
  mutate(family = ifelse(!family %in% c("Poaceae", "Fabaceae"), "Other", family)) |>
  group_by(site, plotID, eventID, family) |>
  summarise(cover = sum(cover, na.rm=T)) |>
  ungroup()
#> `summarise()` has grouped output by 'site', 'plotID', 'eventID'. You can
#> override using the `.groups` argument.
  
ggplot(lf_f, aes(x=eventID, y=cover, fill = family)) +
  geom_boxplot(position="dodge") +
  facet_wrap(~site, scales = "free_x") +
  theme_bw() +
  scale_fill_brewer(palette = "Set1", name = "Family") +
  ylab("Cover") +
  xlab("eventID (Annual Timestep)")