EMODnetMapCleanStyle <- theme_bw() +
  theme(axis.text = element_blank(),
        axis.title = element_blank())

EMODnetMapNormalStyle <- 
  theme(axis.text = element_text(size = 12, color = "blue"),
        axis.title = element_text(size = 14, color = "blue"),
        legend.text = element_text(size = 10),
        legend.position = "right",
        legend.key.width = unit(10,"mm"),
        legend.key.height = unit(10, "mm"),
        plot.background = element_rect(fill = "grey90", colour = "white", size = 0),
        panel.border = element_blank()
)

Introduction

Documentation

Github project: https://github.com/wstolte/mixoplanktonFractions

Data extraction

Temporal scale

Data from 1995 until now are considered for the current product.

Geographical scale

The regions that were selected were assembled from the intersection of the IHO regions and the EEZ from the different countries. These subregions have ID’s that can be used in the WFS query to the EMODnet Biology database.

regions <- sf::st_read(quiet = T, dsn = "../data/derived_data/greater_north_sea-selection_from_eez-iho_v4.geojson")
regions %>% ggplot() +
  geom_sf(fill = "blue", color = "white") +
  geom_sf_text(aes(label = mrgid), size = 2, color = "white") +
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank())

Strategy

For each of these subregions, data extraction was done in a two-step approach.

  1. Extract phytoplankton data by selecting only observations linkted to the trait “phytoplankton”.
  2. The unique datasets obtained in the first steps were inspected, and suitable datasets, which were expected to contain phytoplankton, were selected. The suitable datasets were then downloaded completely (without selection for traits)

Extraction of complete relevant datasets

Data were extracted from EMODnet Biology via WFS, using the requestData.R script. These data are stored locally and taken into this script here.

# all2Data <- read_delim(file.path(dataDir, "all2Data.csv"), delim = ";")
load(file.path(dataDir, "all2Data.Rdata"))

In total, this resulted in 2024331 observations

Per subregion, the number of observations/km2 looks like this:

regionN <- all2Data %>%
  group_by(mrgid) %>% summarize(n = n()) %>% ungroup() %>%
  mutate(mrgid = as.numeric(mrgid))

regions %>% right_join(regionN, by = c(mrgid = "mrgid")) %>%
  mutate(n_km2 = n/area_km2) %>%
  ggplot() +
  geom_sf(aes(fill = n_km2)) +
  scale_fill_viridis_c()

Join with trophic mode table

Schneider et al. (submitted to Journal of Biological Data, data submitted to WoRMS) has defined trophic mode for about 1500 species of plankton. A first subdivision is in

  • phytoplankton - phototrophic, no phagotrophy known
  • mixoplankton - capable of phototrophy and phagotrophy
  • zooplankton - phagotrophic, no phototrophy known.

Mixoplankton is subdivided in different classes, but at this moment, we will not use this.

# read information available in WoRMS traits. 
# This is a local copy

trophy <- read_delim("../data/raw_data/traits.csv", delim = ";", guess_max = 10000)
FALSE 
FALSE -- Column specification --------------------------------------------------------
FALSE cols(
FALSE   ScientificName = col_character(),
FALSE   AphiaID = col_character(),
FALSE   Trophy = col_character(),
FALSE   typeMX = col_character(),
FALSE   source = col_character(),
FALSE   sourceType = col_character(),
FALSE   contributor = col_character()
FALSE )
trophyData <- all2Data %>% 
  mutate(AphiaID = sub("http://marinespecies.org/aphia.php[?]p=taxdetails&id=", "", aphiaid, )) %>% 
  full_join(trophy, by = c(AphiaID = "AphiaID")) %>% 
  select(
    mrgid, datasetID, datecollected, decimallongitude, decimallatitude,
    scientificname, AphiaID, Trophy, typeMX
  ) %>%
  drop_na(Trophy) %>%
  mutate(
    season = case_when(
      month(datecollected) %in% c(6,7,8) ~ "summer", 
      month(datecollected) %in% c(3,4,5) ~ "spring", 
      month(datecollected) %in% c(9,10,11) ~ "autumn", 
      month(datecollected) %in% c(12, 1, 2) ~ "winter"
    )
  ) %>% drop_na(season) %>%
  mutate(season = factor(season, levels = c("winter", "spring", "summer", "autumn")))

