What is the role of water table in the ash fertilization growth response?

R code for the comparative analysis of contrasting drained Scots pine peatland forests in northern Finland.

Libraries

library(broom)
library(colorspace)
library(cowplot)
library(dplR)

This is dplR version 1.7.8.
dplR is part of openDendro https://opendendro.org.
New users can visit https://opendendro.github.io/dplR-workshop/ to get started.

library(dplyr)

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(ggplot2)
library(ggpubr)

Attaching package: 'ggpubr'

The following object is masked from 'package:cowplot':

    get_legend

library(lmtest)

Loading required package: zoo

Attaching package: 'zoo'

The following object is masked from 'package:dplR':

    time<-

The following objects are masked from 'package:base':

    as.Date, as.Date.numeric

library(lubridate)

Attaching package: 'lubridate'

The following object is masked from 'package:cowplot':

    stamp

The following objects are masked from 'package:base':

    date, intersect, setdiff, union

library(mgcv)

Loading required package: nlme

Attaching package: 'nlme'

The following object is masked from 'package:dplyr':

    collapse

This is mgcv 1.9-3. For overview type 'help("mgcv-package")'.

library(patchwork)

Attaching package: 'patchwork'

The following object is masked from 'package:cowplot':

    align_plots

library(readxl)
library(rstatix)

Attaching package: 'rstatix'

The following object is masked from 'package:stats':

    filter

library(tibble)
library(tidyr)
library(utils)

Define colours

# Vesi: Muhos site with shallow WT
vesi1 <- "#00B5E2"
vesi2 <- lighten(vesi1, 0.2)
vesi3 <- lighten(vesi1, 0.5)

# Valo: Sievi site with deep WT
valo1 <- "#FF8200"
valo2 <- lighten(valo1, 0.2)
valo3 <- lighten(valo1, 0.5)

# Vadelma: Highlighting
vadelma1 <- "#CC0082"
vadelma2 <- lighten(vadelma1, 0.5)

# Kallio: Lines
kallio1 <- "#54585A"
kallio2 <- lighten(kallio1, 0.3)

Figure 2: Annual July-August WT levels

Read and prepare WT data:

WT_data <- bind_rows(
  read_excel("DATA WT Muhos.xlsx") |>
    mutate(across(`1`:`50`, as.numeric), 
           site = factor(site, levels = "Ansasaari", labels = "Shallow WT"), 
           treatment = factor(treatment, levels = c("Control", "Fertilization"), labels = c("Control", "Treatment")),  
           date = as.Date(paste(day, month, year, sep = "-"), format = "%d-%m-%Y")) |> 
    select(site, treatment, year, month, date, `1`:`50`) |>  
    pivot_longer(cols = `1`:`50`, names_to = "WT_pipe", values_to = "WT_cm") |>   
    mutate(WT_cm = WT_cm * -1,
           WT_pipe = as.numeric(WT_pipe)),
  read_excel("DATA WT Sievi.xlsx") |>
    mutate(across(`1`:`50`, as.numeric), 
           site = factor(site, levels = "Nevajärvi", labels = "Deep WT"), 
           treatment = factor(treatment, levels = c("Control", "Fertilization"), labels = c("Control", "Treatment")), 
           date = as.Date(paste(day, month, year, sep = "-"), format = "%d-%m-%Y")) |> 
    select(site, treatment, year, month, date, `1`:`50`) |>  
    pivot_longer(cols = `1`:`50`, names_to = "WT_pipe", values_to = "WT_cm") |>   
    mutate(WT_cm = WT_cm * -1,
           WT_pipe = as.numeric(WT_pipe))
) |> 
  filter(month %in% c(7, 8))

Warning: There were 22 warnings in `mutate()`.
The first warning was:
ℹ In argument: `across(`1`:`50`, as.numeric)`.
Caused by warning:
! NAs introduced by coercion
ℹ Run `dplyr::last_dplyr_warnings()` to see the 21 remaining warnings.

WT_summary <- WT_data |>
  group_by(site, treatment, year) |> 
  summarize(WT_cm = mean(WT_cm, na.rm = TRUE, .groups = "drop"))

`summarise()` has grouped output by 'site', 'treatment'. You can override using
the `.groups` argument.

wt_thinn_dry <- WT_summary |> filter(site == "Deep WT", treatment == "Treatment", year == "2009") |> pull(WT_cm)
wt_thinn_wet <- WT_summary |> filter(site == "Shallow WT", treatment == "Treatment", year == "2010") |> pull(WT_cm)
wt_af_dry <- WT_summary |> filter(site == "Deep WT", treatment == "Treatment", year == "2014") |> pull(WT_cm)
wt_af_wet <- WT_summary |> filter(site == "Shallow WT", treatment == "Treatment", year == "2014") |> pull(WT_cm)

Plot:

