Light Penertration Proxies

library(ggOceanMaps)
library(ggOceanMapsData)
library(ncdf4)
library(raster)
library(RColorBrewer)
library(dplyr)
library(plyr)
library(patchwork)
library(knitr)
library(kableExtra)
library(scales)
library(utils)

large <- theme(legend.title = element_text(size=15),
        legend.text = element_text(size=15),
        axis.title = element_text(size=15),
        axis.text = element_text(size=15),
        strip.text = element_text(size=15))

no_legend <- theme(legend.position = "none")
vert_x <- theme(axis.text.x = element_text(angle = 60))
letter_strip <- theme(strip.text = element_text(hjust = 0, vjust = 1.5, margin = margin(0.2,0,0.2,0, "cm"), size=15), strip.background = element_blank())

Spectral.colors <- colorRampPalette(rev(brewer.pal(11, "Spectral")))
ex <- extent(c(-180, 180, 50, 90)) # Eye Index data extant

# Function to convert scientific notation to exponential notation
scientific_10 <- function(x) { parse(text=gsub("e", "%*%10^", scales::label_scientific()(x)))}

#  Stereographic map axis title
hx <- cbind(transform_coord(lon = seq(-135, 180, by = 45), lat = rep(48, 8), proj.out = 3995), label = c("135°W", "90°W", "45°W", "0°", "45°E", "90°E", "135°E", "180°"), hjust = c(rep(0.8, 3), 0, rep(0.2, 3), 0.5))
vx <-  cbind(transform_coord(lon = rep(180, 4), lat = c(60, 70, 80, 90), proj.out = 3995), label = c("60°N", "70°N", "80°N", "90°N"))

# Eye Index data
er <- read.csv("../csv/fulldataset-20240227.csv")
er <- cbind(transform_coord(lon=er$Longitude, lat=er$Latitude, proj.out = 3995), er)
# Show species eye ratio
er$ER <- cut(er$Species, breaks=c(0, 0.0001, 0.25, 0.5, 0.75, 1.1), right=FALSE, labels=c("0", "0-0.25", "0.25-0.5", "0.5-0.75", "0.75-1"))
# Only keep data with > 100 individuals
er <- subset(er, Total>=100)

# Download Arctic land shapfile
url <- "https://github.com/MikkoVihtakari/ggOceanMapsLargeData/raw/master/data/arctic_land.rda"
destfile <- paste(tempdir(), "arctic_land.rda", sep="//")
download.file(url, destfile)

Lee et al. (2005) developed a model for vertical transmittance of solar radiation in the upper ocean. The model uses total absorption, particulate backscatter, water depth, and solar zenith angle as inputs. In this test, I used the total absorption and backscatter at 488 nm from the entire mission composite climatology to calculate the vertical transmittance (Tvis) for the visible spectrum at surface, 25 m, euphotic depth, and seafloor. Tvis is the ratio between downwelling solar radiation at the surface, Evis(0), and at depth, Evis(z). The total absorption, backscatter, and euphotic depth were downloaded from OceanColor Web L3 browser. The seafloor bathymetry were downloaded from ETOPO Global Relief Model. The complete references for Lee’s algorithms can be found here.

Data from OceanColor Web

IOP Total Absorption at 488 nm (MODIS_Aqua Entire Mission Composite climatology)

nc <- nc_open("../OceanColorData/MODIS_Aqua/IOP_a_488/AQUA_MODIS.20020704_20220930.L3m.CU.IOP.a_488.4km.nc")
ng <- ncvar_get(nc, varid = "a_488")
a488 <- raster(t(ng), xmn=-180, xmx=180, ymn=-90, ymx=90, crs = "+proj=longlat +datum=WGS84") %>% crop(ex)
nc_close(nc)

r <- projectRaster(a488, crs = projection(ggOceanMapsData::arctic_bathy))
df <- rasterToPoints(r) %>% as.data.frame

p <- basemap(limits=50, shapefiles = "Arctic")+ 
  geom_raster(data=df, aes(x=x, y=y, fill=layer))+
  scale_fill_gradientn(colours=Spectral.colors(100), trans="log10")

p_a488 <- reorder_layers(p)+
  geom_text(data=hx, aes(x=lon, y=lat, label=label, hjust=hjust), colour="gray30")+
  geom_text(data=vx, aes(x=lon, y=lat, label=label), colour="gray30")+
  geom_point(data=er, aes(x=lon, y=lat, colour=ER))+
  scale_colour_viridis_d()+
  labs(fill=expression("a488"~(m^-1)), colour="Eye Index", title=expression("a"~~~"Total Absorption at 488 nm"))+
  guides(fill = guide_colourbar(order=1, barheight=10), colour = guide_legend(order=2))+
  coord_sf(clip = 'off')+
  theme(legend.key = element_rect(fill = "gray80", colour = NA))

IOP Total Backscatter at 488 nm (MODIS_Aqua Entire Mission Composite climatology)

nc <- nc_open("../OceanColorData/MODIS_Aqua/IOP_bb_488/AQUA_MODIS.20020704_20220930.L3m.CU.IOP.bb_488.4km.nc")
ng <- ncvar_get(nc, varid = "bb_488")
bb488 <- raster(t(ng), xmn=-180, xmx=180, ymn=-90, ymx=90, crs = "+proj=longlat +datum=WGS84") %>% crop(ex)
nc_close(nc)

r <- projectRaster(bb488, crs = projection(ggOceanMapsData::arctic_bathy))
df <- rasterToPoints(r) %>% as.data.frame

p <- basemap(limits=50, shapefiles = "Arctic")+ 
  geom_raster(data=df, aes(x=x, y=y, fill=layer))+
  scale_fill_gradientn(colours=Spectral.colors(100), trans="log10")

p_bb488 <- reorder_layers(p)+
  geom_text(data=hx, aes(x=lon, y=lat, label=label, hjust=hjust), colour="gray30")+
  geom_text(data=vx, aes(x=lon, y=lat, label=label), colour="gray30")+
  geom_point(data=er, aes(x=lon, y=lat, colour=ER))+
  scale_color_viridis_d()+
  labs(fill=expression("bb488"~(m^-1)), colour="Eye Index", title=expression("b"~~~"Total Backscatter at 488 nm"))+
  guides(fill = guide_colourbar(order=1, barheight=10), colour = guide_legend(order=2))+
  coord_sf(clip = 'off')+
  theme(legend.key = element_rect(fill = "gray80", colour = NA))

p_a488+p_bb488

Fig.S1. IOP total absorption at 488 nm and total backscatter at 488 nm based on MODIS_Aqua entire mission composite climatology.