Fraction of mixoplankton (first results)

As a first approach, the fraction of mixoplankton/(phytoplankto + mixoplankton) is plotted on a map:

plottrophy <- trophyData %>%
  group_by(decimallongitude, decimallatitude, Trophy, season) %>%
  summarize(n = n()) %>%
  filter(Trophy %in% c("mixotroph", "phototroph")) %>%
  spread(Trophy, n) %>%
  mutate(fMix = mixotroph/(mixotroph + phototroph)) %>%
  drop_na(fMix)
FALSE `summarise()` has grouped output by 'decimallongitude', 'decimallatitude', 'Trophy'. You can override using the `.groups` argument.
write_delim(plottrophy, file.path(csvDir, "trophydata.csv"), delim = ";")

require(rworldxtra)
FALSE Loading required package: rworldxtra
data(countriesHigh)
world <- st_as_sf(countriesHigh) %>% 
  filter(REGION == "Europe") %>%
  filter(ADMIN %in% c("Denmark", "Germany", "Netherlands", "Belgium", "United Kingdom", "France", "Norway", "Sweden", "Ireland")) %>%
  select(ADMIN) %>%
  st_crop(xmin = -15, xmax = 17, ymin = 45, ymax = 65)
FALSE although coordinates are longitude/latitude, st_intersection assumes that they are planar
ggplot() +
  stat_summary_2d(data = plottrophy, aes(x = decimallongitude, y = decimallatitude, z=fMix), fun=mean, binwidth = c(1,0.5)) +
  geom_sf(data = world, fill = "transparent", color = "grey", alpha = 0) +
  scale_fill_viridis_c(option = "viridis") +
  coord_sf(xlim = c(-15, 17), ylim = c(45, 65), expand = T, clip = "on") +
  facet_wrap(~ season) +
  EMODnetMapCleanStyle

This first analysis shows highest fraction of mixoplankton in summer. Even in summer, the fraction of mixoplankton is still low along the Dutch, Belgian and French coast. This is in line with Schneider et al. (2020a). It is also in line with a general idea that mixoplankton, at least the constitutive mixoplankton, appear in post-bloom clear water conditions, which in the North Sea are associated with offshore, possibly stratified waters.

Limitations of this analysis - discussion

Incomplete classification

At the moment, only 1500 species could be automatically classified as either phytoplankton or mixoplankton. This means that we are missing a number of species in this analysis, because the data contain more than 5000 species names. Part of those may be heterotrophic, which does not influence the maps above. If we assume that the fraction of mixoplankton is not very different in these unknown species as compared to the known species, the bias is not very strong. In the current analysis, the number of observations where species could not be classified was lower than the number of phytoplankton observations (figure ]@ref(fig:.

all2Data %>% 
  mutate(AphiaID = sub("http://marinespecies.org/aphia.php[?]p=taxdetails&id=", "", aphiaid, )) %>% 
  full_join(trophy, by = c(AphiaID = "AphiaID")) %>% 
  select(
    mrgid, datasetID, datecollected, decimallongitude, decimallatitude,
    scientificname, AphiaID, Trophy, typeMX
  ) %>%
  group_by(Trophy) %>%
  summarize(`nr of observations` = n()) %>%
  ggplot(aes(Trophy, `nr of observations`)) + geom_col(aes(fill = Trophy))
Distribution of observations over known mixotrophs, phagotrophs, phototrophs and unknown trophy according to Schneider et al., (2020b).

Distribution of observations over known mixotrophs, phagotrophs, phototrophs and unknown trophy according to Schneider et al., (2020b).

Bias in species per dataset

Not all datasets address the same species. For example, the Plankton Recorder data contain species that are retained on a 60 µm filter, which excludes a number of small species. A possible way to get a better grip of the effect of this bias is to only regards species that occur in all (or most) datasets. This has not been undertaken yet.

References