WT_data |> 
  filter(year < 2021 & year > 2006) |>
  mutate(site.treatment = factor(interaction(site, treatment), levels = c( "Shallow WT.Control", "Shallow WT.Treatment", "Deep WT.Control", "Deep WT.Treatment"))) |> 
  ggplot(aes(x = year, y = WT_cm, group = site.treatment, color = site.treatment, shape = site.treatment, size = site.treatment)) +
  geom_vline(xintercept = 2014, linetype = "longdash", color = vadelma2) + # AF line
  stat_summary(fun = mean, geom = "point") +
  stat_summary(fun = mean, geom = "line", size = 0.5) +
  scale_color_manual(values = c("Shallow WT.Control" = vesi2, "Shallow WT.Treatment" = vesi1, "Deep WT.Control" = valo2, "Deep WT.Treatment" = valo1),
                     labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_shape_manual(values = c("Shallow WT.Control" = 0, "Shallow WT.Treatment" = 15, "Deep WT.Control" = 5, "Deep WT.Treatment" = 18),
                     labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_size_manual(values = c("Shallow WT.Control" = 1.5, "Shallow WT.Treatment" = 1.5, "Deep WT.Control" = 1.5, "Deep WT.Treatment" = 2),
                    labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
   labs(x = "Year", y = "Water table (cm)", color = "", shape = "", size = "") +
  scale_x_continuous(breaks = 2006:2023,  labels = ifelse(2006:2023 %% 2 == 0, 2006:2023, "")) +
  scale_y_continuous(breaks = seq(-80, 0, by = 10)) +
  coord_cartesian(xlim = c(2007, 2020), ylim = c(-85, 0)) +
  geom_point(x = 2009, y = wt_thinn_dry, shape = 6, size = 5, color = vadelma1) + # thinning point deep WT
  annotate("text", x = 2009, y = wt_thinn_dry, label = "thinning", hjust = 0.5, vjust = -1.2, size = 2.5, color = vadelma1) + # thinning text deep WT
  geom_point(x = 2010, y = wt_thinn_wet, shape = 6, size = 5, color = vadelma1) + # thinning point shallow WT
  annotate("text", x = 2010, y = wt_thinn_wet, label = "thinning", hjust = 0.5, vjust = -1.2, size = 2.5, color = vadelma1) + # thinning texts shallow WT
  geom_point(x = 2014, y = wt_af_dry, shape = 1, size = 5, color = vadelma1) + # AF point deep WT
  annotate("text", x = 2014, y = wt_af_dry, label = "AF", hjust = 1.2, vjust = -1.1, size = 2.5, color = vadelma1) + # AF text deep WT
  geom_point(x = 2014, y = wt_af_wet, shape = 1, size = 5, color = vadelma1) + # AF point shallow WT
  annotate("text", x = 2014, y = wt_af_wet, label = "AF", hjust = 1.2, vjust = -1.1, size = 2.5, color = vadelma1) + # AF text shallow WT
  theme_bw(base_size = 10) +
  theme(legend.position = "right", legend.text = element_text(size = 9), legend.margin = margin(0, 0, 0, 0, "cm"),
        axis.ticks.x = element_line())

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

Warning: Removed 3080 rows containing non-finite outside the scale range
(`stat_summary()`).
Removed 3080 rows containing non-finite outside the scale range
(`stat_summary()`).

Statistics:

WT_data |>
  filter(year > 2005 & year < 2021, !is.na(WT_cm)) |> 
  mutate(time = factor(ifelse(year < 2010, "pre", ifelse(year > 2013, "postAF", "postT")))) |> 
  group_by(site, treatment, time) |> 
  summarise(n = n(), mean_WT = mean(WT_cm), sd_WT = sd(WT_cm), .groups = "drop") |> 
  arrange(site, treatment, desc(time)) |> 
  group_by(site, time) |> 
  select(site, treatment, time, mean_WT) |> 
  pivot_wider(names_from = treatment, values_from = mean_WT) |> 
  mutate(difference = Treatment - Control)

# A tibble: 6 × 5
# Groups:   site, time [6]
  site       time   Control Treatment difference
  <fct>      <fct>    <dbl>     <dbl>      <dbl>
1 Shallow WT pre      -27.3     -23.9       3.37
2 Shallow WT postT    -28.7     -22.5       6.23
3 Shallow WT postAF   -33.9     -31.3       2.62
4 Deep WT    pre      -75.8     -72.1       3.67
5 Deep WT    postT    -68.7     -53.3      15.4 
6 Deep WT    postAF   -64.1     -54.6       9.47

Figure 3: Stand-level BA & BAI

Build the stand-level data:

# Read tree ring data:
bai <- read.csv2(file = "DATA d13C Muhos-Sievi.txt", header = FALSE, skip = 2, fill = TRUE, sep = "\t", as.is = TRUE) |> 
  filter(V11 == "RINGWIDTH" &  # Filter rows with "RINGWIDTH" in column 11
           V1 != "C1" & V1 != "C2" & V1 != "F1" & V1 != "F5") |> # Remove C1 & C2 (incomplete) & F1 & F5 (double)  
  rowwise() |> 
  mutate(
    g = list(as.numeric(c_across(V13:last_col()))), # Select ring observations from column 13 onwards
    g = list(na.omit(g)), # Remove NA values from the ring observations
    year = list((V4 - length(g) + 1):V4), # Calculate years
    rw = list(rev(g)) # Reverse ring width (bark to pit)
  ) |> 
  unnest(c(year, rw)) |> # Expand lists of years and ring widths into rows
  select(year, rw, tree = V1, site = V2) |> # Select the necessary columns and rename them
  mutate(
    id = paste(site, tree), # Create ID for grouping
    csum = ave(rw, id, FUN = cumsum), # Cumulative sum of ring widths, mm
    g = csum^2 * pi, # Calculate basal area, mm^2
    d.g = c(0, diff(g)), # Calculate BAI (basal area increment)
    d.g = ifelse(d.g < 0, 0, d.g) # Replace negative BAI values with 0
  ) |> 
  mutate(
    treat = substr(tree, 1, 1), # Extract the first letter of the tree code as treatment
    treatment = case_when(
      treat == "C" ~ "Control",
      treat == "F" ~ "Treatment",
      TRUE ~ NA_character_
    ),
    treatment = factor(treatment, levels = c("Control", "Treatment")),
    site = ifelse(site == "Muhos-Ansasaari", "Shallow WT", "Deep WT"),
    site = factor(site, levels = c("Shallow WT", "Deep WT"))
    ) |>
  na.omit() |> 
  group_by(site, treatment, year) |> 
  summarize(bai = mean(d.g, na.rm = TRUE), .groups = "drop")

# calculate growth proportions
bai_normalized <- bind_rows(
  bai |>
    filter(site == "Shallow WT", treatment == "Treatment", year >= 2009 & year <= 2010) |>  
    arrange(year) |> 
    mutate(total_growth = sum(bai, na.rm = TRUE),
           growth_proportion = bai / total_growth,
           relative_increase = cumsum(growth_proportion) 
           ) |>
    filter(relative_increase != 1) |> 
    complete(site = "Shallow WT", treatment = "Treatment", year = 2008:2009),
  bai |>
    filter(site == "Shallow WT", treatment == "Treatment", year >= 2011 & year <= 2019) |> 
    arrange(year) |> 
    mutate(total_growth = sum(bai, na.rm = TRUE),
           growth_proportion = bai / total_growth,
           relative_increase = cumsum(growth_proportion) 
           ) |> 
    filter(relative_increase != 1) |> 
    complete(site = "Shallow WT", treatment = "Treatment", year = 2010:2018),
  bai |>
    filter(site == "Shallow WT", treatment == "Control", year >= 2009 & year <= 2019) |> 
    arrange(year) |> 
    mutate(total_growth = sum(bai, na.rm = TRUE),
           growth_proportion = bai / total_growth,
           relative_increase = cumsum(growth_proportion) 
           ) |>
    filter(relative_increase != 1) |> 
    complete(site = "Shallow WT", treatment = "Control", year = 2008:2018),
  bai |>
    filter(site == "Deep WT", treatment == "Control", year >= 2007 & year <= 2011) |> 
    arrange(year) |> 
    mutate(total_growth = sum(bai, na.rm = TRUE),
           growth_proportion = bai / total_growth, 
           relative_increase = cumsum(growth_proportion) 
           ) |>
    filter(relative_increase != 1) |> 
    complete(site = "Deep WT", treatment = "Control", year = 2006:2010),
  bai |>
    filter(site == "Deep WT", treatment == "Treatment", year >= 2007 & year <= 2009) |> 
    arrange(year) |> 
    mutate(total_growth = sum(bai, na.rm = TRUE),
           growth_proportion = bai / total_growth, 
           relative_increase = cumsum(growth_proportion) 
           ) |>
    filter(relative_increase != 1) |> 
    complete(site = "Deep WT", treatment = "Treatment", year = 2006:2008),
  bai |> 
    filter(site == "Deep WT", treatment == "Treatment", year > 2009 & year <= 2011) |> 
    arrange(year) |> 
    mutate(total_growth = sum(bai, na.rm = TRUE),
           growth_proportion = bai / total_growth, 
           relative_increase = cumsum(growth_proportion) 
           ) |> 
    filter(relative_increase != 1) |> 
    complete(site = "Deep WT", treatment = "Treatment", year = 2009:2010),
  bai |> 
    filter(site == "Deep WT", year >= 2012 & year <= 2019) |>
    group_by(treatment) |>
    arrange(treatment, year) |> 
    mutate(total_growth = sum(bai, na.rm = TRUE),
           growth_proportion = bai / total_growth,
           relative_increase = cumsum(growth_proportion) 
           ) |> 
    filter(relative_increase != 1) |> 
    complete(site = "Deep WT", year = 2011:2018),
  bai |> 
    filter(year >= 2020 & year <= 2022) |> 
    group_by(site, treatment) |> 
    arrange(site, treatment, year) |> 
    complete(year = 2020:2024) |> 
    fill(bai, .direction = "updown") |> 
    mutate(bai = case_when(year == 2023 ~ mean(bai), 
                         year == 2024 ~ mean(bai), 
                         TRUE ~ bai)
           ) |> 
    mutate(total_growth = sum(bai, na.rm = TRUE),
           growth_proportion = bai / total_growth,
           relative_increase = cumsum(growth_proportion) 
           ) |>
    filter(year != 2024) |> 
    complete(year = 2019:2024)
) |> 
  ungroup() |> 
  replace_na(list(relative_increase = 0)) |> 
  select(site, treatment, year, relative_increase) |> 
  arrange(site, treatment, year)

# Read forest inventory data:
forest <- bind_rows(
  read.csv2("DATA stand inventory Muhos.csv", sep = ",", na="") |> 
    filter(tulostusryhma == 1) |> 
    mutate(site = "Shallow WT",
           treatment = factor(ifelse(koealatunnus == "1", "Treatment", ifelse(koealatunnus == "2", "Control", NA_character_)), levels = c("Control", "Treatment")), # Muhos: 1 = Fertilization, 2 = Control
           year = as.numeric(ajankohtavuosi),
           G = as.numeric(G),
           N = as.numeric(N)
    ) |>  
    select(site, treatment, year, G, N), 
  read.csv2("DATA stand inventory Sievi.csv", sep = ",", na="") |> 
    filter(tulostusryhma == 1) |> 
    mutate(site = "Deep WT",
           treatment = factor(ifelse(koealatunnus == "1", "Control", ifelse(koealatunnus == "2", "Treatment", NA_character_)), levels = c("Control", "Treatment")), # Sievi: 1 = Control, 2 = Fertilization
           year = as.numeric(ajankohtavuosi),
           G = as.numeric(G),
           N = as.numeric(N)
    ) |> 
  select(site, treatment, year, G, N)) |> 
  arrange(site, treatment, year) |> 
  na.omit()   

# Calculate scaled BA:
growth <- forest |>
  group_by(site, treatment) |>
  complete(year = 2006:2024) |>
  arrange(site, treatment, year) |> 
  fill(c(G, N), .direction = "down") |> 
  na.omit() |> 
  mutate(total_growth = case_when(
      site == "Shallow WT" & treatment == "Treatment" & year >= 2008 & year <= 2010 ~ G[year == 2010] - G,
      site == "Shallow WT" & treatment == "Treatment" & year >= 2010 & year <= 2019 ~ G[year == 2019] - G, # replace Ansasaari F 2011-2018
      site == "Shallow WT" & treatment == "Control" & year >= 2008 & year <= 2019 ~ G[year == 2019] - G,
      site == "Deep WT" & treatment == "Treatment" & year >= 2006 & year <= 2009 ~ G[year == 2009] - G,
      site == "Deep WT" & treatment == "Treatment" & year >= 2009 & year <= 2011 ~ G[year == 2011] - G, # replace Nevajärvi F 2010
      site == "Deep WT" & year >= 2011 & year <= 2019 ~ G[year == 2019] - G,
      site == "Deep WT" & treatment == "Control" & year >= 2006 & year <= 2011 ~ G[year == 2011] - G,
      year >= 2019 & year <= 2024 ~ G[year == 2024] - G,
      year == 2024 ~ 0,
      TRUE ~ NA_real_
      ) 
  ) |> 
  ungroup() |> 
  # Combine with growth proportion:  
  left_join(bai_normalized, by = c("site", "treatment", "year")) |> 
  mutate(scaled_g = G + (total_growth * relative_increase),
         treatment = factor(treatment, levels = c("Control", "Treatment"))
         ) |> 
  # Remove double:
  filter(!(year == 2019 & site == "Deep WT" & treatment == "Treatment" & round(relative_increase, 0) == 1)) 

# Replace Shallow WT Fertilized 2011-2018:
growth <- growth |> 
    left_join(bai |> 
              filter(site == "Shallow WT" & treatment == "Treatment" & year >= 2011 & year <= 2018) |> 
              arrange(desc(year)) |> 
              mutate(N = forest |>  filter(site == "Shallow WT", treatment == "Treatment", year == 2019) |>  pull(N),
                     ansasFbefore2018 = N * bai / 1e6,
                     ansasFbefore2018 = cumsum(ansasFbefore2018)
                     ) |> 
              select(site, treatment, year, ansasFbefore2018), 
              by = c("site", "treatment", "year") 
              ) |> 
  mutate(scaled_g = case_when(site == "Shallow WT" & treatment == "Treatment" & year >= 2011 & year <= 2018 ~ scaled_g - ansasFbefore2018, TRUE ~ scaled_g)
         ) 

# Add years before 2008 (Shallow WT) & 2006 (Deep WT):

g_ansasF2008 <- growth |> 
  filter(site == "Shallow WT", treatment == "Treatment", year == 2008) |> 
  pull(scaled_g)
g_ansasC2008 <- growth |> 
  filter(site == "Shallow WT", treatment == "Control", year == 2008) |> 
  pull(scaled_g)
g_nevajF2006 <- growth |> 
  filter(site == "Deep WT", treatment == "Treatment", year == 2006) |> 
  pull(scaled_g)
g_nevajC2006 <- growth |> 
  filter(site == "Deep WT", treatment == "Control", year == 2006) |> 
  pull(scaled_g)
g_nevajF2011 <- growth |> 
  filter(site == "Deep WT", treatment == "Treatment", year == 2011) |> 
  pull(scaled_g)

# Combine the data:
growth <- growth |> 
  complete(year = 2000:2024, site = site, treatment = treatment) |> 
  mutate(scaled_g = case_when(
    site == "Shallow WT" & treatment == "Treatment" & year <= 2008 ~ g_ansasF2008,
    site == "Shallow WT" & treatment == "Control" & year <= 2008 ~ g_ansasC2008,
    site == "Deep WT" & treatment == "Treatment" & year <= 2006 ~ g_nevajF2006,
    site == "Deep WT" & treatment == "Control" & year <= 2006 ~ g_nevajC2006,
    site == "Deep WT" & treatment == "Treatment" & year == 2010 ~ g_nevajF2011,
    TRUE ~ scaled_g 
    ) ) |> 
  left_join(bai |> 
              filter(site == "Shallow WT" & year < 2008 & year >= 2000) |>
              arrange(treatment, desc(year)) |> 
              group_by(treatment) |> 
              mutate(N = case_when(
                treatment == "Treatment" ~ forest |> filter(site == "Shallow WT", treatment == "Treatment", year == 2008) |> pull(N),
                treatment == "Control" ~ forest |> filter(site == "Shallow WT", treatment == "Control", year == 2008) |> pull(N),
                TRUE ~ NA_real_),
                before2008 = N * bai / 1e6,
                before2008 = cumsum(before2008)
                ) |> 
              ungroup() |> 
              select(site, treatment, year, before2008),
            by = c("site", "treatment", "year")
            ) |> 
  left_join(bai |> 
              filter(site == "Deep WT" & year < 2006 & year >= 2000) |> 
              arrange(treatment, desc(year)) |> 
              group_by(treatment) |> 
              mutate(N = case_when(
                treatment == "Treatment" ~ forest |> filter(site == "Deep WT", treatment == "Treatment", year == 2006) |> pull(N),
                treatment == "Control" ~ forest |> filter(site == "Deep WT", treatment == "Control", year == 2006) |> pull(N),
                TRUE ~ NA_real_), 
                before2006 = N * bai / 1e6,
                before2006 = cumsum(before2006)
                ) |> 
              ungroup() |> 
              select(site, treatment, year, before2006),
            by = c("site", "treatment", "year")) |> 
  left_join(bai |> 
              filter(site == "Deep WT" & treatment == "Treatment" & year == 2010) |>               mutate(N = forest |> filter(site == "Deep WT", treatment == "Treatment", year == 2019) |> pull(N),
                     nevajF2010 = N * bai / 1e6) |> 
              select(site, treatment, year, nevajF2010),
            by = c("site", "treatment", "year") ) |> 
  mutate(scaled_g = case_when(
    site == "Shallow WT" & year < 2008  ~ scaled_g - before2008,
    site == "Deep WT" & year < 2006  ~ scaled_g - before2006,
    site == "Deep WT" & treatment == "Treatment" & year == 2010 ~ scaled_g - nevajF2010,
    TRUE ~ scaled_g)) |> 
  select(site, treatment, year, scaled_g)

# Calculate scaled BAI:
bai_scaled <- growth |> 
  arrange(site, treatment, year) |>
  group_by(site, treatment) |>
  mutate(bai = scaled_g - lag(scaled_g) ) |> 
  ungroup() |> 
  mutate(bai = case_when(
    site == "Shallow WT" & treatment == "Treatment" & year == 2011 ~ 0.22245082, # value from bai df
    site == "Deep WT" & treatment == "Treatment" & year == 2010 ~ 0.23707775, # value from bai df
    TRUE ~ bai
  )) |> 
  select(site, treatment, year, bai)

# For plotting: BA & BAI thinning & AF
g_thinn_dry1 <- growth |> 
  filter(site == "Deep WT", treatment == "Treatment", year == "2009") |> 
  pull(scaled_g)
g_thinn_dry2 <- growth |> 
  filter(site == "Deep WT", treatment == "Treatment", year == "2010") |> 
  pull(scaled_g)
g_thinn_wet1 <- growth |> 
  filter(site == "Shallow WT", treatment == "Treatment", year == "2010") |> 
  pull(scaled_g)
g_thinn_wet2 <- growth |> 
  filter(site == "Shallow WT", treatment == "Treatment", year == "2011") |> 
  pull(scaled_g)
g_af_dry <- growth |> 
  filter(site == "Deep WT", treatment == "Treatment", year == "2014") |> 
  pull(scaled_g)
g_af_wet <- growth |> 
  filter(site == "Shallow WT", treatment == "Treatment", year == "2014") |> 
  pull(scaled_g)
dg_thinn_dry <- bai_scaled |> 
  filter(site == "Deep WT", treatment == "Treatment", year == "2009") |> 
  pull(bai)
dg_thinn_wet <- bai_scaled |> 
  filter(site == "Shallow WT", treatment == "Treatment", year == "2010") |> 
  pull(bai)
dg_af_dry <- bai_scaled |> 
  filter(site == "Deep WT", treatment == "Treatment", year == "2014") |> 
  pull(bai)
dg_af_wet <- bai_scaled |> 
  filter(site == "Shallow WT", treatment == "Treatment", year == "2014") |> 
  pull(bai)

Thinning % according to BA:

bind_rows(
  growth |>
    filter(year == 2010 | year == 2011) |>
    filter(site == "Shallow WT" & treatment == "Treatment") |>
    pivot_wider(id_cols = c(site), names_from = year, values_from = scaled_g) |>
    mutate(change_perc = (`2011` - `2010`) / `2010`),    growth |>
    filter(year == 2009 | year == 2010) |>
    filter(site == "Deep WT" & treatment == "Treatment") |>
    pivot_wider(id_cols = c(site), names_from = year, values_from = scaled_g) |>
    mutate(change_perc = (`2010` - `2009`) / `2009`) ) |>
  select(site, change_perc)

# A tibble: 2 × 2
  site       change_perc
  <chr>            <dbl>
1 Shallow WT      -0.177
2 Deep WT         -0.291

Post-AF BA differences between Control & Treatment:

growth |>
  filter(year == 2014 | year == 2024) |>
  pivot_wider(id_cols = c(site, treatment), names_from = year, values_from = scaled_g) |>
  mutate(post_af_growth = `2024` - `2014`) |>
  pivot_wider(id_cols = site, names_from = treatment, values_from = post_af_growth) |>
  mutate(growth_diff = Treatment - Control)

# A tibble: 2 × 4
  site       Control Treatment growth_diff
  <chr>        <dbl>     <dbl>       <dbl>
1 Deep WT       2.82      3.76       0.940
2 Shallow WT    3.60      5.77       2.17 

Plot:

# Plot BA shallow WT 
p_3a_BA <- growth |> 
  filter(year >= 2005 & site == "Shallow WT") |> 
  mutate(thinning = case_when(
    site == "Shallow WT" & treatment == "Treatment" & year <= 2010 ~ "Shallow WT.Treatment.pre",
    site == "Shallow WT" & treatment == "Treatment" & year >= 2011 ~ "Shallow WT.Treatment.post",
    site == "Shallow WT" & treatment == "Control" ~ "Shallow WT.Control",
    TRUE ~ NA_character_
    ) ) |>
  ggplot(aes(x = year, group = thinning, color = treatment)) + 
  geom_vline(xintercept = 2014, linetype = "longdash", color = vadelma2) + # AF line
  stat_summary(aes(y = scaled_g), geom = "line") +
  annotate("segment", x = 2010, y = g_thinn_wet1, xend = 2011, yend = g_thinn_wet2, linetype = "dotted", color = vesi1) +  # post-thinning line
  stat_summary(aes(y = scaled_g, shape = treatment), size = 2, geom = "point") +
  geom_point(data = forest |> 
               filter(site == "Shallow WT") |> 
               mutate(site_treatment = factor(interaction(site, treatment), levels = c("Shallow WT.Treatment", "Shallow WT.Control", "Deep WT.Treatment", "Deep WT.Control")),
                      thinning = case_when(
                        site_treatment == "Shallow WT.Treatment" & year <= 2010 ~ paste0(site_treatment, ".pre"),
                        site_treatment == "Shallow WT.Treatment" & year >= 2012 ~ paste0(site_treatment, ".post"),
                        site_treatment == "Shallow WT.Control" ~ paste0(site_treatment),
                        TRUE ~ NA_character_
                      ) ),
             aes(x = year, y = G), shape = 0, size = 3) +
  scale_color_manual(values = c("Control" = vesi2, "Treatment" = vesi1), labels = c("Control", "Treatment")) +
  scale_shape_manual(values = c("Control" = 0, "Treatment" = 15), labels = c("Control", "Treatment")) +
  labs(title = "a)", x = "Year", y = expression(atop("Basal area", (m^2 ~ ha^{-1}))), color = "Shallow WT", shape = "Shallow WT") +
  scale_x_continuous(breaks = seq(2005, 2024, by = 1), labels = function(x) ifelse(x %% 5 == 0, x, "")) +
  coord_cartesian(x = c(2005, 2024), y = c(10, 25)) +
scale_x_continuous(breaks = seq(2005, 2024, by = 1), labels = function(x) ifelse(x %% 5 == 0, x, "")) +
  geom_point(x = 2014, y = g_af_wet, shape = 1, size = 5, color = vadelma1) + # AF point
  annotate("text", x = 2014, y = g_af_wet, label = "AF", hjust = 1.2, vjust = 2.4, size = 3, color = vadelma1) +  # AF text
  geom_point(x = 2010, y = g_thinn_wet1, shape = 6, size = 4, color = vadelma1) + # thinning point
  annotate("text", x = 2010, y = g_thinn_wet1, label = "thinning", hjust = 0.5, vjust = -1, size = 2.8, color = vadelma1) + # thinning text
  theme_bw(base_size = 10) +
  theme(legend.position = "top", legend.title = element_text(size = 10), legend.text  = element_text(size = 10))

Scale for x is already present.
Adding another scale for x, which will replace the existing scale.

# Plot BA deep WT
p_3b_BA <- growth |> 
  filter(year >= 2005 & site == "Deep WT") |> 
  mutate(thinning = case_when(
    site == "Deep WT" & treatment == "Treatment" & year <= 2009 ~ "Deep WT.Treatment.pre",
    site == "Deep WT" & treatment == "Treatment" & year >= 2010 ~ "Deep WT.Treatment.post",
    site == "Deep WT" & treatment == "Control" ~ "Deep WT.Control",
    TRUE ~ NA_character_
    ) ) |>
  ggplot(aes(x = year, group = thinning, color = treatment)) + 
  geom_vline(xintercept = 2014, linetype = "longdash", color = vadelma2) + # AF line
  stat_summary(aes(y = scaled_g), geom = "line") +
  annotate("segment", x = 2009, y = g_thinn_dry1, xend = 2010, yend = g_thinn_dry2, linetype = "dotted", color = valo1) +  # post-thinning line
  stat_summary(aes(y = scaled_g, shape = treatment, size = treatment), geom = "point") +
  geom_point(data = forest |> 
               filter(site == "Deep WT") |> 
               mutate(site_treatment = factor(interaction(site, treatment), levels = c("Shallow WT.Treatment", "Shallow WT.Control", "Deep WT.Treatment", "Deep WT.Control")),
                      thinning = case_when(
                        site_treatment == "Deep WT.Treatment" & year <= 2010 ~ paste0(site_treatment, ".pre"),
                        site_treatment == "Deep WT.Treatment" & year >= 2012 ~ paste0(site_treatment, ".post"),
                        site_treatment == "Deep WT.Control" ~ paste0(site_treatment),
                        TRUE ~ NA_character_
                      ) ),
             aes(x = year, y = G), shape = 5, size = 3) + 
   scale_color_manual(values = c("Control" = valo2, "Treatment" = valo1), labels = c("Control", "Treatment")) +
  scale_shape_manual(values = c("Control" = 5, "Treatment" = 18), labels = c("Control", "Treatment")) +
  scale_size_manual(values = c("Control" = 2, "Treatment" = 2.5), labels = c("Control", "Treatment")) +
  labs(title = "b)", x = "Year", y = expression(atop("Basal area", (m^2 ~ ha^{-1}))), color = "Deep WT", shape = "Deep WT", size = "Deep WT") +
  geom_point(x = 2014, y = g_af_dry, shape = 1, size = 5, color = vadelma1) + # AF point
  annotate("text", x = 2014, y = g_af_dry, label = "AF", hjust = 1.2, vjust = -1, size = 2.8, color = vadelma1) +  # AF text
  geom_point(x = 2009, y = g_thinn_dry1, shape = 6, size = 4, color = vadelma1) + # thinning point
  annotate("text", x = 2009, y = g_thinn_dry1, label = "thinning", hjust = 0.5, vjust = 2.65, size = 2.8, color = vadelma1) +  # thinning text
  scale_x_continuous(breaks = seq(2005, 2024, by = 1), labels = function(x) ifelse(x %% 5 == 0, x, "")) +
  theme_bw(base_size = 10) +
  theme(legend.position = "top",  legend.title = element_text(size = 10), legend.text  = element_text(size = 10))

# Plot BAI shallow WT
p_3a_BAI <- bai_scaled |> 
  filter(site == "Shallow WT" & year >= 2005) |> 
  ggplot(aes(x = year, group = treatment, color = treatment)) +
  geom_vline(xintercept = 2014, linetype = "longdash", color = vadelma2) + #AF line 
  stat_summary(aes(y = bai), geom = "line") +
  stat_summary(aes(y = bai, shape = treatment), size = 2, geom = "point") +
  scale_color_manual(values = c("Control" = vesi2, "Treatment" = vesi1), labels = c("Shallow WT, Control", "Shallow WT, Treatment")) +
  scale_shape_manual(values = c("Control" = 0, "Treatment" = 15), labels = c("Shallow WT, Control", "Shallow WT, Treatment")) +
  labs(title = "", x = "Year", y = expression(atop("Basal area increment", (m^2 ~ ha^{-1} ~ yr^{-1}))), color = "", shape = "") +
  geom_point(x = 2010, y = dg_thinn_wet, shape = 6, size = 4, color = vadelma1) + # thinning point
  annotate("text", x = 2010, y = dg_thinn_wet, label = "thinning", hjust = 0.38, vjust = -1, size = 2.8, color = vadelma1) + # thinning text
  geom_point(x = 2014, y = dg_af_wet, shape = 1, size = 5, color = vadelma1) + # AF point
  annotate("text", x = 2014, y = dg_af_wet, label = "AF", hjust = 1.2, vjust = 2.4, size = 2.8, color = vadelma1) + # AF text
  scale_x_continuous(breaks = seq(2005, 2024, by = 1), labels = function(x) ifelse(x %% 5 == 0, x, "")) +
  coord_cartesian(x = c(2005, 2024), y = c(0, 1.14)) +
  theme_bw(base_size = 10) +
  theme(legend.position = "top")

# Plot BAI deep WT
p_3b_BAI <- bai_scaled |> 
  filter(site == "Deep WT" & year >= 2005) |> 
  ggplot(aes(x = year, group = treatment, color = treatment)) +
  geom_vline(xintercept = 2014, linetype = "longdash", color = vadelma2) + # AF line
  stat_summary(aes(y = bai), geom = "line") +
  stat_summary(aes(y = bai, shape = treatment, size = treatment), geom = "point") +
  scale_color_manual(values = c("Treatment" = valo1, "Control" = valo2), labels = c("Deep WT, Treatment", "Deep WT, Control")) +
  scale_shape_manual(values = c("Treatment" = 18, "Control" = 5), labels = c("Deep WT, Treatment", "Deep WT, Control")) +
  scale_size_manual(values = c("Treatment" = 2.5, "Control" = 2), labels = c("Deep WT, Treatment", "Deep WT, Control")) +
  labs(title = "", x = "Year", y = expression(atop("Basal area increment", (m^2 ~ ha^{-1} ~ yr^{-1}))), color = "", shape = "", size = "") +
  geom_point(x = 2009, y = dg_thinn_dry, shape = 6, size = 4, color = vadelma1) + # thinning point
  annotate("text", x = 2009, y = dg_thinn_dry, label = "thinning", hjust = 0.65, vjust = -1, size = 2.8, color = vadelma1) + # thinning text
  geom_point(x = 2014, y = dg_af_dry, shape = 1, size = 5, color = vadelma1) + # AF point
  annotate("text", x = 2014, y = dg_af_dry, label = "AF", hjust = 1.2, vjust = 2.4, size = 2.8, color = vadelma1) + # AF text
  coord_cartesian(x = c(2005, 2024), y = c(0, 1.14)) +
  scale_x_continuous(breaks = seq(2005, 2024, by = 1), labels = function(x) ifelse(x %% 5 == 0, x, "")) +
  theme_bw(base_size = 10) +
  theme(legend.position = "top")

legend_shallow <- get_legend(p_3a_BA)

No summary function supplied, defaulting to `mean_se()`
No summary function supplied, defaulting to `mean_se()`

legend_deep <- get_legend(p_3b_BA)

No summary function supplied, defaulting to `mean_se()`
No summary function supplied, defaulting to `mean_se()`

p_growth_shallow <- plot_grid(
  plot_grid(legend_shallow), 
  p_3a_BA + theme(legend.position = "none"),
  p_3a_BAI + theme(legend.position = "none"),
  ncol = 1, 
  rel_heights = c(0.1, 1, 1)
  ) 

No summary function supplied, defaulting to `mean_se()`
No summary function supplied, defaulting to `mean_se()`
No summary function supplied, defaulting to `mean_se()`
No summary function supplied, defaulting to `mean_se()`

p_growth_deep <- plot_grid(
  plot_grid(legend_deep), 
  p_3b_BA + theme(legend.position = "none"),
  p_3b_BAI + theme(legend.position = "none"),
  ncol = 1, 
  rel_heights = c(0.1, 1, 1)
) 

No summary function supplied, defaulting to `mean_se()`
No summary function supplied, defaulting to `mean_se()`
No summary function supplied, defaulting to `mean_se()`
No summary function supplied, defaulting to `mean_se()`

p_growth_shallow + p_growth_deep

Figure 4: BAI lme-models

Create models:

# Data for shallow WT
bai_data_ansas <- bai_scaled |> 
  filter(site == "Shallow WT" & year >= 2001) |> 
  mutate(
    site = ifelse(site == "Shallow WT", "Shallow WT", "Deep WT"),
    yr_treat = factor(ifelse(treatment == "Treatment", pmax(year - 2014, -1), -1)),
    thinned = factor(treatment == "Treatment" & year > 2010),
    year = as.numeric(year),
    treatment = relevel(treatment, ref = "Control")
  ) |> 
  na.omit()

# Model for shallow WT
bai_model_ansas <- lme(bai ~ treatment + thinned + yr_treat, random = ~1|year, data = bai_data_ansas)

# Data for deep WT
bai_data_nevaj <- bai_scaled |> 
  filter(site == "Deep WT" & year >= 2001) |> 
  mutate(
    site = ifelse(site == "Deep WT", "Deep WT", "Shallow WT"),
    yr_treat = factor(ifelse(treatment == "Treatment", pmax(year - 2014, -1), -1)),
    thinned = factor(treatment == "Treatment" & year > 2009),
    year = as.numeric(year),
    treatment = relevel(treatment, ref = "Control")
  ) |> 
  na.omit()


# Model for deep WT
bai_model_nevaj <- lme(bai ~ treatment + thinned + yr_treat, random = ~1|year, data = bai_data_nevaj)

Print data for supplementary data:

# Data for shallow WT
as.data.frame(round(summary(bai_model_ansas)$tTable, 3)) |> 
  mutate(p.adj.signif = sapply(`p-value`, function(p) {
    if (p < 0.001) { return("***") } 
    else if (p < 0.01) { return("**")} 
    else if (p < 0.05) { return("*") } 
    else if (p < 0.10) { return(".")}
    else { return("ns") } }))

                    Value Std.Error DF t-value p-value p.adj.signif
(Intercept)         0.540     0.042 23  12.751   0.000          ***
treatmentTreatment  0.076     0.074 11   1.035   0.323           ns
thinnedTRUE        -0.328     0.134 11  -2.438   0.033            *
yr_treat0           0.055     0.236 11   0.234   0.819           ns
yr_treat1           0.121     0.236 11   0.514   0.618           ns
yr_treat2           0.146     0.236 11   0.617   0.550           ns
yr_treat3           0.404     0.236 11   1.711   0.115           ns
yr_treat4           0.399     0.236 11   1.694   0.118           ns
yr_treat5           0.331     0.236 11   1.403   0.188           ns
yr_treat6           0.277     0.236 11   1.175   0.265           ns
yr_treat7           0.334     0.236 11   1.414   0.185           ns
yr_treat8           0.428     0.236 11   1.815   0.097            .
yr_treat9           0.379     0.236 11   1.607   0.136           ns
yr_treat10          0.379     0.236 11   1.607   0.136           ns

# Data for deep WT
as.data.frame(round(summary(bai_model_nevaj)$tTable, 3)) |> 
  mutate(p.adj.signif = sapply(`p-value`, function(p) {
    if (p < 0.001) { return("***") } 
    else if (p < 0.01) { return("**")} 
    else if (p < 0.05) { return("*") } 
    else if (p < 0.10) { return(".") }
    else { return("ns") } }))

                    Value Std.Error DF t-value p-value p.adj.signif
(Intercept)         0.427     0.037 23  11.693   0.000          ***
treatmentTreatment -0.062     0.068 11  -0.925   0.375           ns
thinnedTRUE        -0.013     0.107 11  -0.122   0.905           ns
yr_treat0           0.154     0.199 11   0.776   0.454           ns
yr_treat1           0.169     0.199 11   0.850   0.414           ns
yr_treat2           0.089     0.199 11   0.447   0.663           ns
yr_treat3          -0.021     0.199 11  -0.108   0.916           ns
yr_treat4          -0.030     0.199 11  -0.150   0.883           ns
yr_treat5          -0.029     0.199 11  -0.144   0.888           ns
yr_treat6          -0.112     0.199 11  -0.562   0.586           ns
yr_treat7           0.033     0.199 11   0.167   0.870           ns
yr_treat8           0.156     0.199 11   0.783   0.450           ns
yr_treat9           0.078     0.199 11   0.391   0.704           ns
yr_treat10          0.078     0.199 11   0.391   0.704           ns

# Variance of random year effect, shallow WT
as.numeric(getVarCov(bai_model_ansas))

[1] 0.007326456

# Residual variance, shallow WT
bai_model_ansas$sigma^2

[1] 0.0356448

# Variance of random year effect, deep WT
bai_model_nevaj$sigma^2

[1] 0.02841263

# Residual variance, deep WT
as.numeric(getVarCov(bai_model_nevaj))

[1] 0.003592523

Plot:

p_bai_lme_shal <- as.data.frame(summary(bai_model_ansas)$tTable) |> 
  rownames_to_column("Term") |> 
  filter(grepl("^yr_treat", Term)) |> 
  mutate(Year = as.numeric(gsub("yr_treat", "", Term)),
         Lower = Value - 1.96 * Std.Error,
         Upper = Value + 1.96 * Std.Error,
         Site = "Shallow WT") |> 
  ggplot(aes(x = Year)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = kallio1) +
  geom_point(aes(y = Value, color = Site, shape = Site, size = Site)) +
  geom_line(aes(y = Value), color = vesi1) +
  geom_line(aes(y = Lower), color = vesi3) +
  geom_line(aes(y = Upper), color = vesi3) +
  scale_color_manual(values = c("Shallow WT" = vesi1)) +
  scale_shape_manual(values = c("Shallow WT" = 15)) +
  scale_size_manual(values = c("Shallow WT" = 2)) +
  labs(title = "a)", x = "Year after AF", y = expression(atop(BAI~treatment~effect, "(m"^2*" ha"^-1*" year"^-1*")")), color = "", shape = "", size = "") +
  scale_x_continuous(breaks = 0:10) +
  coord_cartesian(y = c(-0.5, 1)) +
  annotate("text", x = 8, y = 0.54, label = ".", hjust = 0.5, vjust = 0, size = 10, color = vadelma1)+
  theme_bw(base_size = 10) +
  theme(legend.position = "top", legend.margin = margin(0,0,0,0, "cm"), legend.text  = element_text(size = 10), plot.margin = margin(0,0.1,0,0, "cm")) 

p_bai_lme_deep <- as.data.frame(summary(bai_model_nevaj)$tTable) |> 
  rownames_to_column("Term") |> 
  filter(grepl("^yr_treat", Term)) |> 
  mutate(Year = as.numeric(gsub("yr_treat", "", Term)),
         Lower = Value - 1.96 * Std.Error,
         Upper = Value + 1.96 * Std.Error,
         Site = "Deep WT") |> 
  ggplot(aes(x = Year)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = kallio1) +
  geom_point(aes(y = Value, color = Site, shape = Site, size = Site)) +
  geom_line(aes(y = Value), color = valo1) +
  geom_line(aes(y = Lower), color = valo3) +
  geom_line(aes(y = Upper), color = valo3) +
  scale_color_manual(values = c("Deep WT" = valo1)) +
  scale_shape_manual(values = c("Deep WT" = 18)) +
  scale_size_manual(values = c("Deep WT" = 3)) +
  labs(title = "b)", x = "Year after AF", y = "", color = "", shape = "", size = "") +
  scale_x_continuous(breaks = 0:10) +
  coord_cartesian(y = c(-0.5, 1)) +
  theme_bw(base_size = 10) +
  theme(legend.position = "top", legend.margin = margin(0,0,0,0, "cm"), legend.text  = element_text(size = 10), plot.margin = margin(0,0.1,0,0, "cm"))

p_bai_lme_shal + p_bai_lme_deep

Figure 5: BAI/WT correlation

Create data:

# Treatment plot with shallow WT 
WT_BAI_ansasF_data <- WT_summary |>
  filter(site == "Shallow WT", treatment == "Treatment") |> 
  right_join(bai |> filter(site == "Shallow WT", treatment == "Treatment") |> select(site, year, treatment, bai), by = join_by(site == site, year == year, treatment == treatment)) |> 
  na.omit()

WT_BAI_ansasF_model <- lm(bai ~ WT_cm, data = WT_BAI_ansasF_data)

WT_BAI_ansasF_intercept <- round(coef(WT_BAI_ansasF_model)[1], 2)
WT_BAI_ansasF_slope <- round(coef(WT_BAI_ansasF_model)[2], 2)
WT_BAI_ansasF_r2 <- round(summary(WT_BAI_ansasF_model)$r.squared, 3)

WT_BAI_ansasF_text <- paste0("Treatment with shallow WT:\ny = ", WT_BAI_ansasF_intercept, ifelse(WT_BAI_ansasF_slope >= 0, " + ", " - "), abs(WT_BAI_ansasF_slope), "x\nR² = ", WT_BAI_ansasF_r2)

# Control plot with shallow WT
WT_BAI_ansasC_data <- WT_summary |>
  filter(site == "Shallow WT", treatment == "Control") |> 
  right_join(bai |> filter(site == "Shallow WT", treatment == "Control") |> select(site, year, treatment, bai), by = join_by(site == site, year == year, treatment == treatment)) |> 
  na.omit()

WT_BAI_ansasC_model <- lm(bai ~ WT_cm, data = WT_BAI_ansasC_data)

WT_BAI_ansasC_intercept <- round(coef(WT_BAI_ansasC_model)[1], 2)
WT_BAI_ansasC_slope <- round(coef(WT_BAI_ansasC_model)[2], 2)
WT_BAI_ansasC_r2 <- round(summary(WT_BAI_ansasC_model)$r.squared, 3)

WT_BAI_ansasC_text <- paste0("Control with shallow WT:\ny = ", WT_BAI_ansasC_intercept, ifelse(WT_BAI_ansasC_slope >= 0, " + ", " - "), abs(WT_BAI_ansasC_slope), "x\nR² = ", WT_BAI_ansasC_r2)

# Treatment plot with deep WT
WT_BAI_nevajF_data <- WT_summary |>
  filter(site == "Deep WT" & treatment == "Treatment") |> 
  right_join(bai |> filter(site == "Deep WT" & treatment == "Treatment") |> select(site, year, treatment, bai), by = join_by(site == site, year == year, treatment == treatment)) |> 
  na.omit()

WT_BAI_nevajF_model <- lm(bai ~ WT_cm, data = WT_BAI_nevajF_data)

WT_BAI_nevajF_intercept <- round(coef(WT_BAI_nevajF_model)[1], 2)
WT_BAI_nevajF_slope <- round(coef(WT_BAI_nevajF_model)[2], 2)
WT_BAI_nevajF_r2 <- round(summary(WT_BAI_nevajF_model)$r.squared, 3)

WT_BAI_nevajF_text <- paste0("Treatment with deep WT:\ny = ", WT_BAI_nevajF_intercept, ifelse(WT_BAI_nevajF_slope >= 0, " + ", " - "), abs(WT_BAI_nevajF_slope), "x\nR² = ", WT_BAI_nevajF_r2)

# Control plot with deept WT
WT_BAI_nevajC_data <- WT_summary |>
  filter(site == "Deep WT" & treatment == "Control") |> 
  right_join(bai |> filter(site == "Deep WT" & treatment == "Control") |> select(site, year, treatment, bai), by = join_by(site == site, year == year, treatment == treatment)) |> 
  na.omit()

WT_BAI_nevajC_model <- lm(bai ~ WT_cm, data = WT_BAI_nevajC_data)

WT_BAI_nevajC_intercept <- round(coef(WT_BAI_nevajC_model)[1], 2)
WT_BAI_nevajC_slope <- round(coef(WT_BAI_nevajC_model)[2], 3)
WT_BAI_nevajC_r2 <- round(summary(WT_BAI_nevajC_model)$r.squared, 3)

WT_BAI_nevajC_text <- paste0("Control with deep WT:\ny = ", WT_BAI_nevajC_intercept, ifelse(WT_BAI_nevajC_slope >= 0, " + ", " - "), abs(WT_BAI_nevajC_slope), "x\nR² = ", WT_BAI_nevajC_r2)

r & p-values:

# Treatment plot with shallow WT
round(cor(WT_BAI_ansasF_data$WT_cm, WT_BAI_ansasF_data$bai, use = "complete.obs"), 3)

[1] -0.653

round(cor.test(WT_BAI_ansasF_data$WT_cm, WT_BAI_ansasF_data$bai)$p.value, 3)

[1] 0.008

# Control plot with shallow WT
round(cor(WT_BAI_ansasC_data$WT_cm, WT_BAI_ansasC_data$bai, use = "complete.obs"), 3)

[1] 0.368

round(cor.test(WT_BAI_ansasC_data$WT_cm, WT_BAI_ansasC_data$bai)$p.value, 3)

[1] 0.177

# Treatment plot with deep WT
round(cor(WT_BAI_nevajF_data$WT_cm, WT_BAI_nevajF_data$bai, use = "complete.obs"), 3)

[1] 0.53

round(cor.test(WT_BAI_nevajF_data$WT_cm, WT_BAI_nevajF_data$bai)$p.value, 3)

[1] 0.035

# Control plot with deept WT
round(cor(WT_BAI_nevajC_data$WT_cm, WT_BAI_nevajC_data$bai, use = "complete.obs"), 3)

[1] 0.005

round(cor.test(WT_BAI_nevajC_data$WT_cm, WT_BAI_nevajC_data$bai)$p.value, 3)

[1] 0.986

Plot:

WT_summary |> 
  right_join(bai |> select(site, year, treatment, bai), by = join_by(site == site, year == year, treatment == treatment)) |> 
  na.omit() |> 
  mutate(site_treatment = factor(interaction(site, treatment), levels = c("Shallow WT.Control", "Shallow WT.Treatment", "Deep WT.Control", "Deep WT.Treatment"))) |> 
  ggplot(aes(x = WT_cm, y = bai)) +
  geom_point(aes(group = site_treatment, color = site_treatment, shape = site_treatment, size = site_treatment)) +
  geom_smooth(aes(group = site_treatment, color = site_treatment, fill = site_treatment), method = "lm", se = TRUE) +
  geom_smooth(aes(group = 1), method = "loess", se = TRUE, color = vadelma1, fill = NA) +
  scale_color_manual(values = c("Shallow WT.Treatment" = vesi1, "Shallow WT.Control" = vesi2, "Deep WT.Treatment" = valo1, "Deep WT.Control" = valo2),
                     labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_shape_manual(values = c("Shallow WT.Treatment" = 15, "Shallow WT.Control" = 0, "Deep WT.Treatment" = 18, "Deep WT.Control" = 5),
                     labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_size_manual(values = c("Shallow WT.Treatment" = 2, "Shallow WT.Control" = 2, "Deep WT.Treatment" = 3, "Deep WT.Control" = 2),
                    labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_fill_manual(values = c("Shallow WT.Treatment" = vesi1, "Shallow WT.Control" = vesi3, "Deep WT.Treatment" = valo1, "Deep WT.Control" = valo3),
                    labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  labs(title = "", x = "Water Table (cm)", y = expression("BAI (mm"^2*" tree"^-1*" year"^-1*")"), color = "", shape = "", size = "", fill = "") +
  coord_cartesian(x = c(-83, -18), y = c(0, 800)) +
  annotate("text", x = -78, y = 650, label = WT_BAI_nevajC_text, hjust = 0, vjust = 0, size = 3, color = valo2) +
  annotate("text", x = -78, y = 50, label = WT_BAI_nevajF_text, hjust = 0, vjust = 0, size = 3, color = valo1) +  
  annotate("text", x = -40, y = 50, label = WT_BAI_ansasC_text, hjust = 0, vjust = 0, size = 3, color = vesi2) +
  annotate("text", x = -40, y = 630, label = WT_BAI_ansasF_text, hjust = 0, vjust = 0, size = 3, color = vesi1) +
  theme_bw(base_size = 10) +
  theme(legend.position = "right", legend.text  = element_text(size = 9))

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

Figure 6. δ¹³C timeline

Read data:

d13c <- read_excel("d13c_SuoHitu.xlsx", sheet = "R", na = "") |>
  mutate(
    treatment = factor(treatment, levels = c("Control", "Fertilization"), labels = c("Control", "Treatment")),
    site = factor(site, levels = c("Ansasaari", "Nevajärvi"), labels = c("Shallow WT", "Deep WT"))
  )

d13c_thinn_wet <- d13c |> filter(site == "Shallow WT", treatment == "Treatment", year == 2010) |> pull(d13c)
d13c_af_wet <- d13c |> filter(site == "Shallow WT", treatment == "Treatment", year == 2014) |> pull(d13c)
d13c_thinn_dry <- d13c |> filter(site == "Deep WT", treatment == "Treatment", year == 2009) |> pull(d13c)
d13c_af_dry <- d13c |> filter(site == "Deep WT", treatment == "Treatment", year == 2014) |> pull(d13c)

Draw plot:

p_d13c_shal <- d13c |> 
  filter(year >= 2005 & site == "Shallow WT") |> 
  ggplot(aes(x = year, y = d13c, group = treatment, color = treatment, shape = treatment)) +
  geom_vline(xintercept = 2010, linetype = "dotted", color = vadelma2) +
  geom_vline(xintercept = 2014, linetype = "longdash", color = vadelma2) +
  stat_summary(fun = mean, geom = "point", size = 2) +
  stat_summary(fun = mean, geom = "line") +
  scale_size_continuous(range = c(0.7,1), guide = "none") +
  labs(title = "a)", x = "Year", y = expression(italic(delta)^13*C~("\u2030")), color = "Shallow WT", shape = "Shallow WT") +
  scale_color_manual(values = c("Control" = vesi2, "Treatment" = vesi1), labels = c("Control", "Treatment")) +
  scale_shape_manual(values = c("Control" = 0, "Treatment" = 15), labels = c("Control", "Treatment")) +
  coord_cartesian(y = c(-26, -22.3)) +
  scale_x_continuous(breaks = seq(2005, 2022, by = 1), labels = function(x) ifelse(x %% 5 == 0, x, "")) +
  scale_y_continuous(breaks = seq(-26, -22, by = 1)) +
  geom_point(x = 2010, y = d13c_thinn_wet, shape = 6, size = 3.5, color = vadelma1) +
  annotate("text", x = 2010, y = d13c_thinn_wet, label = "thinning", hjust = 0.5, vjust = 2, size = 2.8, color = vadelma1) +
  geom_point(x = 2014, y = d13c_af_wet, shape = 1, size = 3.5, color = vadelma1) + 
  annotate("text", x = 2014, y = d13c_af_wet, label = "AF", hjust = 0.5, vjust = 2, size = 2.8, color = vadelma1) + 
  theme_bw(base_size = 10) +
  theme(legend.position = "top", legend.margin = margin(0,0,0,0, "cm"), legend.title = element_text(size = 10), legend.text  = element_text(size = 10),
        axis.title.x = element_blank(), axis.ticks.x = element_line(), axis.ticks.length = unit(0.1, "cm"), 
        panel.grid.major.y = element_line(color = kallio2, linewidth = 0.25, linetype = "dotted"),
        plot.margin = margin(0, 1, 2, 1, "mm"))

p_d13c_deep <- d13c |> 
  filter(year >= 2005 & site == "Deep WT") |> 
  ggplot(aes(x = year, y = d13c, group = treatment, color = treatment, shape = treatment)) +
  geom_vline(xintercept = 2009, linetype = "dotted", color = vadelma2) +   # thinning
  geom_vline(xintercept = 2014, linetype = "longdash", color = vadelma2) + # fertilization
  stat_summary(aes(size = treatment), fun = mean, geom = "point") +
  stat_summary(fun = mean, geom = "line") +
  stat_summary(fun.data = mean_se, geom = "errorbar") +
  labs(title = "b)", x = "Year", y = expression(italic(delta)^13*C~("\u2030")), color = "Deep WT", shape = "Deep WT", size = "Deep WT") +
  scale_color_manual(values = c("Control" = valo2, "Treatment" = valo1), labels = c("Control", "Treatment")) +
  scale_shape_manual(values = c("Control" = 5, "Treatment" = 18), labels = c("Control", "Treatment")) +
  scale_size_manual(values = c("Control" = 2, "Treatment" = 2.5), labels = c("Control", "Treatment")) +
  coord_cartesian(y = c(-26, -22.3)) +
  scale_x_continuous(breaks = seq(2005, 2022, by = 1), labels = function(x) ifelse(x %% 5 == 0, x, "")) +
  scale_y_continuous(breaks = seq(-26, 22, by = 1)) +
  geom_point(x = 2009, y = d13c_thinn_dry, shape = 6, size = 3.5, color = vadelma1) +
  annotate("text", x = 2009, y = d13c_thinn_dry, label = "thinning", hjust = 0.5, vjust = -1, size = 2.8, color = vadelma1) +  
  geom_point(x = 2014, y = d13c_af_dry, shape = 1, size = 3.5, color = vadelma1) +
  annotate("text", x = 2014, y = d13c_af_dry, label = "AF", hjust = 0.5, vjust = -1, size = 2.8, color = vadelma1) +
  theme_bw(base_size = 10) +
  theme(legend.position = "top", legend.margin = margin(0,0,0,0, "cm"), legend.title = element_text(size = 10), legend.text = element_text(size = 10),
        axis.title.x = element_blank(), axis.title.y = element_blank(), axis.ticks.x = element_line(), axis.text.y = element_blank(), axis.ticks.length = unit(0.1, "cm"),
        panel.grid.major.y = element_line(color = kallio2, linewidth = 0.25, linetype = "dotted"),
        plot.margin = margin(0, 1, 2, 1, "mm"))

p_d13c_shal + p_d13c_deep + plot_layout(nrow = 1) 

Figure 7. δ¹³C lme model

Model for Muhos:

d13c_data_ansas <- d13c |> 
  filter(site == "Shallow WT") |> 
  mutate(
    yr_treat = factor(ifelse(treatment == "Treatment", pmax(year - 2014, -1), -1)),
    thinned = factor(treatment == "Treatment" & year > 2010),
    year = factor(year),
    treatment = relevel(treatment, ref = "Control")
  )

d13c_model_ansas <- lme(d13c ~ treatment + yr_treat + thinned, random = ~1|year, data = d13c_data_ansas)

as.data.frame(round(summary(d13c_model_ansas)$tTable, 3)) |> 
  mutate(p.adj.signif = sapply(`p-value`, function(p) {
           if (p < 0.001) { return("***") } 
           else if (p < 0.01) { return("**")} 
           else if (p < 0.05) { return("*") } 
           else if (p < 0.1) { return(".") } 
           else { return("ns") } })) 

                     Value Std.Error DF  t-value p-value p.adj.signif
(Intercept)        -24.613     0.141 17 -175.168   0.000          ***
treatmentTreatment  -0.057     0.060  7   -0.951   0.373           ns
yr_treat0           -0.210     0.169  7   -1.245   0.253           ns
yr_treat1            0.757     0.169  7    4.484   0.003           **
yr_treat2            0.276     0.169  7    1.635   0.146           ns
yr_treat3            0.290     0.169  7    1.719   0.129           ns
yr_treat4            0.322     0.169  7    1.909   0.098            .
yr_treat5            0.913     0.169  7    5.412   0.001           **
yr_treat6            0.403     0.169  7    2.389   0.048            *
yr_treat7            0.418     0.169  7    2.476   0.042            *
yr_treat8            0.721     0.169  7    4.276   0.004           **
thinnedTRUE          0.001     0.103  7    0.012   0.991           ns

d13c_model_ansas$sigma^2

[1] 0.01084408

as.numeric(getVarCov(d13c_model_ansas))

[1] 0.3445328

Sievi:

d13c_data_nevaj <- d13c |> 
  filter(site == "Deep WT") |> 
  mutate(
    yr_treat = factor(ifelse(treatment == "Treatment", pmax(year - 2014, -1), -1)),
    thinned = factor(treatment == "Treatment" & year > 2009),
    year = factor(year),
    treatment = relevel(treatment, ref = "Control")
  )

d13c_model_nevaj <- lme(d13c ~ treatment + yr_treat + thinned, random = ~1|year, data = d13c_data_nevaj)

as.data.frame(round(summary(d13c_model_nevaj)$tTable, 3)) |> 
  mutate(p.adj.signif = sapply(`p-value`, function(p) {
    if (p < 0.001) { return("***") } 
    else if (p < 0.01) { return("**")} 
    else if (p < 0.05) { return("*") }
    else if (p < 0.1) { return(".") } 
    else { return("ns") } })) 

                     Value Std.Error DF  t-value p-value p.adj.signif
(Intercept)        -24.342     0.139 17 -175.116   0.000          ***
treatmentTreatment   0.489     0.265  7    1.842   0.108           ns
yr_treat0            0.154     0.622  7    0.248   0.812           ns
yr_treat1           -0.126     0.622  7   -0.203   0.845           ns
yr_treat2            0.387     0.622  7    0.623   0.553           ns
yr_treat3           -0.308     0.622  7   -0.496   0.635           ns
yr_treat4            0.957     0.622  7    1.539   0.168           ns
yr_treat5            0.040     0.622  7    0.065   0.950           ns
yr_treat6           -0.603     0.622  7   -0.970   0.365           ns
yr_treat7           -0.156     0.622  7   -0.250   0.809           ns
yr_treat8            0.268     0.622  7    0.431   0.679           ns
thinnedTRUE         -0.246     0.373  7   -0.660   0.530           ns

d13c_model_nevaj$sigma^2

[1] 0.2317788

as.numeric(getVarCov(d13c_model_nevaj))

[1] 0.1160163

Plot:

p_d13c_lme_shal <- as.data.frame(summary(d13c_model_ansas)$tTable) |> 
  rownames_to_column("Term") |> 
  filter(grepl("^yr_treat", Term)) |> 
  mutate(Year = as.numeric(gsub("yr_treat", "", Term)),
         Lower = Value - 1.96 * Std.Error,
         Upper = Value + 1.96 * Std.Error,
         Site = "Shallow WT") |> 
  ggplot(aes(x = Year)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = kallio1) +
  geom_point(aes(y = Value, color = Site, shape = Site, size = Site)) +
  geom_line(aes(y = Value), color = vesi1) +
  geom_line(aes(y = Lower), color = vesi3) +
  geom_line(aes(y = Upper), color = vesi3) +
  labs(title = "a)", x = "Year after AF", y = expression(italic(delta)^13*C~treatment~effect~("\u2030")), color = "", shape = "", size = "") +
  scale_color_manual(values = c("Shallow WT" = vesi1)) +
  scale_shape_manual(values = c("Shallow WT" = 15)) +
  scale_size_manual(values = c("Shallow WT" = 2)) +
  scale_x_continuous(breaks = 0:10) +
  coord_cartesian(y = c(-1.6, 1.6)) +
  theme_bw(base_size = 10) +
  theme(legend.position = "top", legend.margin = margin(0,0,0,0, "cm"), legend.text  = element_text(size = 10),
        plot.margin = margin(0,0.1,0,0, "cm")) +
  annotate("text", x = 1, y = 0.75, label = "**", hjust = 0.5, vjust = 0, size = 6, color = vadelma1) +
  annotate("text", x = 4, y = 0.5, label = ".", hjust = 0.5, vjust = 0, size = 11, color = vadelma1) +
  annotate("text", x = 5, y = 0.9, label = "**", hjust = 0.5, vjust = 0, size = 6, color = vadelma1) +
  annotate("text", x = 6, y = 0.4, label = "*", hjust = 0.5, vjust = 0, size = 6, color = vadelma1) +
  annotate("text", x = 7, y = 0.42, label = "*", hjust = 0.5, vjust = 0, size = 6, color = vadelma1) +
  annotate("text", x = 8, y = 0.75, label = "**", hjust = 0.5, vjust = 0, size = 6, color = vadelma1)

p_d13c_lme_deep <- as.data.frame(summary(d13c_model_nevaj)$tTable) |> 
  rownames_to_column("Term") |> 
  filter(grepl("^yr_treat", Term)) |> 
  mutate(Year = as.numeric(gsub("yr_treat", "", Term)),
         Lower = Value - 1.96 * Std.Error,
         Upper = Value + 1.96 * Std.Error,
         Site = "Deep WT") |> 
  ggplot(aes(x = Year)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = kallio1) +
  geom_point(aes(y = Value, color = Site, shape = Site, size = Site)) +
  geom_line(aes(y = Value), color = valo1) +
  geom_line(aes(y = Lower), color = valo3) +
  geom_line(aes(y = Upper), color = valo3) +
  labs(title = "b)", x = "Year after AF", y = "", color = "", shape = "", size = "") +
  scale_color_manual(values = c("Deep WT" = valo1)) +
  scale_shape_manual(values = c("Deep WT" = 18)) +
  scale_size_manual(values = c("Deep WT" = 3)) +
  scale_x_continuous(breaks = 0:10) +
  coord_cartesian(y = c(-1.6, 1.6)) +
  theme_bw(base_size = 10) +
  theme(legend.position = "top", legend.margin = margin(0,0,0,0, "cm"), legend.text  = element_text(size = 10),
        plot.margin = margin(0,0,0,0, "cm"),
        axis.text.y = element_blank()) 

p_d13c_lme_shal + p_d13c_lme_deep

Figure 8. δ¹³C & WT correlations

Calculate the correlations:

# Muhos Treatment
WT_ansasF_data <- WT_summary |>
  right_join(d13c |> select(site, year, treatment, d13c), by = join_by(site == site, year == year, treatment == treatment)) |> 
  filter(site == "Shallow WT", treatment == "Treatment") |> 
  na.omit()

WT_ansasF_model <- lm(d13c ~ WT_cm, data = WT_ansasF_data)

WT_ansasF_intercept <- round(coef(WT_ansasF_model)[1], 2)
WT_ansasF_slope <- round(coef(WT_ansasF_model)[2], 2)
WT_ansasF_r2 <- round(summary(WT_ansasF_model)$r.squared, 3)

WT_ansasF_text <- paste0("Shallow WT, Treatment:\ny = ", WT_ansasF_intercept, ifelse(WT_ansasF_slope >= 0, " + ", " - "), abs(WT_ansasF_slope), "x\nR² = ", WT_ansasF_r2)

# Muhos Control
WT_ansasC_data <- WT_summary |>
  right_join(d13c |> select(site, year, treatment, d13c), by = join_by(site == site, year == year, treatment == treatment)) |> 
  filter(site == "Shallow WT", treatment == "Control") |> 
  na.omit()

WT_ansasC_model <- lm(d13c ~ WT_cm, data = WT_ansasC_data)

WT_ansasC_intercept <- round(coef(WT_ansasC_model)[1], 2)
WT_ansasC_slope <- round(coef(WT_ansasC_model)[2], 2)
WT_ansasC_r2 <- round(summary(WT_ansasC_model)$r.squared, 3)

WT_ansasC_text <- paste0("Shallow WT, Control:\ny = ", WT_ansasC_intercept, ifelse(WT_ansasC_slope >= 0, " + ", " - "), abs(WT_ansasC_slope), "x\nR² = ", WT_ansasC_r2)

# Sievi Treatment
WT_nevajF_data <- WT_summary |>
  right_join(d13c |> mutate(site = ifelse(site == "Deep WT", "Deep WT", "Shallow WT"), treatment = ifelse(treatment == "Control", "Control", "Treatment")) |> select(site, year, treatment, d13c), by = join_by(site == site, year == year, treatment == treatment)) |> 
  filter(site == "Deep WT" & treatment == "Treatment") |> 
  na.omit()

WT_nevajF_model <- lm(d13c ~ WT_cm, data = WT_nevajF_data)

WT_nevajF_intercept <- round(coef(WT_nevajF_model)[1], 2)
WT_nevajF_slope <- round(coef(WT_nevajF_model)[2], 2)
WT_nevajF_r2 <- round(summary(WT_nevajF_model)$r.squared, 3)

WT_nevajF_text <- paste0("Deep WT, Fertilizer:\ny = ", WT_nevajF_intercept, ifelse(WT_nevajF_slope >= 0, " + ", " - "), abs(WT_nevajF_slope), "x\nR² = ", WT_nevajF_r2)

# Sievi Control
WT_nevajC_data <- WT_summary |>
  right_join(d13c |> mutate(site = ifelse(site == "Deep WT", "Deep WT", "Shallow WT")) |> select(site, year, treatment, d13c), by = join_by(site == site, year == year, treatment == treatment)) |> 
  filter(site == "Deep WT" & treatment == "Control") |> 
  na.omit()

WT_nevajC_model <- lm(d13c ~ WT_cm, data = WT_nevajC_data)

WT_nevajC_intercept <- round(coef(WT_nevajC_model)[1], 2)
WT_nevajC_slope <- round(coef(WT_nevajC_model)[2], 3)
WT_nevajC_r2 <- round(summary(WT_nevajC_model)$r.squared, 3)

WT_nevajC_text <- paste0("Deep WT, Control:\ny = ", WT_nevajC_intercept, ifelse(WT_nevajC_slope >= 0, " + ", " - "), abs(WT_nevajC_slope), "x\nR² = ", WT_nevajC_r2)

Calculate r & p-values:

round(cor(WT_ansasF_data$WT_cm, WT_ansasF_data$d13c, use = "complete.obs"), 3)

[1] -0.668

round(cor.test(WT_ansasF_data$WT_cm, WT_ansasF_data$d13c)$p.value, 3)

[1] 0.006

round(cor(WT_ansasC_data$WT_cm, WT_ansasC_data$d13c, use = "complete.obs"), 3)

[1] -0.538

round(cor.test(WT_ansasC_data$WT_cm, WT_ansasC_data$d13c)$p.value, 3)

[1] 0.039

round(cor(WT_nevajF_data$WT_cm, WT_nevajF_data$d13c, use = "complete.obs"), 3)

[1] -0.368

round(cor.test(WT_nevajF_data$WT_cm, WT_nevajF_data$d13c)$p.value, 3)

[1] 0.16

round(cor(WT_nevajC_data$WT_cm, WT_nevajC_data$d13c, use = "complete.obs"), 3)

[1] -0.138

round(cor.test(WT_nevajC_data$WT_cm, WT_nevajC_data$d13c)$p.value, 3)

[1] 0.609

Draw plot:

p_WTd13C <- WT_summary |> 
  right_join(d13c |> mutate(treatment = ifelse(treatment == "Control", "Control", "Treatment")) |> select(site, year, treatment, d13c), by = join_by(site == site, year == year, treatment == treatment)) |> 
  na.omit() |> 
  mutate(site_treatment = factor(interaction(site, treatment), levels = c("Shallow WT.Control", "Shallow WT.Treatment", "Deep WT.Control", "Deep WT.Treatment"))) |> 
  ggplot(aes(x = WT_cm, y = d13c)) +
  geom_point(aes(group = site_treatment, color = site_treatment, shape = site_treatment, size = site_treatment)) +
  geom_smooth(aes(group = site_treatment, color = site_treatment, fill = site_treatment), method = "lm", se = TRUE) +
  geom_smooth(aes(group = 1), method = "loess", se = TRUE, color = vadelma1, fill = NA) +
  annotate("text", x = -82, y = -26.1, label = WT_nevajC_text, hjust = 0, vjust = 0, size = 3, color = valo2) +
  annotate("text", x = -69, y = -26.67, label = WT_nevajF_text, hjust = 0, vjust = 0, size = 3, color = valo1) +  
  annotate("text", x = -50, y = -26.1, label = WT_ansasC_text, hjust = 0, vjust = 0, size = 3, color = vesi2) +
  annotate("text", x = -37, y = -26.67, label = WT_ansasF_text, hjust = 0, vjust = 0, size = 3, color = vesi1) +
  scale_color_manual(values = c("Shallow WT.Treatment" = vesi1, "Shallow WT.Control" = vesi2, "Deep WT.Treatment" = valo1, "Deep WT.Control" = valo2, "Treatment" = vesi2, "Control" = vesi3),
                     labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_shape_manual(values = c("Shallow WT.Treatment" = 15, "Shallow WT.Control" = 0, "Deep WT.Treatment" = 18, "Deep WT.Control" = 5),
                     labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_size_manual(values = c("Shallow WT.Treatment" = 2, "Shallow WT.Control" = 2,  "Deep WT.Treatment" = 3, "Deep WT.Control" = 2),
                    labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_fill_manual(values = c("Shallow WT.Treatment" = vesi1, "Shallow WT.Control" = vesi2, "Deep WT.Treatment" = valo1, "Deep WT.Control" = valo2),
                    labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  labs(title = "", x = "Water Table (cm)", y = expression(italic(delta)^13*C~("\u2030")), color = "", shape = "", size = "", fill = "") +
  coord_cartesian(x = c(-83, -18), y = c(-26.6, -22.7)) +
  theme_bw(base_size = 10) +
  theme(legend.position = "right", legend.text  = element_text(size = 10))


p_WTd13C + plot_layout(ncol = 1) 

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

Figure 9. Needles

Read data:

needle <- read_excel("Neulasravinteet_Ansa_Neva.xlsx", sheet = "R all", na="") |> 
  filter(year != 2016) |> 
  mutate(site = factor(site, levels = c("Ansasaari", "Nevajärvi"), labels = c("Shallow WT", "Deep WT")),
         treatment = factor(treatment, levels = c("Control", "Fertilization"), labels = c("Control", "Treatment")),
         site_treatment = factor(interaction(site, treatment), levels = c("Shallow WT.Control", "Shallow WT.Treatment", "Deep WT.Control", "Deep WT.Treatment")),
         plot = factor(ifelse(year == 2014, "Before", "After"), levels = c("Before", "After")),
         mass_g100 = as.numeric(mass_g100),
         N = `N_%ka`,
         C = `C_%ka`,
         P_gkg = as.numeric(`P_mgkg ka` / 1000), 
         K_gkg = as.numeric(`K_mgkg ka` / 1000),
         B_mgkg = as.numeric(`B_mgkg ka`),
         length = as.numeric(length)) |> 
  select(site, treatment, year, mass_g100, length, N, C, P_gkg, K_gkg, B_mgkg, site_treatment, plot) 

Warning: There were 2 warnings in `mutate()`.
The first warning was:
ℹ In argument: `mass_g100 = as.numeric(mass_g100)`.
Caused by warning:
! NAs introduced by coercion
ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.

Needle length percentual differences between plots

needle |> 
  filter(length > 0) |> 
  group_by(site, treatment, year) |> 
  summarize(length = mean(length)) |> 
  pivot_wider(names_from = treatment, values_from = length) |> 
  mutate(difference = (Treatment - Control) / Control)

`summarise()` has grouped output by 'site', 'treatment'. You can override using
the `.groups` argument.

# A tibble: 3 × 5
# Groups:   site [2]
  site        year Control Treatment difference
  <fct>      <dbl>   <dbl>     <dbl>      <dbl>
1 Shallow WT  2019      44      54.5      0.239
2 Shallow WT  2020      46      53        0.152
3 Deep WT     2020      42      50        0.190

Mass of 100 needles, percentual differences between plots

needle |> 
  filter(mass_g100 > 0) |> 
  group_by(site, treatment, year) |> 
  summarize(mass = mean(mass_g100)) |> 
  pivot_wider(names_from = treatment, values_from = mass) |> 
  mutate(difference = (Treatment - Control) / Control)

`summarise()` has grouped output by 'site', 'treatment'. You can override using
the `.groups` argument.

# A tibble: 5 × 5
# Groups:   site [2]
  site        year Control Treatment difference
  <fct>      <dbl>   <dbl>     <dbl>      <dbl>
1 Shallow WT  2019    1.94      2.82     0.450 
2 Shallow WT  2020    2.04      3.10     0.513 
3 Shallow WT  2023    1.63      2.15     0.319 
4 Deep WT     2020    1.38      2        0.449 
5 Deep WT     2023    1.4       1.46     0.0429

Draw plot:

# Length
p_needleleng <- needle |>
  filter(length > 0 & year == 2020) |> 
  ggplot(aes(x = site, y = length, group = site_treatment, fill = site_treatment, color = site_treatment)) +
  stat_summary(fun = mean, geom = "bar", position = position_dodge(0.7), width = 0.5) +
  scale_size_continuous(range = c(0.7,1), guide = "none") +
  labs(title = "a)", x = "", y = "Needle length (mm)", fill = "", color="") +
  scale_fill_manual(values = c("Shallow WT.Control" = "white", "Shallow WT.Treatment" = vesi1, "Deep WT.Control" = "white", "Deep WT.Treatment" = valo1), 
                    labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_color_manual(values = c("Shallow WT.Control" = vesi1, "Shallow WT.Treatment" = vesi1, "Deep WT.Control" = valo1, "Deep WT.Treatment" = valo1), 
                     labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  facet_wrap(~ year, nrow = 1) +
  theme_bw(base_size = 10) +
  theme(legend.position = "none", 
        panel.grid.major.y = element_line(color = kallio2, linetype = "dotted"),
        axis.ticks.x = element_blank(), axis.text.x = element_blank(),
        strip.text = element_text(size = 9), strip.background = element_blank(), 
        plot.margin = margin(1,4,0,0, "mm"))


# Mass (g) of 100 needles
p_needlemass <- needle |>
  filter(mass_g100 > 0 & year != 2019) |> 
  ggplot(aes(x = site, y = mass_g100, group = site_treatment, fill = site_treatment, color = site_treatment)) +
  stat_summary(fun = mean, geom = "bar", position = position_dodge(0.7), width = 0.5) +
  scale_size_continuous(range = c(0.7,1), guide = "none") +
  labs(title = "b)", x = "", y = "The mass of 100 needles\n(g dw)", fill = "", color="") +
  scale_fill_manual(values = c("Shallow WT.Control" = "white", "Shallow WT.Treatment" = vesi1, "Deep WT.Control" = "white", "Deep WT.Treatment" = valo1), 
                    labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_color_manual(values = c("Shallow WT.Control" = vesi1, "Shallow WT.Treatment" = vesi1, "Deep WT.Control" = valo1, "Deep WT.Treatment" = valo1), 
                     labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  facet_wrap(~ year, nrow = 1) +
  theme_bw(base_size = 10) +
  theme(legend.text  = element_text(size = 10), 
        panel.grid.major.y = element_line(color = kallio2, linetype = "dotted"),
        axis.ticks.x = element_blank(), axis.text.x = element_blank(),
        strip.text = element_text(size = 9), strip.background = element_blank(),
        plot.margin = margin(1,4,0,0, "mm"))

p_needleab <- plot_grid(p_needleleng + p_needlemass + plot_layout(widths = c(1,2)))


# (Reinikainen et al. 1998)
# N 1.30 - 1.80 %
# P 1.6 - 2.2 g/kg
# K 4.5 - 6.5 g/kg
# B 10 - 30 mg/kg

p_needle_N <- needle |> 
  ggplot(aes(x = site, y = N)) +
  geom_hline(yintercept = c(1.3, 1.8), linetype = "dashed", color = vadelma1) +
  stat_summary(aes(group = site_treatment, fill = site_treatment, color = site_treatment), fun = mean, geom = "bar", position = position_dodge(0.7), width = 0.5) +
  scale_size_continuous(range=c(0.7,1), guide="none") +
  labs(title = "c)", x = "", y = "Nitrogen (%)", fill = "", color="") +
  scale_fill_manual(values = c("Shallow WT.Control" = "white", "Shallow WT.Treatment" = vesi1, "Deep WT.Control" = "white", "Deep WT.Treatment" = valo1), 
                    labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_color_manual(values = c("Shallow WT.Control" = vesi1, "Shallow WT.Treatment" = vesi1, "Deep WT.Control" = valo1, "Deep WT.Treatment" = valo1), 
                     labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.2))) +
  facet_wrap(~ year, nrow = 1) +
  theme_bw(base_size = 10) +
  theme(legend.position = "none", 
        panel.grid.major.y = element_line(color = kallio2, linetype = "dotted"),
        axis.text.x = element_blank(), axis.ticks.x = element_blank(),
        strip.text = element_text(size = 9), strip.background = element_blank(),
        plot.margin = margin(0,1,0,0, "mm")) 


p_needle_P <- needle |> 
  ggplot(aes(x = site, y = P_gkg)) +
  geom_hline(yintercept = c(1.6, 2.2), linetype = "dashed", color = vadelma1) +
  stat_summary(aes(group = site_treatment, fill = site_treatment, color = site_treatment), fun = mean, geom = "bar", position = position_dodge(0.7), width = 0.5) +
  scale_size_continuous(range=c(0.7,1), guide="none") +
  labs(title = "d)", x = "", y = "Phosphorus (g/kg)", fill = "", color="") +
  scale_fill_manual(values = c("Shallow WT.Control" = "white", "Shallow WT.Treatment" = vesi1, "Deep WT.Control" = "white", "Deep WT.Treatment" = valo1), 
                    labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_color_manual(values = c("Shallow WT.Control" = vesi1, "Shallow WT.Treatment" = vesi1, "Deep WT.Control" = valo1, "Deep WT.Treatment" = valo1), 
                     labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.2))) +
  facet_wrap(~ year, nrow = 1) +
  theme_bw(base_size = 10) +
  theme(legend.position = "none", 
        panel.grid.major.y = element_line(color = kallio2, linetype = "dotted"),
        axis.text.x = element_blank(), axis.ticks.x = element_blank(),
        strip.text = element_text(size = 9), strip.background = element_blank(),
        plot.margin = margin(0,1,0,0, "mm")) 


p_needle_K <- needle |> 
  ggplot(aes(x = site, y = K_gkg)) +
  geom_hline(yintercept = c(4.5, 6.5), linetype = "dashed", color = vadelma1) +
  stat_summary(aes(group = site_treatment, fill = site_treatment, color = site_treatment), fun = mean, geom = "bar", position = position_dodge(0.7), width = 0.5) +
  scale_size_continuous(range=c(0.7,1), guide="none") +
  labs(title = "e)", x = "", y = "Potassium (g/kg)", fill = "", color="") +
  scale_fill_manual(values = c("Shallow WT.Control" = "white", "Shallow WT.Treatment" = vesi1, "Deep WT.Control" = "white", "Deep WT.Treatment" = valo1), 
                    labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_color_manual(values = c("Shallow WT.Control" = vesi1, "Shallow WT.Treatment" = vesi1, "Deep WT.Control" = valo1, "Deep WT.Treatment" = valo1), 
                     labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.2))) +
  facet_wrap(~ year, nrow = 1) +
  theme_bw(base_size = 10) +
  theme(legend.position = "none", 
        panel.grid.major.y = element_line(color = kallio2, linetype = "dotted"),
        axis.text.x = element_blank(), axis.ticks.x = element_blank(),
        strip.text = element_text(size = 9), strip.background = element_blank(), 
        plot.margin = margin(0,1,0,0, "mm")) 


p_needle_B <- needle |> 
  ggplot(aes(x = site, y = B_mgkg)) +
  geom_hline(yintercept = c(10, 30), linetype = "dashed", color = vadelma1) +
  stat_summary(aes(group = site_treatment, fill = site_treatment, color = site_treatment), fun = mean, geom = "bar", position = position_dodge(0.7), width = 0.5) +
  scale_size_continuous(range=c(0.7,1), guide="none") +
  labs(title = "f)", x = "", y = "Boron (mg/kg)", fill = "", color="") +
  scale_fill_manual(values = c("Shallow WT.Control" = "white", "Shallow WT.Treatment" = vesi1, "Deep WT.Control" = "white", "Deep WT.Treatment" = valo1), 
                    labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_color_manual(values = c("Shallow WT.Control" = vesi1, "Shallow WT.Treatment" = vesi1, "Deep WT.Control" = valo1, "Deep WT.Treatment" = valo1), 
                     labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.2))) +
  facet_wrap(~ year, nrow = 1) +
  theme_bw(base_size = 10) +
  theme(legend.position = "none", 
        panel.grid.major.y = element_line(color = kallio2, linetype = "dotted"),
        axis.text.x = element_blank(), axis.ticks.x = element_blank(),
        strip.text = element_text(size = 9), strip.background = element_blank(),
        plot.margin = margin(0, 1, 0, 0, "mm")) 


p_needlecdef <- p_needle_N + p_needle_P + p_needle_K + p_needle_B + plot_layout(ncol = 1)

p_needles <- p_needleab + p_needlecdef + plot_layout(ncol = 1, heights = c(1, 3))

p_needles

Figure 1S. Peat properties

Read data:

peat <- read_excel("turpeen_ravinteisuus_suohitu.xls", sheet = "R", na="") |> 
  mutate(site = factor(site, levels = c("Ansasaari", "Nevajärvi"), labels = c("Shallow WT", "Deep WT")),
         treatment = factor(treatment, levels = c("Control", "Fertilized"), labels = c("Control", "Treatment")),
         site_treatment = factor(interaction(site, treatment), levels = c("Shallow WT.Control", "Shallow WT.Treatment", "Deep WT.Control", "Deep WT.Treatment")),
         plot = factor(plot),
         depth = factor(depth, levels = c("0-10", "10-20"), labels = c("0 - 10 cm", "10 - 20 cm")),
         C = as.numeric(`C_%ka`),
         N = as.numeric(`N_%ka`),
         B = as.numeric(B_mgkgka) / 1000000 * 100,
         K = as.numeric(K_mgkgka) / 1000000 * 100, 
         P = as.numeric(P_mgkgka) / 1000000 * 100) |> 
  filter(plot == 1) |> 
  select(site, sample, treatment, site_treatment, plot, depth, N, P, K, B)

Draw plot:

# plot N
peat_N <- peat |> 
  ggplot(aes(x = site, y = N, group = site_treatment, fill = site_treatment, color = site_treatment)) + 
  stat_summary(fun = mean, geom = "bar", position = position_dodge(0.7), width = 0.5) +
  stat_summary(fun.data = mean_se, geom = "errorbar", position = position_dodge(0.7), width = 0.2, color = "black") +
  labs(title = "a)", x = "", y = "Peat N content (% dw)", fill = "", color = "") +
  scale_fill_manual(values = c("Shallow WT.Control" = "white", "Shallow WT.Treatment" = vesi1, "Deep WT.Control" = "white", "Deep WT.Treatment" = valo1), 
                    labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_color_manual(values = c("Shallow WT.Control" = vesi1, "Shallow WT.Treatment" = vesi1, "Deep WT.Control" = valo1, "Deep WT.Treatment" = valo1), 
                     labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  facet_wrap(~ depth, nrow = 2) +
  theme_bw(base_size = 10) +
  theme(legend.position = "top", legend.text  = element_text(size = 10), 
        panel.grid.major.y = element_line(color = kallio2, linetype = "dotted"),
        axis.ticks.x = element_blank(), axis.text.x = element_blank(),  axis.title.y = element_text(size = 11),
        strip.text = element_text(size = 11), strip.background = element_blank(),
        plot.margin = margin(1,4,0,0, "mm")) +
  guides(col = guide_legend(nrow = 2, ncol = 2))


# plot P
peat_P <- peat |> 
  ggplot(aes(x = site, y = P, group = site_treatment, fill = site_treatment, color = site_treatment)) + 
  stat_summary(fun = mean, geom = "bar", position = position_dodge(0.7), width = 0.5) +
  stat_summary(fun.data = mean_se, geom = "errorbar", position = position_dodge(0.7), width = 0.2, color = "black") +
  labs(title = "b)", x = "", y = expression("Peat P content (mg kg"^-1*" dw)"), fill = "") +
  scale_fill_manual(values = c("Shallow WT.Control" = "white", "Shallow WT.Treatment" = vesi1, "Deep WT.Control" = "white", "Deep WT.Treatment" = valo1), 
                    labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_color_manual(values = c("Shallow WT.Control" = vesi1, "Shallow WT.Treatment" = vesi1, "Deep WT.Control" = valo1, "Deep WT.Treatment" = valo1), 
                     labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  facet_wrap(~ depth, nrow = 2) +
  theme_bw(base_size = 10) +
  theme(legend.position = "none",
        panel.grid.major.y = element_line(color = kallio2, linetype = "dotted"),
        axis.ticks.x = element_blank(), axis.text.x = element_blank(),  axis.title.y = element_text(size = 11),
        strip.text = element_text(size = 11), strip.background = element_blank(),
        plot.margin = margin(1,4,0,0, "mm"))

# plot K
peat_K <- peat |> 
  ggplot(aes(x = site, y = K, group = site_treatment, fill = site_treatment, color = site_treatment)) + 
  stat_summary(fun = mean, geom = "bar", position = position_dodge(0.7), width = 0.5) +
  stat_summary(fun.data = mean_se, geom = "errorbar", position = position_dodge(0.7), width = 0.2, color = "black") +
  labs(title = "c)", x = "", y = expression("Peat K content (mg kg"^-1*" dw)"), fill = "") +
  scale_fill_manual(values = c("Shallow WT.Control" = "white", "Shallow WT.Treatment" = vesi1, "Deep WT.Control" = "white", "Deep WT.Treatment" = valo1), 
                    labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_color_manual(values = c("Shallow WT.Control" = vesi1, "Shallow WT.Treatment" = vesi1, "Deep WT.Control" = valo1, "Deep WT.Treatment" = valo1), 
                     labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  facet_wrap(~ depth, nrow = 2) +
  theme_bw(base_size = 10) +
  theme(legend.position = "none",
        panel.grid.major.y = element_line(color = kallio2, linetype = "dotted"),
        axis.ticks.x = element_blank(), axis.text.x = element_blank(),  axis.title.y = element_text(size = 11),
        strip.text = element_text(size = 11), strip.background = element_blank(),
        plot.margin = margin(1,4,0,0, "mm"))


# plot B
peat_B <- peat |> 
  ggplot(aes(x = site, y = B, group = site_treatment, fill = site_treatment, color = site_treatment)) + 
  stat_summary(fun = mean, geom = "bar", position = position_dodge(0.7), width = 0.5) +
  stat_summary(fun.data = mean_se, geom = "errorbar", position = position_dodge(0.7), width = 0.2, color = "black") +
  labs(title = "d)", x = "", y = expression("Peat B content (mg kg"^-1*" dw)"), fill = "") +
  scale_fill_manual(values = c("Shallow WT.Control" = "white", "Shallow WT.Treatment" = vesi1, "Deep WT.Control" = "white", "Deep WT.Treatment" = valo1), 
                    labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_color_manual(values = c("Shallow WT.Control" = vesi1, "Shallow WT.Treatment" = vesi1, "Deep WT.Control" = valo1, "Deep WT.Treatment" = valo1), 
                     labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  facet_wrap(~ depth, nrow = 2) +
  theme_bw(base_size = 10) +
  theme(legend.position = "none",
        panel.grid.major.y = element_line(color = kallio2, linetype = "dotted"),
        axis.ticks.x = element_blank(), axis.text.x = element_blank(), axis.title.y = element_text(size = 11),
        strip.text = element_text(size = 11), strip.background = element_blank(),
        plot.margin = margin(1,4,0,0, "mm"))

legend_top <- get_legend(peat_N)

peat_N <- peat_N + theme(legend.position = "none")

plot_grid(
  plot_grid(legend_top), 
  peat_N + peat_P + peat_K + peat_B + plot_layout(nrow = 1),
  ncol = 1, 
  rel_heights = c(0.1, 1)
) 

Figure S2. Needle chemistry - molar ratios

Draw plot:

p_CN <- needle |> 
  group_by(site, treatment, year) |> 
  mutate(N_gkg = N * 10,
         C_gkg = C * 10,
         C_mol = C_gkg / 12.011,
         N_mol = N_gkg / 14.007,
         CN = C_mol / N_mol) |> 
  ggplot(aes(x = site, y = CN)) +
  stat_summary(aes(group = site_treatment, fill = site_treatment, color = site_treatment), fun = mean, geom = "bar", position = position_dodge(0.7), width = 0.5) +
  scale_size_continuous(range=c(0.7,1), guide="none") +
  labs(title = "a)", x = "", y = "C:N", fill = "", color="") +
  scale_fill_manual(values = c("Shallow WT.Control" = "white", "Shallow WT.Treatment" = vesi1, "Deep WT.Control" = "white", "Deep WT.Treatment" = valo1), 
                    labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_color_manual(values = c("Shallow WT.Control" = vesi1, "Shallow WT.Treatment" = vesi1, "Deep WT.Control" = valo1, "Deep WT.Treatment" = valo1), 
                     labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.2))) +
  facet_wrap(~ year, nrow = 1) +
  theme_bw(base_size = 10) +
  theme(legend.position = "none", 
        panel.grid.major.y = element_line(color = kallio2, linetype = "dotted"),
        axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.title.y = element_text(size = 11),
        strip.text = element_text(size = 11), strip.background = element_blank(),
        plot.margin = margin(0, 1, 0, 0, "mm")) 


p_NP <- needle |> 
  group_by(site, treatment, year) |> 
  mutate(N_gkg = N * 10,
         N_mol = N_gkg / 14.007,
         P_mol = P_gkg / 30.974,
         NP = N_mol / P_mol) |> 
  ggplot(aes(x = site, y = NP)) +
  stat_summary(aes(group = site_treatment, fill = site_treatment, color = site_treatment), fun = mean, geom = "bar", position = position_dodge(0.7), width = 0.5) +
  scale_size_continuous(range=c(0.7,1), guide="none") +
  labs(title = "b)", x = "", y = "N:P", fill = "", color="") +
  scale_fill_manual(values = c("Shallow WT.Control" = "white", "Shallow WT.Treatment" = vesi1, "Deep WT.Control" = "white", "Deep WT.Treatment" = valo1), 
                    labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_color_manual(values = c("Shallow WT.Control" = vesi1, "Shallow WT.Treatment" = vesi1, "Deep WT.Control" = valo1, "Deep WT.Treatment" = valo1), 
                     labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.2))) +
  facet_wrap(~ year, nrow = 1) +
  theme_bw(base_size = 10) +
  theme(legend.position = "none", 
        panel.grid.major.y = element_line(color = kallio2, linetype = "dotted"),
        axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.title.y = element_text(size = 11),
        strip.text = element_text(size = 11), strip.background = element_blank(),
        plot.margin = margin(0, 1, 0, 0, "mm")) 

p_NK <- needle |> 
  group_by(site, treatment, year) |> 
  mutate(N_gkg = N * 10,
         N_mol = N_gkg / 14.007,
         K_mol = K_gkg / 39.0983,
         NK = N_mol / K_mol) |> 
  ggplot(aes(x = site, y = NK)) +
  stat_summary(aes(group = site_treatment, fill = site_treatment, color = site_treatment), fun = mean, geom = "bar", position = position_dodge(0.7), width = 0.5) +
  scale_size_continuous(range=c(0.7,1), guide="none") +
  labs(title = "c)", x = "", y = "N:K", fill = "", color="") +
  scale_fill_manual(values = c("Shallow WT.Control" = "white", "Shallow WT.Treatment" = vesi1, "Deep WT.Control" = "white", "Deep WT.Treatment" = valo1), 
                    labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_color_manual(values = c("Shallow WT.Control" = vesi1, "Shallow WT.Treatment" = vesi1, "Deep WT.Control" = valo1, "Deep WT.Treatment" = valo1), 
                     labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.2))) +
  facet_wrap(~ year, nrow = 1) +
  theme_bw(base_size = 10) +
  theme(legend.position = "none", 
        panel.grid.major.y = element_line(color = kallio2, linetype = "dotted"),
        axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.title.y = element_text(size = 11),
        strip.text = element_text(size = 11), strip.background = element_blank(),
        plot.margin = margin(0, 1, 0, 0, "mm")) 

p_KP <- needle |> 
  group_by(site, treatment, year) |> 
  mutate(K_mol = K_gkg / 39.0983,
         P_mol = P_gkg / 30.974,
         KP = K_mol / P_mol) |> 
  ggplot(aes(x = site, y = KP)) +
  stat_summary(aes(group = site_treatment, fill = site_treatment, color = site_treatment), fun = mean, geom = "bar", position = position_dodge(0.7), width = 0.5) +
  scale_size_continuous(range=c(0.7,1), guide="none") +
  labs(title = "d)", x = "", y = "K:P", fill = "", color="") +
  scale_fill_manual(values = c("Shallow WT.Control" = "white", "Shallow WT.Treatment" = vesi1, "Deep WT.Control" = "white", "Deep WT.Treatment" = valo1), 
                    labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_color_manual(values = c("Shallow WT.Control" = vesi1, "Shallow WT.Treatment" = vesi1, "Deep WT.Control" = valo1, "Deep WT.Treatment" = valo1), 
                     labels = c("Control with shallow WT", "Treatment with shallow WT", "Control with deep WT", "Treatment with deep WT")) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.2))) +
  facet_wrap(~ year, nrow = 1) +
  theme_bw(base_size = 10) +
  theme(legend.position = "none", 
        panel.grid.major.y = element_line(color = kallio2, linetype = "dotted"),
        axis.text.x = element_blank(), axis.ticks.x = element_blank(), axis.title.y = element_text(size = 11),
        strip.text = element_text(size = 11), strip.background = element_blank(),
        plot.margin = margin(0, 1, 0, 0, "mm"))


plot_grid(legend_top) + p_CN + p_NP + p_NK + p_KP + plot_layout(ncol = 1, heights = c(0.4, 1, 1, 1, 1))