Euphotic depth (MODIS_Aqua Entire Mission Composite climatology) from Lee Algorithm (ZLEE)

nc <- nc_open("../OceanColorData/MODIS_Aqua/Zeu_lee/A20021852022059.L3m_CU_ZLEE_Zeu_lee_4km.nc")
ng <- ncvar_get(nc, varid = "Zeu_lee")
Zeu <- raster(t(ng), xmn=-180, xmx=180, ymn=-90, ymx=90, crs = "+proj=longlat +datum=WGS84") %>% crop(ex)
nc_close(nc)

r <- projectRaster(Zeu, crs = projection(ggOceanMapsData::arctic_bathy))
df <- rasterToPoints(r) %>% as.data.frame

p <- basemap(limits=50, shapefiles = "Arctic")+ 
  geom_raster(data=df, aes(x=x, y=y, fill=layer))+
  scale_fill_gradientn(colours=Spectral.colors(100) %>% rev)

p_zeu <- reorder_layers(p)+
  geom_text(data=hx, aes(x=lon, y=lat, label=label, hjust=hjust), colour="gray30")+
  geom_text(data=vx, aes(x=lon, y=lat, label=label), colour="gray30")+
  geom_point(data=er, aes(x=lon, y=lat, colour=ER))+
  scale_color_viridis_d()+
  labs(fill="Zeu (m)", colour="Eye Index", title=expression("a"~~~"Euphotic depth"))+
  guides(fill = guide_colourbar(order=1, barheight=10), colour = guide_legend(order=2))+
  coord_sf(clip = 'off')+
  theme(legend.key = element_rect(fill = "gray80", colour = NA))

Photosynthetically available radiation (MODIS_Aqua Entire Mission Composite climatology)

nc <- nc_open("../OceanColorData/MODIS_Aqua/PAR/AQUA_MODIS.20020704_20220930.L3m.CU.PAR.par.4km.nc")
ng <- ncvar_get(nc, varid = "par")
par <- raster(t(ng), xmn=-180, xmx=180, ymn=-90, ymx=90, crs = "+proj=longlat +datum=WGS84")  %>% crop(ex)
nc_close(nc)

r <- projectRaster(par, crs = projection(ggOceanMapsData::arctic_bathy))
df <- rasterToPoints(r) %>% as.data.frame
df$layer[df$layer>quantile(df$layer, .99)] <- quantile(df$layer, .99) # make 99 percentile as maximum value

p <- basemap(limits=50, shapefiles = "Arctic")+ 
  geom_raster(data=df, aes(x=x, y=y, fill=layer))+
  scale_fill_gradientn(colours=Spectral.colors(100))

p_par <- reorder_layers(p)+
  geom_text(data=hx, aes(x=lon, y=lat, label=label, hjust=hjust), colour="gray30")+
  geom_text(data=vx, aes(x=lon, y=lat, label=label), colour="gray30")+
  geom_point(data=er, aes(x=lon, y=lat, colour=ER))+
  scale_color_viridis_d()+
  labs(fill=expression("PAR"~(E~m^-2*d^-1)), colour="Eye Index", title=expression("b"~~~"Photosynthetically available radiation"))+
  guides(fill = guide_colourbar(order=1, barheight=10), colour = guide_legend(order=2))+
  coord_sf(clip = 'off')+
  theme(legend.key = element_rect(fill = "gray80", colour = NA))

p_zeu+p_par

ETOPO 2022 Global Relief Model

nc <- nc_open("../ETOPO2022/ETOPO_2022_v1_60s_N90W180_bed.nc")
ng <- ncvar_get(nc, varid = "z")
bathy <- raster(t(ng[,10800:1]), xmn=-180, xmx=180, ymn=-90, ymx=90, crs = "+proj=longlat +datum=WGS84")  %>% crop(ex)
nc_close(nc)

bathy <- calc(bathy, fun=function(x){x[x>0]<-NA; return(-x)})

bathy <- projectRaster(bathy, a488)

p <- basemap(limits=50, shapefiles = "Arctic", bathymetry=TRUE)

p_bathy <- reorder_layers(p)+
  geom_text(data=hx, aes(x=lon, y=lat, label=label, hjust=hjust), colour="gray30")+
  geom_text(data=vx, aes(x=lon, y=lat, label=label), colour="gray30")+
  geom_point(data=er, aes(x=lon, y=lat, colour=ER))+
  scale_color_viridis_d()+
  labs(fill="Depth (m)", colour="Eye Index", title=expression("a"))+
  coord_sf(clip = 'off')

p <- basemap(limits=50, shapefiles = "Arctic", bathymetry=FALSE)

p_region <- reorder_layers(p)+
  geom_text(data=hx, aes(x=lon, y=lat, label=label, hjust=hjust), colour="gray30")+
  geom_text(data=vx, aes(x=lon, y=lat, label=label), colour="gray30")+
  geom_point(data=er, aes(x=lon, y=lat, colour=Region))+
  scale_color_viridis_d()+
  labs(fill="Depth (m)", colour="Region", title=expression("b"))+
  coord_sf(clip = 'off')

p_bathy+p_region

Figure 1. Distribution map showing the location and eye index (EIS) and region of each sample from the Arctic Ostracods Database in the Arctic Ocean.

Lee et al. (2005) Journal of Geophysical Research, 110(C2)

Table 2. Values of the Model Parameters for Kvis(z)

X0 = -0.057
X1 = 0.482
X2 = 4.221
sigma0 = 0.183
sigma1 = 0.702
sigma2 = -2.567
alpha0 = 0.090
alpha1 = 1.465
alpha2 = -0.667

Vertical transmittance (Tvis) at surface, 25-m depth, and euphotic depth at zero solar zenith angle

# Function for two-parameter model to estimate vertical transmittance (Tvis)
Tvis_fun <- function(a488, bb488, z, theta){
  # Equation 9a
  K1 = (X0 + X1 * a488^0.5 + X2 * bb488) * (1 + alpha0 * sin(theta * pi/180))
  # Equation 9b
  K2 = (sigma0 + sigma1 * a488 + sigma2 * bb488) * (alpha1 +alpha2 * cos(theta * pi/180))
  # Equation 7
  # Attenuation coefficient
  Kvis = K1 + K2/(1 + z)^0.5
  # Equation 10
  # Vertical transmittance
  Tvis <- exp(-Kvis * z)
  return(Tvis)
}

Tvis at surface

Tvis <- Tvis_fun(a488, bb488, 0, 0)

r <- projectRaster(Tvis, crs = projection(ggOceanMapsData::arctic_bathy))
df <- rasterToPoints(r) %>% as.data.frame

