Mixotrophic nanoflagellates from metabarcoding
Suzana G Leles
2026-03-23
Overall goal: interrogate metabarcoding data from metaPR2 in combination with the mixoplankton database to evaluate the relative abundance of mixotrophic nanoflagellates for Lana’s modeling work with mixoCOBALT.
Objective 1: Evaluate the abundance of mixotrophic nanoflagellates relative to diatoms.
Initial set up
Import and merge datasets
Select dataset
# Filtering by gene region and sampling strategy
d2 <- d2 %>% filter(d2$gene_region == "V4", fraction_name == "total")
# Exporting ASV entries and their respective species and genus
unique_combinations_asv_species_genus <- d2 %>%
distinct(asv_code, species, genus)
# Search for Ochromonas
data_filtered <- unique_combinations_asv_species_genus %>%
filter(str_detect(species, "Ochromo"))
# Working with a subset of the variables
deco <- data.frame(asv = d2$asv_code,
species = d2$species,
genus = d2$genus,
nreads = d2$n_reads,
relab = d2$n_reads_pct_corrected,
station = d2$file_code,
type = d2$mixoplankton_functional_type,
size = d2$size.class,
type_size = d2$type_size,
lat = d2$latitude,
long = d2$longitude,
season = d2$season,
date = d2$date,
depth = d2$depth_level,
sizefraction = d2$fraction_name,
temp = d2$temperature,
sal = d2$salinity,
nitrate = d2$NO3)
# Number of samples
nlevels(factor(deco$station))## [1] 1376# Number of unique ASVs
nlevels(factor(deco$asv))## [1] 20595# Quick look at the dataset
head(deco)## asv species genus nreads
## 1 5b23c637cf Dino-Group-I-Clade-1_X_sp. Dino-Group-I-Clade-1_X 74
## 2 49ea4db573 Choreotrichida_XX_sp. Choreotrichida_XX 78
## 3 72de8a9444 Prorocentrum_sp. Prorocentrum 148
## 4 916dee5b0d Picozoa_XXXXX_sp. Picozoa_XXXXX 112
## 5 aaf18167c8 Pleurostomatida_X_sp. Pleurostomatida_X 60
## 6 8ce3d7c533 Hydrurales_XX_sp. Hydrurales_XX 7
## relab station type size type_size lat
## 1 1.043665e-03 ERR2172161 parasite <NA> parasite_NA 81.9255
## 2 1.100079e-03 ERR2172161 protozooplankton <NA> protozooplankton_NA 81.9255
## 3 2.087329e-03 ERR2172161 not assessed <NA> not assessed_NA 81.9255
## 4 1.579601e-03 ERR2172161 protozooplankton <NA> protozooplankton_NA 81.9255
## 5 8.462146e-04 ERR2172161 protozooplankton <NA> protozooplankton_NA 81.9255
## 6 9.872504e-05 ERR2172161 not assessed <NA> not assessed_NA 81.9255
## long season date depth sizefraction temp sal nitrate
## 1 131.129 autumn 2012-09-04 surface total -1.5 30.6 NA
## 2 131.129 autumn 2012-09-04 surface total -1.5 30.6 NA
## 3 131.129 autumn 2012-09-04 surface total -1.5 30.6 NA
## 4 131.129 autumn 2012-09-04 surface total -1.5 30.6 NA
## 5 131.129 autumn 2012-09-04 surface total -1.5 30.6 NA
## 6 131.129 autumn 2012-09-04 surface total -1.5 30.6 NA# Functional types
levels(factor(deco$type))## [1] "CM" "eSNCM" "GNCM" "not assessed"
## [5] "parasite" "phytoplankton" "protozooplankton" "pSNCM"# Change this since phytoplankton are actually all diatoms
deco <- deco %>% mutate(type = fct_recode(type, "diatoms" = "phytoplankton"))
# Check size classes
levels(factor(deco$size))## [1] "meso" "meso-macro" "micro" "micro-meso" "nano"
## [6] "nano-micro" "not recorded" "pico"# Separate CM into nanoCM and otherCM
data_eco <- deco %>% mutate(typecobalt = case_when(type == "CM" & size == "nano" ~ "nanoCM", type == "CM" ~ "otherCM", TRUE ~ type))
data_eco %>% count(type, size, typecobalt)## type size typecobalt n
## 1 CM micro otherCM 4351
## 2 CM nano nanoCM 9281
## 3 CM nano-micro otherCM 66
## 4 CM pico otherCM 265
## 5 eSNCM meso eSNCM 607
## 6 eSNCM meso-macro eSNCM 39
## 7 eSNCM micro eSNCM 647
## 8 eSNCM micro-meso eSNCM 846
## 9 eSNCM nano eSNCM 52
## 10 eSNCM not recorded eSNCM 15
## 11 GNCM micro GNCM 903
## 12 not assessed <NA> not assessed 121450
## 13 parasite <NA> parasite 144083
## 14 diatoms <NA> diatoms 42487
## 15 protozooplankton <NA> protozooplankton 164035
## 16 pSNCM micro pSNCM 218
## 17 pSNCM nano pSNCM 25Interrogating only mixotrophic nanoflagellates and diatoms
# Compute the total number of reads for nanoCM and diatoms by station
data_nanoCM_diatoms <- data_eco %>% filter(typecobalt %in% c("diatoms", "nanoCM"))
levels(factor(data_nanoCM_diatoms$typecobalt))## [1] "diatoms" "nanoCM"head(data_nanoCM_diatoms)## asv species genus nreads relab
## 1 0b1827b046 Dolichomastix_tenuilepis Dolichomastix 93 0.0013116326
## 2 d94a0e224a Karlodinium_veneficum Karlodinium 75 0.0010577683
## 3 d246da9ef2 Fragilariopsis_cylindrus Fragilariopsis 212 0.0029899583
## 4 dc4597b06a Pterosperma_cristatum Pterosperma 492 0.0069389597
## 5 844d4850e1 Chaetoceros_neogracilis Chaetoceros 91 0.0012834255
## 6 0d307e6ac2 Teleaulax_amphioxeia Teleaulax 55 0.0007756967
## station type size type_size lat long season date
## 1 ERR2172161 CM nano CM_nano 81.9255 131.129 autumn 2012-09-04
## 2 ERR2172161 CM nano CM_nano 81.9255 131.129 autumn 2012-09-04
## 3 ERR2172161 diatoms <NA> phytoplankton_NA 81.9255 131.129 autumn 2012-09-04
## 4 ERR2172161 CM nano CM_nano 81.9255 131.129 autumn 2012-09-04
## 5 ERR2172161 diatoms <NA> phytoplankton_NA 81.9255 131.129 autumn 2012-09-04
## 6 ERR2172161 CM nano CM_nano 81.9255 131.129 autumn 2012-09-04
## depth sizefraction temp sal nitrate typecobalt
## 1 surface total -1.5 30.6 NA nanoCM
## 2 surface total -1.5 30.6 NA nanoCM
## 3 surface total -1.5 30.6 NA diatoms
## 4 surface total -1.5 30.6 NA nanoCM
## 5 surface total -1.5 30.6 NA diatoms
## 6 surface total -1.5 30.6 NA nanoCM# Calculate the relative number of reads for each species by station
data_nanoCM_diatoms <- data_nanoCM_diatoms %>% group_by(station) %>%
mutate(total_reads_group = sum(nreads, na.rm = TRUE),
rel_reads = nreads/total_reads_group) %>%
ungroup()
# Checking the relative abundance sums up to 1 - it does!
test <- data_nanoCM_diatoms %>%
group_by(station) %>%
summarize(sum_rel_reads = sum(rel_reads, na.rm = TRUE))
max(test$sum_rel_reads)## [1] 1min(test$sum_rel_reads)## [1] 1# Calculate the relative number of reads for each type by station
dsum <- data_nanoCM_diatoms %>%
group_by(station, typecobalt) %>%
summarise(rel_reads_type = sum(rel_reads, na.rm = TRUE),
temp = first(temp), sal = first(sal), nitrate = first(nitrate), lat = first(lat), long = first(long),
.groups = "drop")
# Checking the relative abundance sums up to 1 - it does!
test2 <- dsum %>%
group_by(station) %>%
summarize(sum_rel_reads = sum(rel_reads_type, na.rm = TRUE))
max(test2$sum_rel_reads)## [1] 1min(test2$sum_rel_reads)## [1] 1Relative abundance of mixotrophic nanoflagellates and diatoms across all samples
# Check how many samples we have
n_distinct(dsum$station)## [1] 1372# Make a barplot showing nanoCM and diatom relative abundance by station
ggplot(dsum, aes(x = station, y = rel_reads_type, fill = typecobalt)) +
geom_col() +
labs(x = "Station", y = "Relative abundance") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
Global map: mixo nanos relative to diatoms
# Get world map data
world <- map_data("world")
# Now plot the relative abundance of mixotrophic nanoflagellates (nanoCM) in the global ocean
p <- ggplot() +
geom_map(data = world, map = world, aes(x = long, y = lat, map_id = region),
colour = "darkgray", fill = "darkgray") +
scale_y_continuous(expand = c(0,0)) +
scale_x_continuous(expand = c(0,0)) +
theme(axis.title=element_blank(),
axis.text=element_blank(),
axis.ticks=element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_rect(fill = "white", size = 0.5, colour = 'darkgray'),
legend.position = "right",
legend.title = element_text(size = 16),
legend.text = element_text(size = 16),
legend.key.width = unit(.5, "cm"),
legend.key.height = unit(.95, "cm"),
strip.background = element_rect(colour="darkgray", fill="white"),
strip.text = element_text(size = 19),
plot.title = element_text(hjust = 0, face = "bold", size = 22)
)
col_com <- c("#762A83", "#C2A5CF", "#FFEE99", "#ACD39E", "#1B7837")
# Plot the ratio between mixotrophic nanoflagellates and diatoms abundance in the global ocean
p + geom_point(data = dsum %>% filter(typecobalt == "nanoCM"), aes(x = long, y = lat, fill = rel_reads_type), size = 4, shape = 21) +
scale_fill_viridis_c(option = "viridis", guide = guide_colorbar(barwidth = 15, barheight = 1), name = "Relative abundance of mixo nano") +
theme(legend.key = element_blank()) +
theme(legend.title = element_text(size = 15, colour = "black"),
legend.position = "top")
# Now in a 2D plot temperature vs latitude
ggplot(dsum %>% filter(typecobalt == "nanoCM"), aes(x = temp, y = lat, fill = rel_reads_type)) +
geom_point(shape = 21, size = 4) +
scale_fill_viridis_c(option = "viridis", guide = guide_colorbar(barwidth = 15, barheight = 1), name = "Relative abundance of mixo nano") +
ylab("Latitude") +
xlab(expression("Temperature ( " * degree * "C)")) +
mytheme +
theme(legend.key = element_blank()) +
theme(legend.title = element_text(size = 15, colour = "black"),
legend.position = "top")
# Now in a 2D plot temperature versus nitrate
ggplot(dsum %>% filter(typecobalt == "nanoCM"), aes(x = temp, y = nitrate, fill = rel_reads_type)) +
geom_point(shape = 21, size = 4) +
scale_fill_viridis_c(option = "viridis", guide = guide_colorbar(barwidth = 15, barheight = 1), name = "Relative abundance of mixo nano") +
ylab(expression("Nitrate (" * mu * "M)")) +
xlab(expression("Temperature ( " * degree * "C)")) +
mytheme +
theme(legend.key = element_blank()) +
theme(legend.title = element_text(size = 15, colour = "black"),
legend.position = "top")