p <- basemap(limits=50, shapefiles = "Arctic")+ 
  geom_raster(data=df, aes(x=x, y=y, fill=layer))+
  scale_fill_gradientn(colours=Spectral.colors(100))

p_tvis0 <- reorder_layers(p)+
  geom_text(data=hx, aes(x=lon, y=lat, label=label, hjust=hjust), colour="gray30")+
  geom_text(data=vx, aes(x=lon, y=lat, label=label), colour="gray30")+
  geom_point(data=er, aes(x=lon, y=lat, colour=ER))+
  scale_color_viridis_d()+
  labs(fill=expression(italic(T[vis](0))), colour="Eye Index", title=expression("a"~~~"Vertical transmittance at surface"))+
  guides(fill = guide_colourbar(order=1, barheight=10), colour = guide_legend(order=2))+
  coord_sf(clip = 'off')+
  theme(legend.key = element_rect(fill = "gray80", colour = NA))

Tvis at 25 m

Tvis <- Tvis_fun(a488, bb488, 25, 0)

r <- projectRaster(Tvis, crs = projection(ggOceanMapsData::arctic_bathy))
df <- rasterToPoints(r) %>% as.data.frame

p <- basemap(limits=50, shapefiles = "Arctic")+ 
  geom_raster(data=df, aes(x=x, y=y, fill=layer))+
  scale_fill_gradientn(colours=Spectral.colors(100))

p_tvis25 <- reorder_layers(p)+
  geom_text(data=hx, aes(x=lon, y=lat, label=label, hjust=hjust), colour="gray30")+
  geom_text(data=vx, aes(x=lon, y=lat, label=label), colour="gray30")+
  geom_point(data=er, aes(x=lon, y=lat, colour=ER))+
  scale_color_viridis_d()+
  labs(fill=expression(italic(T[vis](25))), colour="Eye Index", title=expression("b"~~~"Vertical transmittance at 25-m"))+
  guides(fill = guide_colourbar(order=1, barheight=10), colour = guide_legend(order=2))+
  coord_sf(clip = 'off')+
  theme(legend.key = element_rect(fill = "gray80", colour = NA))

Tvis at 25 m * PAR

r <- projectRaster(Tvis*par, crs = projection(ggOceanMapsData::arctic_bathy))
df <- rasterToPoints(r) %>% as.data.frame

p <- basemap(limits=50, shapefiles = "Arctic")+ 
  geom_raster(data=df, aes(x=x, y=y, fill=layer))+
  scale_fill_gradientn(colours=Spectral.colors(100))

p_par25 <- reorder_layers(p)+
  geom_text(data=hx, aes(x=lon, y=lat, label=label, hjust=hjust), colour="gray30")+
  geom_text(data=vx, aes(x=lon, y=lat, label=label), colour="gray30")+
  geom_point(data=er, aes(x=lon, y=lat, colour=ER))+
  scale_color_viridis_d()+
  labs(fill=expression("PAR"~(25)~("E"~m^-2*d^-1)), colour="Eye Index", title=expression("c"~~~"Photosynthetically available radiation at 25-m"))+
  guides(fill = guide_colourbar(order=1, barheight=10), colour = guide_legend(order=2))+
  coord_sf(clip = 'off')+
  theme(legend.key = element_rect(fill = "gray80", colour = NA))

Tvis at euphotic depth

Tvis <- Tvis_fun(a488, bb488, Zeu, 0)
# remove the top 1% of data
Tvis[Tvis>quantile(Tvis, 0.99)] <- NA

r <- projectRaster(Tvis, crs = projection(ggOceanMapsData::arctic_bathy))
df <- rasterToPoints(r) %>% as.data.frame

p <- basemap(limits=50, shapefiles = "Arctic")+ 
  geom_raster(data=df, aes(x=x, y=y, fill=layer))+
  scale_fill_gradientn(colours=Spectral.colors(100))

p_tvis_zeu <- reorder_layers(p)+
  geom_text(data=hx, aes(x=lon, y=lat, label=label, hjust=hjust), colour="gray30")+
  geom_text(data=vx, aes(x=lon, y=lat, label=label), colour="gray30")+
  geom_point(data=er, aes(x=lon, y=lat, colour=ER))+
  scale_color_viridis_d()+
  labs(fill=expression(italic(T[vis](Z[eu]))), colour="Eye Index", title=expression("d"~~~"Vertical transmittance at Zeu"))+
  guides(fill = guide_colourbar(order=1, barheight=10), colour = guide_legend(order=2))+
  coord_sf(clip = 'off')+
  theme(legend.key = element_rect(fill = "gray80", colour = NA))

(p_tvis0+p_tvis25)/(p_par25+p_tvis_zeu)

p <- p_tvis25+labs(title=expression("a"~~~"Vertical transmittance at 25-m"))
p+p_par

Fig. 5 Geographical distribution of light availability proxies. Heat maps indicate (A) Tvis (Vertical transmittance) at 25-m water depth and (B) PAR (Photosynthetically available radiation). Note that sites in the high Arctic with white background indicates permanently ice-covered sites and thus no light availability data is available for those sites. Eye index (EIS) shown on the maps.

Tvis at bottom

Tvis <- Tvis_fun(a488, bb488, bathy, 0)

r <- projectRaster(Tvis, crs = projection(ggOceanMapsData::arctic_bathy))
df <- rasterToPoints(r) %>% as.data.frame

p <- basemap(limits=50, shapefiles = "Arctic")+ 
  geom_raster(data=df, aes(x=x, y=y, fill=layer))+
  scale_fill_gradientn(colours=Spectral.colors(100), trans= trans_new("5rt", function(x)x^(0.2), function(x)x^5))

p_tvis_bottom <- reorder_layers(p)+
  geom_text(data=hx, aes(x=lon, y=lat, label=label, hjust=hjust), colour="gray30")+
  geom_text(data=vx, aes(x=lon, y=lat, label=label), colour="gray30")+
  geom_point(data=er, aes(x=lon, y=lat, colour=ER))+
  scale_color_brewer(type="seq", palette="YlOrBr", direction=-1)+
  labs(fill=expression(italic(T[vis](z))), colour="Eye Index", title=expression("a"~~~"Vertical transmittance at bottom"))+
  guides(fill = guide_colourbar(order=1, barheight=10), colour = guide_legend(order=2))+
  coord_sf(clip = 'off')+
  theme(legend.key = element_rect(fill = "gray80", colour = NA))

Tvis at bottom * PAR

r <- projectRaster(Tvis*par, crs = projection(ggOceanMapsData::arctic_bathy))
df <- rasterToPoints(r) %>% as.data.frame

p <- basemap(limits=50, shapefiles = "Arctic")+ 
  geom_raster(data=df, aes(x=x, y=y, fill=layer))+
  scale_fill_gradientn(colours=Spectral.colors(100), trans= trans_new("5rt", function(x)x^(0.2), function(x)x^5))

p_par_bottom <- reorder_layers(p)+
  geom_text(data=hx, aes(x=lon, y=lat, label=label, hjust=hjust), colour="gray30")+
  geom_text(data=vx, aes(x=lon, y=lat, label=label), colour="gray30")+
  geom_point(data=er, aes(x=lon, y=lat, colour=ER))+
  scale_color_brewer(type="seq", palette="YlOrBr", direction=-1)+
  labs(fill=expression("PAR"~(italic(z))~("E"~m^-2*d^-1)), colour="Eye Index", title=expression("b"~~~"Photosynthetically available radiation at bottom"))+
  guides(fill = guide_colourbar(order=1, barheight=10), colour = guide_legend(order=2))+
  coord_sf(clip = 'off')+
  theme(legend.key = element_rect(fill = "gray80", colour = NA))

p_tvis_bottom+p_par_bottom

Figure 5 Geographical distribution of light availability proxies. Heat maps indicate (A) Vertical transmittance and (B) Photosynthetically available radiation at sampling depth [Tvis (z) and PAR (z), respectively]. Note that for sites in the high Arctic the white background indicates permanently ice-covered sites and thus no light availability data is available for those sites. Eye index shown on the maps.

Extract a488, bb488, Zeu and calculate TVis

We extracted the gridded a488 and bb488 by sampling coordinates. For the site outside of the grid coverage, we extract the mean grid values within a 100-km radius buffer. The remaining sites not matched by spatial grids or buffers are under the permanent sea ice cover (between 84.3 and 89.98 latitudes); therefore, their Tvis are set to zero.

# Data distribution
loc <- er[, c("Longitude", "Latitude")]
coordinates(loc) <- c("Longitude", "Latitude")
projection(loc) <- "+proj=longlat +datum=WGS84"

# IOP Total Absorption at 488 nm (m^-1)
ta <- extract(a488, loc)
na <- is.na(ta)
mi <- extract(a488, loc[na,], buffer=100000, fun=mean, na.rm=TRUE)
ta0 <- ta
ta0[na] <- mi

# IOP Total Backscatter at 488 nm (m^-1)
tb <- extract(bb488, loc)
na <- is.na(tb)
mi <- extract(bb488, loc[na,], buffer=100000, fun=mean, na.rm=TRUE)
tb0 <- tb
tb0[na] <- mi

# Euphotic depth (m)
ed <- extract(Zeu, loc)
na <- is.na(ed)
mi <- extract(Zeu, loc[na,], buffer=100000, fun=mean, na.rm=TRUE)
ed0 <- ed
ed0[na] <- mi

# Photosynthetically available radiation (einstein m^-2 day^-1)
pa <- extract(par, loc)
na <- is.na(pa)
mi <- extract(par, loc[na,], buffer=100000, fun=mean, na.rm=TRUE)
pa0 <- pa
pa0[na] <- mi

#TVIS at sampling depth
tv <- Tvis_fun(ta, tb, er$Water.Depth, 0)
tv0 <- Tvis_fun(ta0, tb0, er$Water.Depth, 0)
# All the missing values are under permanent ice cover
#coordinates(loc[is.na(tv0),])[,2] %>% as.vector %>% range
tv0[is.na(tv0)] <- 0

out <- cbind(er, "a488"= ta, "bb488" = tb, "Zeu" = ed, "PAR" = pa, "Tvis" = tv, "a488_buffer"= ta0, "bb488_buffer" = tb0, "Zeu_buffer" = ed0, "PAR_buffer" = pa0, "Tvis_buffer" = tv0)
out$Tvis_PAR <- with(out, Tvis_buffer*PAR_buffer)

with(out, cbind(Water.Depth, Latitude))%>%cor(use="pairwise.complete.obs")

##             Water.Depth  Latitude
## Water.Depth   1.0000000 0.6740313
## Latitude      0.6740313 1.0000000

with(subset(out, Tvis_PAR>0), cbind(Water.Depth, Latitude, log_Tvis_PAR=log(Tvis_PAR), log_Tvis=log(Tvis_buffer)))%>%cor(use="pairwise.complete.obs")

##              Water.Depth   Latitude log_Tvis_PAR   log_Tvis
## Water.Depth    1.0000000  0.5776202   -0.9680071 -0.9673197
## Latitude       0.5776202  1.0000000   -0.5781727 -0.5755277
## log_Tvis_PAR  -0.9680071 -0.5781727    1.0000000  0.9999782
## log_Tvis      -0.9673197 -0.5755277    0.9999782  1.0000000

with(subset(out, Tvis_PAR>0), cbind(Tvis_buffer, Tvis_PAR))%>%cor(use="pairwise.complete.obs")

##             Tvis_buffer  Tvis_PAR
## Tvis_buffer   1.0000000 0.9887171
## Tvis_PAR      0.9887171 1.0000000

write.csv(out, file= "../csv/fulldataset_Tvis.csv", row.names=FALSE)

library(ggplot2)
library(reshape2)
library(car)
library(mgcv)
library(doBy)
library(lme4)

ml <- melt(out, id.vars=c("Number", "Water.Depth", "Latitude", "Longitude", "Region", "Tvis", "Zeu", "PAR", "Total"), measure.vars = c("Abundance", "Species", "Genus"), value.name = "Eye_Ratio")
ml$Tvis_PAR <- with(ml, Tvis*PAR)
ml$variable <- factor(ml$variable, labels=c("Individual", "Species", "Genus"))

ml0 <- melt(out, id.vars=c("Number", "Water.Depth", "Latitude", "Longitude", "Region", "Tvis_buffer", "Zeu_buffer", "PAR_buffer", "Total"), measure.vars = c("Abundance", "Species", "Genus"), value.name = "Eye_Ratio")
ml0$Tvis_PAR <- with(ml0, Tvis_buffer*PAR_buffer)
ml0$variable <- factor(ml0$variable, labels=c("Individual", "Species", "Genus"))
ml0$withEye <- with(ml0, round(Total*Eye_Ratio, 0))
ml0$noEye <- with(ml0, Total-withEye)

Fitting water depth to eye indices

Generalized Linear Mixed-Effects Models (GLMM)

# Eye indices as response variables
# Water depth as fixed factor (scaled due to large depth differences)
# Region as random factor

s1 <- lapply(splitBy(~variable, ml0)[2], FUN=function(x)glmer(Eye_Ratio~Water.Depth + (1 | Region), family = binomial, data=x) %>% summary)

s2 <- lapply(splitBy(~variable, ml0)[2], FUN=function(x)glmer(Eye_Ratio~Tvis_buffer+(1|Region), family = binomial, data=x) %>% summary)

s3 <- lapply(splitBy(~variable, ml)[2], FUN=function(x)glmer(Eye_Ratio~Tvis_PAR+(1|Region), family = binomial, data=x) %>% summary)

s4 <- lapply(splitBy(~variable, ml0)[2], FUN=function(x)glmer(Eye_Ratio~Latitude+(1|Region), family = binomial, data=x) %>% summary)

rbind(s1$Species$coefficients, NA, s2$Species$coefficients, NA,s3$Species$coefficients, NA,s4$Species$coefficients)%>%kable(caption="Table x. GLMM with binomial distribution on eye indices with water depth, Tvis(z), PAR(z), or latitude as a fixed factor and region as a random factor") %>% kable_classic(full_width = F, html_font = "Cambria", font_size=16)

Table x. GLMM with binomial distribution on eye indices with water depth, Tvis(z), PAR(z), or latitude as a fixed factor and region as a random factor

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	-1.5936893	0.6462953	-2.465884	0.0136676
Water.Depth	-0.1568411	0.0668435	-2.346391	0.0189562
	NA	NA	NA	NA
(Intercept)	-5.6985606	0.7230835	-7.880916	0.0000000
Tvis_buffer	7.8985666	1.6047239	4.922072	0.0000009
	NA	NA	NA	NA
(Intercept)	-5.7732118	0.8049629	-7.172023	0.0000000
Tvis_PAR	0.3447811	0.0757513	4.551486	0.0000053
	NA	NA	NA	NA
(Intercept)	4.9242499	4.7519433	1.036260	0.3000807
Latitude	-0.1316256	0.0693939	-1.896790	0.0578557

# According to Crawley (2012) The R Book Chapter 16: Propositional data. https://onlinelibrary.wiley.com/doi/10.1002/9781118448908.ch16
# Combine numbers of specimens with eyes and without eye as two-vector response variable.
# Water depth as fixed factor (scaled due to large depth differences)
# Region as random factor

s1 <- lapply(splitBy(~variable, ml0)[2], FUN=function(x){
  y <- cbind(x$withEye, x$noEye)
  glmer(y~Water.Depth + (1 | Region), family = binomial, data=x) %>% summary
  })

s2 <- lapply(splitBy(~variable, ml0)[2], FUN=function(x){
  y <- cbind(x$withEye, x$noEye)
  glmer(y~Tvis_buffer + (1 | Region), family = binomial, data=x) %>% summary
  })

s3 <- lapply(splitBy(~variable, ml0)[2], FUN=function(x){
  y <- cbind(x$withEye, x$noEye)
  glmer(y~Tvis_PAR + (1 | Region), family = binomial, data=x) %>% summary
  })

s4 <- lapply(splitBy(~variable, ml0)[2], FUN=function(x){
  y <- cbind(x$withEye, x$noEye)
  glmer(y~Latitude + (1 | Region), family = binomial, data=x) %>% summary
  })

rbind(s1$Species$coefficients, NA, s2$Species$coefficients, NA,s3$Species$coefficients, NA,s4$Species$coefficients)%>%kable(caption="Table x. GLMM with binomial distribution on the numbers of specimen with eyes and without eye using water depth, Tvis(z), PAR(z), or latitude as a fixed factor and region as a random factor") %>% kable_classic(full_width = F, html_font = "Cambria", font_size=16)

Table x. GLMM with binomial distribution on the numbers of specimen with eyes and without eye using water depth, Tvis(z), PAR(z), or latitude as a fixed factor and region as a random factor

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	-1.0159913	0.1220916	-8.321547	0
Water.Depth	-0.0037792	0.0000750	-50.376018	0
	NA	NA	NA	NA
(Intercept)	-1.6816979	0.1532962	-10.970250	0
Tvis_buffer	2.0325379	0.0450720	45.095357	0
	NA	NA	NA	NA
(Intercept)	-1.5920158	0.1290458	-12.336831	0
Tvis_PAR	0.0679283	0.0016582	40.964064	0
	NA	NA	NA	NA
(Intercept)	7.9757850	0.1794886	44.436160	0
Latitude	-0.1323085	0.0014590	-90.686313	0

ml0$Lab <- factor(ml0$variable, labels = c(NA, "a", NA))

# Menzies, George & Rowe (1968) Nature 217:93-95
# Isopod eye ratio @ 65 degree N
Water.Depth = c(2, 20, 200, 1000, 2000, 4000, 6000)
Eye_Ratio <- c(1, .65, .46, .17, .08, .02, 0)
iso <- data.frame(Water.Depth, Eye_Ratio, variable = "Species", Lab = "a")

p_depth_eye1 <- ggplot(data=subset(ml0, variable=="Species"), aes(x=Water.Depth, y=Eye_Ratio))+
  geom_point(aes(colour=Region))+
  stat_smooth(data=subset(ml0, variable=="Species"), method="glm", formula = y~x, method.args = list(family = "binomial"))+
  geom_point(data=iso, aes(x=Water.Depth, y=Eye_Ratio), colour = "red")+
  geom_path(data=iso, aes(x=Water.Depth, y=Eye_Ratio), colour = "red", linetype=2)+
  coord_flip()+
  scale_x_reverse()+
  scale_colour_viridis_d()+
  labs(x="Water depth (m)", y="Eye Index")+
  facet_wrap(~Lab, scales="free", nrow=2)+
  theme_bw() %+replace% large %+replace% vert_x %+replace% letter_strip 

ml0$Lab <- factor(ml0$variable, labels = c(NA, "b", NA))
iso$Lab <- "b"

p_depth_eye2 <- ggplot(data=subset(ml0, variable=="Species"), aes(x=Water.Depth, y=Eye_Ratio))+
  geom_point(aes(colour=Region))+
  stat_smooth(method="glm", formula = y~x, method.args = list(family = "binomial"))+
  geom_point(data=iso, aes(x=Water.Depth, y=Eye_Ratio), colour = "red")+
  geom_path(data=iso, aes(x=Water.Depth, y=Eye_Ratio), colour = "red", linetype=2)+
  coord_flip(xlim = c(500, 0))+
  scale_x_reverse()+
  scale_colour_viridis_d()+
  labs(x="Water depth (m)", y="Eye Index")+
  facet_wrap(~Lab, scales="free", nrow=2)+
  theme_bw() %+replace% large %+replace% vert_x %+replace% letter_strip 

p1 <- p_depth_eye1+labs(y="Eye Index") + no_legend
p2 <- p_depth_eye2+theme(axis.title.y = element_blank()) 
p1+p2

Figure 3. Bathymetric distribution of eye index values in the Arctic Ocean. Panel a shows the entire depth range, and panel b only shows the samples from the shallow region (0-500 m). The blue line shows significant relationship with 95% confidence interval (shaded area) based on the Generalized Linear Mixed-Effect Model (GLMM) with binomial distribution. Red symbols and dashed lines show isopod eye index from Menzies et al. (1968).

eii_r <- subset(ml0, variable =="Species")
eii_r$Lab <- eii_r$Region
eii_r$Lab <- factor(eii_r$Lab, levels = c("Chukchi", "Beaufort", "Baffin", "Greenland", "Barents", "Kara", "Laptev", "Siberian"), labels = c("a   Chukchi", "b   Beaufort", "c   Baffin", "d   Greenland", "e   Barents", "f   Kara", "g   Laptev", "h   Siberian"))
ggplot(data=eii_r, aes(x=Water.Depth, y=Eye_Ratio))+
  geom_point()+
  coord_flip(xlim = c(500, 0))+
  scale_x_reverse()+
  scale_colour_viridis_d()+
  labs(x="Water depth (m)", y="Eye Index")+
  facet_wrap(~Lab, nrow=2)+
  theme_bw() %+replace% large %+replace% vert_x %+replace% letter_strip 

Figure 4. Bathymetric distribution of eye index values computed by species method (EIS) in different regions (A–F) of the Arctic Ocean. Only samples from relatively shallow regions (0-500 m) are presented to show the regional difference of eye index in shallow water among the various Arctic Ocean coastal seas.

Scatter plots with Tvis at Sampling Depth

s <- lapply(splitBy(~variable, ml)[2], FUN=function(x)glmer(Eye_Ratio~Tvis+(1|Region), family = binomial, data=x) %>% summary)

lapply(s, FUN=function(x)cbind(Variable=rownames(x$coef), data.frame(x$coef)))%>%ldply%>%kable(caption="Table x. GLMM with binomial distribution on eye indices with water depth as fixed factor and regions as random factor") %>% kable_classic(full_width = F, html_font = "Cambria", font_size=16)

ml$Lab <- factor(ml$variable, labels = c("A   Individual", "B   Species", NA))

ggplot(data=subset(ml, variable=="Species"), aes(x=Tvis, y=Eye_Ratio))+
  geom_point(aes(colour=Region))+
  stat_smooth(data=subset(ml, variable=="Species"), method="glm", formula = y~x, method.args = list(family = "binomial"))+
  scale_colour_viridis_d()+
  labs(x=expression("Vertical transmittance"~italic(T[vis](z))~"at sampling depth"), y="Eye Index")+
  #facet_wrap(~Lab, scales="free")+
  theme_bw() %+replace% large %+replace% vert_x %+replace% letter_strip

ml0$Lab <- factor(ml0$variable, labels = c(NA, "a", NA))

p_Tvis_eye <- ggplot(data=subset(ml0, variable=="Species"), aes(x=Tvis_buffer, y=Eye_Ratio))+
  geom_point(aes(colour=Region))+
  stat_smooth(data=subset(ml0, variable=="Species"), method="glm", formula = y~x, method.args = list(family = "binomial"))+
  scale_colour_viridis_d()+
  coord_flip()+
  scale_x_log10(labels = scientific_10)+
  labs(x=expression("Vertical transmittance"~italic(T[vis](z))~"at sampling depth with 100-km Buffer"), y="Eye Index")+
  facet_wrap(~Lab, scales="free")+
  theme_bw() %+replace% large %+replace% vert_x %+replace% letter_strip

Scatter plots with PAR at Sampling Depth

ml$Lab <- factor(ml$variable, labels = c("A   Individual", "B   Species", NA))

# E m-2 d-1 = 6.97 x 10^12 μmol quanta m-2 s-1
ggplot(data=subset(ml, variable=="Species"), aes(x=Tvis*PAR*6.97E12, y=Eye_Ratio))+
  geom_point(aes(colour=Region))+
  stat_smooth(data=subset(ml, variable!="Genus"), method="glm", formula = y~x, family = binomial, method.args = list(family = "binomial"))+
  coord_flip()+
  scale_x_log10(labels = scientific_10)+
  scale_colour_viridis_d()+
  labs(x=expression("PAR"~(italic(z))~(mu*mol~"quanta"~m^-2*s^-1)), y="Eye Index")+
  #facet_wrap(~Lab, scales="free")+
  theme_bw() %+replace% large %+replace% letter_strip

ml0$Lab <- factor(ml0$variable, labels = c(NA, "b", NA))

# E m-2 d-1 = 6.97 x 10^12 μmol quanta m-2 s-1
p_PAR_eye <- ggplot(data=subset(ml0, variable=="Species"), aes(x=Tvis_PAR*6.97E12, y=Eye_Ratio))+
  geom_point(aes(colour=Region))+
  stat_smooth(data=subset(ml0, variable=="Species"), method="glm", formula = y~x, family = binomial, method.args = list(family = "binomial"))+
  coord_flip()+
  scale_x_log10(labels = scientific_10)+
  scale_colour_viridis_d()+
  labs(x="Photosynthetically available radiation (μmol quanta m^-2 s^-1) at sampling depth with 100-km Buffer", y="Eye Index")+
  facet_wrap(~Lab, scales="free")+
  theme_bw() %+replace% large %+replace% letter_strip

p1 <- p_Tvis_eye+labs(x=expression("Vertical transmittance"~italic(T[vis](z))), y="Eye Index")
p2 <- p_PAR_eye+labs(x=expression("PAR"~(italic(z))~(mu*mol~"quanta"~m^-2*s^-1)), y="Eye Index") + no_legend
p1/p2

Figure 6 Eye index and light availability proxies of the vertical transmittance and photosynthetically available radiation at sampling depth [Tvis (z) and PAR (z), respectively]. The blue line shows significant relationship with 95% confidence interval (shaded area) based on Generalized Linear Mixed-Effect Model (GLMM) with binomial distribution.

ml0$Lab <- factor(ml0$variable, labels = c("A   Individual", "B   Species", NA))

ggplot(data=subset(ml0, variable=="Species"), aes(x=Latitude, y=Eye_Ratio))+
  geom_point(aes(colour=Region))+
  stat_smooth(method="glm", formula = y~x, method.args = list(family = "binomial"))+
  scale_colour_viridis_d()+
  labs(x=expression("Latitude"~(degree*N)), y="Eye Index")+
  #facet_wrap(~Lab, scales="free")+
  theme_bw() %+replace% large %+replace% letter_strip

Figure 7 Latitudinal gradient of eye index. The blue line shows significant relationship with 95% confidence interval (shaded area) based on Generalized Linear Mixed-Effect Model (GLMM) with binomial distribution.

# Eye indices as response variables
# Water depth and latitude as fixed factors 
# Region as random factor

s <- lapply(splitBy(~variable, ml0)[2], FUN=function(x)glmer(Eye_Ratio~Water.Depth + Latitude + (1 | Region), family = binomial, data=x) %>% summary)

lapply(s, FUN=function(x)cbind(Variable=rownames(x$coef), data.frame(x$coef)))%>%ldply%>%kable(caption="Table x. GLMM with binomial distribution on eye indices with water depth and latitude as fixed factors and regions as random factor") %>% kable_classic(full_width = F, html_font = "Cambria", font_size=16)

Table x. GLMM with binomial distribution on eye indices with water depth and latitude as fixed factors and regions as random factor

.id	Variable	Estimate	Std..Error	z.value	Pr…z..
Species	(Intercept)	13.3146169	6.6440766	2.003983	0.0450719
Species	Water.Depth	-0.1293521	0.0608561	-2.125541	0.0335415
Species	Latitude	-0.2226279	0.0995347	-2.236685	0.0253069

ml0$Lab <- factor(ml0$variable, labels = c(NA, "a", NA))

p1 <- ggplot(data=subset(ml0, variable=="Species"), aes(x=Water.Depth, y=Eye_Ratio))+
  geom_point(aes(colour=Region))+
  stat_smooth(data=subset(ml0, variable=="Species"), method="glm", formula = y~x, method.args = list(family = "binomial"))+
  coord_flip()+
  scale_x_reverse()+
  scale_colour_viridis_d()+
  labs(x="Water depth (m)", y="Eye Index")+
  facet_wrap(~Lab, scales="free")+
  theme_bw() %+replace% large %+replace% letter_strip

ml0$Lab <- factor(ml0$variable, labels = c(NA, "b", NA))

p2 <- ggplot(data=subset(ml0, variable=="Species"), aes(x=Latitude, y=Eye_Ratio))+
  geom_point(aes(colour=Region))+
  stat_smooth(data=subset(ml0, variable=="Species"), method="glm", formula = y~x, method.args = list(family = "binomial"))+
  coord_flip()+
  scale_colour_viridis_d()+
  labs(x=expression("Latitude"~(degree*N)), y="Eye Index")+
  facet_wrap(~Lab, scales="free")+
  theme_bw() %+replace% large %+replace% letter_strip %+replace% no_legend

p1/p2

Fig. x. Eye index (EIS) as functions of the sampling depth and latitude. EIS is computed by species method.

# Eye indices as response variables
# Latitude (scaled due to large depth differences) and PAR(z) as fixed factors 
# Region as random factor

s <- lapply(splitBy(~variable, subset(ml0, Tvis_PAR>0))[2], FUN=function(x)glmer(Eye_Ratio~Tvis_PAR + Latitude + (1 | Region), family = binomial, data=x) %>% summary)

lapply(s, FUN=function(x)cbind(Variable=rownames(x$coef), data.frame(x$coef)))%>%ldply%>%kable(caption="Table x. GLMM with binomial distribution on eye indices with photosynthetically available radiation (PAR) at sampling depth and latitude as fixed factors and regions as random factor") %>% kable_classic(full_width = F, html_font = "Cambria", font_size=16)

Table x. GLMM with binomial distribution on eye indices with photosynthetically available radiation (PAR) at sampling depth and latitude as fixed factors and regions as random factor

.id	Variable	Estimate	Std..Error	z.value	Pr…z..
Species	(Intercept)	0.1665264	7.2754968	0.0228887	0.9817391
Species	Tvis_PAR	0.2686358	0.0680225	3.9492163	0.0000784
Species	Latitude	-0.0815454	0.1036276	-0.7869087	0.4313353

ml0$Lab <- factor(ml0$variable, labels = c(NA, "a", NA))

# E m-2 d-1 = 6.97 x 10^12 μmol quanta m-2 s-1
p1 <- ggplot(data=subset(ml0, variable=="Species"), aes(x=Tvis_PAR*6.97E12, y=Eye_Ratio))+
  geom_point(aes(colour=Region))+
  stat_smooth(data=subset(ml0, variable=="Species"), method="glm", formula = y~x, method.args = list(family = "binomial"))+
  coord_flip()+
  scale_x_log10(labels = scientific_10)+
  scale_colour_viridis_d()+
  labs(x=expression("PAR"~(italic(z))~(mu*mol~"quanta"~m^-2*s^-1)), y="Eye Index")+
  facet_wrap(~Lab, scales="free")+
  theme_bw() %+replace% large %+replace% letter_strip

ml0$Lab <- factor(ml0$variable, labels = c(NA, "b", NA))

p2 <- ggplot(data=subset(ml0, variable=="Species"), aes(x=Latitude, y=Eye_Ratio))+
  geom_point(aes(colour=Region))+
  stat_smooth(data=subset(ml0, variable=="Species"), method="glm", formula = y~x, method.args = list(family = "binomial"))+
  coord_flip()+
  scale_colour_viridis_d()+
  labs(x=expression("Latitude"~(degree*N)), y="Eye Index")+
  facet_wrap(~Lab, scales="free")+
  theme_bw() %+replace% large %+replace% letter_strip %+replace% no_legend

p1/p2

Fig. x. Eye index (EIS) as functions of the photosynthetically available radiation (PAR) at sampling depth and latitude. Panel (A), (C) are computed by abundance method; Panel (B), (D) are computed by species method.

summary(glmer(Tvis_buffer~Water.Depth+(1|Region), family = gaussian(link = "log"), data=subset(out, Tvis_buffer>0.001)))$coef %>% kable(caption = "Table x. GLMM with gaussian distribution and logarithmic link function on Tvis. Data with Tvis < 0.1% were deleted before analysis.") %>% kable_classic(full_width = F, html_font = "Cambria", font_size=16)

Table x. GLMM with gaussian distribution and logarithmic link function on Tvis. Data with Tvis

	Estimate	Std. Error	t value	Pr(>|z|)
(Intercept)	-0.4142500	0.4480243	-0.9246152	0.3551661
Water.Depth	-0.1828439	0.0066353	-27.5564365	0.0000000

p_Tvis_depth <- ggplot(data=subset(out, Tvis_buffer!=0), aes(x=Water.Depth, y=Tvis_buffer))+
  geom_point(aes(colour=Region))+
  stat_smooth(method="glm", formula = y~x, method.args = list(family = gaussian(link = "log")))+
  coord_flip()+
  scale_x_reverse()+
  scale_y_continuous(labels = scientific_10)+
  scale_colour_viridis_d()+
  labs(x=expression("Water depth"~(m)), y=expression(italic(T[vis](z))~"at sampling depth"), tag = "a")+
  theme_bw() %+replace% large %+replace% vert_x 

summary(glmer(Tvis_PAR~Water.Depth+(1|Region), family = gaussian(link = "log"), data=subset(out, Tvis_buffer>0.001)))$coef%>% kable(caption = "Table x. GLMM with gaussian distribution and logarithmic link function on PAR(z) at sampling depth. Data with Tvis < 0.1% were deleted before analysis.") %>% kable_classic(full_width = F, html_font = "Cambria", font_size=16)

Table x. GLMM with gaussian distribution and logarithmic link function on PAR(z) at sampling depth. Data with Tvis

	Estimate	Std. Error	t value	Pr(>|z|)
(Intercept)	2.6958175	0.2456224	10.97545	0
Water.Depth	-0.1699915	0.0057960	-29.32909	0

# E m-2 d-1 = 6.97 x 10^12 μmol quanta m-2 s-1
p_PAR_depth <- ggplot(data=subset(out, Tvis_buffer!=0), aes(x=Water.Depth, y=Tvis_PAR*6.97E12))+
  geom_point(aes(colour=Region))+
  stat_smooth(method="glm", formula = y~x, method.args = list(family = gaussian(link = "log")))+
  coord_flip()+
  scale_x_reverse()+
  scale_y_continuous(labels = scientific_10)+
  scale_colour_viridis_d()+
  labs(x=expression("Water depth"~(m)), y=expression("PAR"~(italic(z))~(mu*mol~"quanta"~m^-2*s^-1)), tag="b")+
  theme_bw() %+replace% large %+replace% vert_x 

p1 <- p_Tvis_depth + labs(y=expression(italic(T[vis](z)))) + no_legend
p2 <- p_PAR_depth + labs(y=expression("PAR"~(italic(z))~(mu*mol~"quanta"~m^-2*s^-1))) 
p1+p2

Fig. X. Vertical transmittance (Tvis) and photosynthetically available radiation (PAR) at sampling depth as functions of water depth. Data with Tvis = 0% were removed before conducting GLMM with with gaussian distribution and logarithmic link

p1_log <- p1+scale_y_log10(labels = scientific_10) + stat_smooth(method="glm", formula = y~x)
p2_log <- p2+scale_y_log10(labels = scientific_10) + stat_smooth(method="glm", formula = y~x)
p1_log+p2_log

Fig. 8. Vertical transmittance (Tvis) and photosynthetically available radiation (PAR) at sampling depth as functions of water depth. Data with Tvis = 0% were removed before conducting GLMM with with gaussian distribution and logarithmic link

summary(glmer(Tvis_buffer~Water.Depth+(1|Region), family = gaussian(link = "log"), data=subset(out, Tvis_buffer>0.001)))$coef %>% kable(caption = "Table x. GLMM with gaussian distribution and logarithmic link function on Tvis. Data with Tvis < 0.1% were deleted before analysis.") %>% kable_classic(full_width = F, html_font = "Cambria", font_size=16)

Table x. GLMM with gaussian distribution and logarithmic link function on Tvis. Data with Tvis

	Estimate	Std. Error	t value	Pr(>|z|)
(Intercept)	-0.4142500	0.4480243	-0.9246152	0.3551661
Water.Depth	-0.1828439	0.0066353	-27.5564365	0.0000000

p_Tvis_depth <- ggplot(data=subset(out, Tvis_buffer>0.001), aes(x=Water.Depth, y=Tvis_buffer))+
  geom_point(aes(colour=Region))+
  stat_smooth(method="glm", formula = y~x, method.args = list(family = gaussian(link = "log")))+
  coord_flip()+
  scale_x_reverse()+
  scale_colour_viridis_d()+
  labs(x=expression("Water depth"~(m)), y=expression(italic(T[vis](z))~"at sampling depth"), tag = "a")+
  theme_bw() %+replace% large %+replace% vert_x 

summary(glmer(Tvis_PAR~Water.Depth+(1|Region), family = gaussian(link = "log"), data=subset(out, Tvis_buffer>0.001)))$coef%>% kable(caption = "Table x. GLMM with gaussian distribution and logarithmic link function on PAR(z) at sampling depth. Data with Tvis < 0.1% were deleted before analysis.") %>% kable_classic(full_width = F, html_font = "Cambria", font_size=16)

Table x. GLMM with gaussian distribution and logarithmic link function on PAR(z) at sampling depth. Data with Tvis

	Estimate	Std. Error	t value	Pr(>|z|)
(Intercept)	2.6958175	0.2456224	10.97545	0
Water.Depth	-0.1699915	0.0057960	-29.32909	0

# E m-2 d-1 = 6.97 x 10^12 μmol quanta m-2 s-1
p_PAR_depth <- ggplot(data=subset(out, Tvis_buffer>0.001), aes(x=Water.Depth, y=Tvis_PAR*6.97E12))+
  geom_point(aes(colour=Region))+
  stat_smooth(method="glm", formula = y~x, method.args = list(family = gaussian(link = "log")))+
  coord_flip()+
  scale_x_reverse()+
  scale_y_continuous(labels = scientific_10)+
  scale_colour_viridis_d()+
  labs(x=expression("Water depth"~(m)), y=expression("PAR"~(italic(z))~(mu*mol~"quanta"~m^-2*s^-1)), tag="a")+
  theme_bw() %+replace% large %+replace% vert_x 

p2 <- p_PAR_depth + labs(y=expression("PAR"~(italic(z))~(mu*mol~"quanta"~m^-2*s^-1))) 

p4 <- p2_log + labs(tag="b") + no_legend

p2/p4

Fig. 8 Photosynthetically available radiation at sampling depth [PAR (z)] as functions of water depth. (A) for the top 70 m; (B) for full bathymetric range. Data with Tvis (z) less than 0.1% were removed before conducting GLMM with gaussian distribution and logarithmic link function.