---
title: "Marchand et al. 2024 Nature Ecology & Evolution"
author: "Lorene J. Marchand"
date: "2024-10-30"
output:
  html_document:
    toc: true
    toc_float: true
    toc_depth: 4
---

This script comprise the analysis of the paper entitled : __No winter halt in belowground wood growth of four angisperm deciduous tree species__

Lorène J. Marchand, Jožica Gričar, Paolo Zuccarini, Inge Dox, Bertold Mariën, Melanie Verlinden, Thilo Heinecke, Peter Prislan, Guillaume Marie, Holger Lange, Jan Van den Bulcke, Josep Penuelas, Patrick Fonti, Matteo Campioli

The analysis comprise (i) wood growth dataset, (ii) leaf phenology dataset (iii) climate dataset and (iv) fine root dataset

(i) __Wood growth dataset__  contain all the measurements of xylogenesis of stem and roots for 
* _Mature beech in 2018 (Belgium)_
* _Mature beech and birch in 2020 (Belgium)_
* _Sapling beech, birch, oak and poplar in 2017 (Belgium)_
* _Sapling beech and birch in 2018 (Belgium)_
* _Sapling beech and birch in 2018 (Spain)_
* _Sapling beech and birch in 2018 (Norway)_

The length of the layer of cell in enlargement is defined as __"EC"__ , the length of the layer of cell in secondary cell wall of thickening is defined as __"TH"__,  the total ring length is defined as __"ring"__

DOY_cumul reefers to cumulative DOY starting from the first year of measurement and continuing the following year. 

(ii) __Phenology dataset__ contains the spring leaf, autumn leaf and fine root phenology observation  for 
* _Mature beech in 2018 (Belgium)_
* _Mature beech and birch in 2020 (Belgium)_
* _Sapling beech, birch, oak and poplar in 2017 (Belgium)_
* _Sapling beech and birch in 2018 (Belgium)_
* _Sapling beech and birch in 2018 (Spain)_
* _Sapling beech and birch in 2018 (Norway)_

file __senescence_dataset.csv__ The autumn phenology variable __"Senescence"__ is calculated with the formulate in _Vitasse, Y., Porté, A. J., Kremer, A., Michalet, R. & Delzon, S. Responses of canopy duration to temperature changes in four temperate tree species: relative contributions of spring and autumn leaf phenology. Oecologia 161, 187–198 (2009)_.


file __Bud_leaf_dev.csv__ The spring phenology variable __"Bud_Leaf_dev"__ is calculated with the methods in _Marchand, L. J. et al. Inter-individual variability in spring phenology of temperate deciduous trees depends on species, tree size and previous year autumn phenology. Agric. For. Meteorol. 290, 108031 (2020)_ and in _Zuccarini, P. et al. Drivers and dynamics of foliar senescence in temperate deciduous forest trees at their southern limit of distribution in Europe. Agric. For. Meteorol. 342, 109716 (2023)._


file __fine_root.csv__ the fine root phenology from the file __fine root dataset__ cumulative fine root surface is calculated form the formula in _Mariën, B. et al. On the Below- and Aboveground Phenology in Deciduous Trees: Observing the Fine-Root Lifespan, Turnover Rate, and Phenology of  Fagus sylvatica L.,  Quercus robur L., and  Betula pendula Roth for Two Growing Seasons. Forests 12, 1680 (2021)._


(iii) __Climate dataset__ contains in the file __climate_soil_2017.csv__ the estimation of the soil temperature based on the sensor of another potted experiment with similar soil composition at the same site. Publiched in the following paper: https://onlinelibrary.wiley.com/doi/abs/10.1111/gcb.15501. the file __climate_2020.csv__ and __climate_2018.csv__ contain soil and air temperature of the site of Brasschaat 


__Packages used__ 
```{r echo=TRUE, message=FALSE, warning=FALSE, results='hide'}
library(dplyr)
library(tidyr)
library(kableExtra)
library(mgcv)
library(ggeffects)
library(ggplot2)
library(patchwork)
library(cowplot)
library(rpart)
library(rpart.plot)
```

# Xylogenesis dataset
```{r echo=TRUE, message=FALSE, warning=FALSE, results='hide'}
xylogenesis<-read.table("xylogenesis.csv", header=T, sep=";", dec=",", na.strings = "NA")
xylogenesis$country<-as.factor(xylogenesis$country)
xylogenesis$species<-as.factor(xylogenesis$species)
xylogenesis$site<-as.factor(xylogenesis$site)
xylogenesis$ferti<-as.factor(xylogenesis$ferti)
xylogenesis$plot<-as.factor(xylogenesis$plot)
xylogenesis$tree<-as.factor(xylogenesis$tree)
xylogenesis$organ<-as.factor(xylogenesis$organ)
xylogenesis$number<-as.factor(xylogenesis$number)
xylogenesis$date<-as.POSIXct(xylogenesis$date)
xylogenesis$DOY<-as.integer(xylogenesis$DOY)
xylogenesis$DOY_cumul<-as.integer(xylogenesis$DOY_cumul)
xylogenesis$exp<-as.factor(xylogenesis$exp)
xylogenesis$age<-as.factor(xylogenesis$age)
xylogenesis$lenght_EC<-as.numeric(xylogenesis$lenght_EC)
xylogenesis$lenght_TH<-as.numeric(xylogenesis$lenght_TH)
xylogenesis$lenght_ring<-as.numeric(xylogenesis$lenght_ring)
xylogenesis$code<-as.factor(xylogenesis$code)
```

## Roots 
we define the mean value per samples (4 measures) and create separate dataframe for each experiment
```{r echo=TRUE, message=FALSE, warning=FALSE, results='hide'}
roots<-subset(xylogenesis, organ=="roots")

roots2 <- roots %>% 
  group_by(code) %>%
  dplyr::summarise(sample_size = n(), percent_TH = mean(((lenght_TH)/lenght_ring)*100), percent_EC = mean(((lenght_EC)/lenght_ring))*100)

roots2<-separate(roots2, "code" ,  sep ="_", into=c("country","species","plot","tree","type","number","DOY_cumul"))

roots2$country<-as.factor(roots2$country)
roots2$species<-as.factor(roots2$species)
roots2$type<-as.factor(roots2$type)
roots2$DOY_cumul<-as.numeric(roots2$DOY_cumul)
roots2<-na.omit(roots2)
```

```{r echo=FALSE, message=FALSE, warning=FALSE}
mature<-subset(roots2, type=="mature")
root_mature_2020<-subset(mature, species== "Betula" | (species=="Fagus" & plot=="Park"))
head(root_mature_2020, 3)%>%
  kbl(caption = "Root mature trees 2020") %>%
  kable_classic(full_width = F, html_font = "Calibri")
root_mature_2018<-subset(mature, (species=="Fagus" & plot=="Military"))
head(root_mature_2018, 3)%>%
  kbl(caption = "Root mature trees 2020") %>%
  kable_classic(full_width = F, html_font = "Calibri")

root_sapling_2017<-subset(roots2, type=="sapling" & (number=="R"))
head(root_sapling_2017, 3)%>%
  kbl(caption = "Root sapling 2017") %>%
  kable_classic(full_width = F, html_font = "Calibri")

exp_2018<-subset(roots2, type=="sapling" & (number=="R1" | number=="R2"| number=="R3" ))
root_sapling_2018_belgium<-subset(exp_2018, country=="belgium")
head(root_sapling_2018_belgium, 3)%>%
  kbl(caption = "Root sapling 2018 Belgium") %>%
  kable_classic(full_width = F, html_font = "Calibri")
root_sapling_2018_spain<-subset(exp_2018, country=="spain")
head(root_sapling_2018_spain, 3)%>%
  kbl(caption = "Root sapling 2018 Spain") %>%
  kable_classic(full_width = F, html_font = "Calibri")
root_sapling_2018_norway<-subset(exp_2018, country=="norway")
head(root_sapling_2018_norway, 3)%>%
  kbl(caption = "Root sapling 2018 Norway") %>%
  kable_classic(full_width = F, html_font = "Calibri")

```

## Stems
we define the mean value per samples (3 measures) and create separate dataframe for each experiment
```{r echo=TRUE, message=FALSE, warning=FALSE }
stem<-subset(xylogenesis, organ=="stem")

stem2 <- stem %>% 
  group_by(code) %>%
  dplyr::summarise(sample_size = n(),  percent_TH = mean(((lenght_TH)/lenght_ring)*100), percent_EC = mean(((lenght_EC)/lenght_ring))*100)

stem2<-separate(stem2, "code" ,  sep ="_", into=c("country","species","plot","tree","type","number","DOY_cumul"))

stem2$country<-as.factor(stem2$country)
stem2$species<-as.factor(stem2$species)
stem2$type<-as.factor(stem2$type)
stem2$DOY_cumul<-as.numeric(stem2$DOY_cumul)
stem2$number<-as.factor(stem2$number)
stem2<-na.omit(stem2)
```


```{r echo=FALSE, message=FALSE, warning=FALSE}
mature<-subset(stem2, type=="mature")
stem_mature_2020_21<-subset(mature, species== "Betula" | (species=="Fagus" & plot=="Park"))
stem_mature_2021<-subset(stem_mature_2020_21, DOY_cumul>365)
head(stem_mature_2021, 3)%>%
  kbl(caption = "Stem mature trees 2021 spring") %>%
  kable_classic(full_width = F, html_font = "Calibri")
stem_mature_2020<-subset(stem_mature_2020_21, DOY_cumul<365)
head(stem_mature_2020, 3)%>%
  kbl(caption = "Stem mature trees 2020") %>%
  kable_classic(full_width = F, html_font = "Calibri")

stem_mature2018_19<-subset(mature,(species=="Fagus" & plot=="Military"))
stem_mature_2019<-subset(stem_mature2018_19, DOY_cumul>365)
head(stem_mature_2019, 3)%>%
  kbl(caption = "Stem mature trees 2019 spring") %>%
  kable_classic(full_width = F, html_font = "Calibri")
stem_mature_2018<-subset(stem_mature2018_19, DOY_cumul<365)
head(stem_mature_2018, 3)%>%
  kbl(caption = "Stem mature trees 2018") %>%
  kable_classic(full_width = F, html_font = "Calibri")

sapling<-subset(stem2, type=="sapling")

stem_sapling_2017<-subset(sapling, (number=="S"))
head(stem_sapling_2017, 3)%>%
  kbl(caption = "Stem sapling 2017") %>%
  kable_classic(full_width = F, html_font = "Calibri")
stem_sapling_2018_belgium<-subset(sapling, number=="S1" & country=="belgium")
head(stem_sapling_2018_belgium, 3)%>%
  kbl(caption = "Stem sapling 2018 Belgium") %>%
  kable_classic(full_width = F, html_font = "Calibri")
stem_sapling_2018_spain<-subset(sapling, number=="S1" & country=="spain")
head(stem_sapling_2018_spain, 3)%>%
  kbl(caption = "Stem sapling 2018 Spain") %>%
  kable_classic(full_width = F, html_font = "Calibri")
stem_sapling_2018_norway<-subset(sapling, number=="S1" & country=="norway")
head(stem_sapling_2018_norway, 3)%>%
  kbl(caption = "Stem sapling 2018 Norway") %>%
  kable_classic(full_width = F, html_font = "Calibri")
```

## Temperature vs thickening


```{r echo=TRUE, message=FALSE, warning=FALSE}
Temp_TH<-read.table("Temp_TH.csv", header=T, sep=";", dec=",", na.strings = "NA")
```

We create separate dataframe for each experiment
```{r echo=FALSE, message=FALSE, warning=FALSE}
temp_TH_mature_2020<- subset(Temp_TH, exp=="exp_mat_2020_")
head(temp_TH_mature_2020, 3)%>%
  kbl(caption = "Temperature vs. ring thickening mature trees 2020") %>%
  kable_classic(full_width = F, html_font = "Calibri")

temp_TH_mature_2018<- subset(Temp_TH, exp=="exp_mat_2018_")
head(temp_TH_mature_2020, 3)%>%
  kbl(caption = "Temperature vs. ring thickening  mature trees 2018") %>%
  kable_classic(full_width = F, html_font = "Calibri")

temp_TH_sampling_2017<- subset(Temp_TH, exp=="exp_young_bel_2017")
head(temp_TH_sampling_2017, 3)%>%
  kbl(caption = "Temperature vs. ring thickening  sapling 2017") %>%
  kable_classic(full_width = F, html_font = "Calibri")

temp_TH_sampling_2018_bel<- subset(Temp_TH, exp=="exp_young_bel_2018")
head(temp_TH_sampling_2017, 3)%>%
  kbl(caption = "Temperature vs. ring thickening  sapling 2018 belgium") %>%
  kable_classic(full_width = F, html_font = "Calibri")
```


# Phenology

## Senescence

```{r echo=TRUE, message=FALSE, warning=FALSE}
senescence<-read.table("senescence_dataset.csv", header=T, sep=";", dec=",", na.strings = "NA")
```

We create separate dataframe for each experiment
```{r echo=FALSE, message=FALSE, warning=FALSE}
mature<-subset(senescence, type=="mature")
senescence_mature_2020<-subset(mature, year=="2020")
head(senescence_mature_2020, 3)%>%
  kbl(caption = "Senescence mature trees 2020") %>%
  kable_classic(full_width = F, html_font = "Calibri")

senescence_mature_2018<-subset(mature, year=="2018")
head(senescence_mature_2018, 3)%>%
  kbl(caption = "Senescence mature trees 2020") %>%
  kable_classic(full_width = F, html_font = "Calibri")

sapling<-subset(senescence, type=="sapling")
senescence_sapling_2017<-subset(sapling,year=="2017")
head(senescence_sapling_2017, 3)%>%
  kbl(caption = "Senescence sapling 2017") %>%
  kable_classic(full_width = F, html_font = "Calibri")

senescence_sapling_2018_belgium<-subset(sapling, country=="belgium" & year=="2018")
head(senescence_sapling_2018_belgium, 3)%>%
  kbl(caption = "Senescence sapling 2018 Belgium") %>%
  kable_classic(full_width = F, html_font = "Calibri")
senescence_sapling_2018_spain<-subset(sapling, country=="spain" & year=="2018")
head(senescence_sapling_2018_spain, 3)%>%
  kbl(caption = "Senescence sapling 2018 Spain") %>%
  kable_classic(full_width = F, html_font = "Calibri")
senescence_sapling_2018_norway<-subset(sapling, country=="norway" & year=="2018")
head(senescence_sapling_2018_norway, 3)%>%
  kbl(caption = "Senescence sapling 2018 Norway") %>%
  kable_classic(full_width = F, html_font = "Calibri")

```

## bud and beaf development 

```{r echo=TRUE, message=FALSE, warning=FALSE}
bud<- read.table("Bud_Leaf_dev.csv", header=T, sep=";", dec=",", na.strings = "NA")
```

We create separate dataframe for each experiment
```{r echo=FALSE, message=FALSE, warning=FALSE}
mature<-subset(bud, type=="mature")
bud_mature_2020<-subset(mature, year=="2021")
head(bud_mature_2020, 3)%>%
  kbl(caption = "bud and leaf development mature trees 2021") %>%
  kable_classic(full_width = F, html_font = "Calibri")

bud_mature_2018<-subset(mature, year=="2019")
head(bud_mature_2018, 3)%>%
  kbl(caption = "bud and leaf development mature trees 2019") %>%
  kable_classic(full_width = F, html_font = "Calibri")

sapling<-subset(bud, type=="sapling")
bud_sapling_2017<-subset(bud,year=="2018")
head(bud_sapling_2017, 3)%>%
  kbl(caption = "bud and leaf development sapling 2018") %>%
  kable_classic(full_width = F, html_font = "Calibri")
```

## fine roots 

the cumulative root surface is calculated per day of measurements, for each site, species, tree, and tube. 
```{r echo=TRUE, message=FALSE, warning=FALSE}
fine_root<- read.table("fine_root.csv", header=T, sep=";", dec=",", na.strings = "NA")

fine_root_2 <- na.omit(fine_root)
fine_root_3 <-  fine_root_2 %>%
  group_by(Doy_tot, Site, Species, Tree_ID, Tube) %>%
  dplyr::summarise(cum.Root_surf = sum(TotSurfArea))
```


```{r echo=FALSE, message=FALSE, warning=FALSE}
fine_root_4<-subset(fine_root_3, Doy_tot > 590 & Doy_tot < 790)
new_Doy<-c(rep(230, time=18), rep(265, time =20), rep(296, time=19), rep(329, time = 11), rep(352, time = 18), rep (422, time=19))
fine_root_4$new_DOY<-new_Doy

fine_root_4$Species<-as.factor(fine_root_4$Species)

head(fine_root_4, 3)%>%
  kbl(caption = "fine roots dataset 2020") %>%
  kable_classic(full_width = F, html_font = "Calibri")
```

outup of the GAM model: 
```{r echo=FALSE, message=FALSE, warning=FALSE}
m_fine_root<-gam(cum.Root_surf ~ s(new_DOY,by = Species, k=5)+Species,  data=fine_root_4, family=nb, method="REML")


summary(m_fine_root)
par(mfrow = c(2, 2))
gam.check(m_fine_root)
est_fine_root<-ggpredict(m_fine_root, terms=c("new_DOY", "Species"))
```


# Climate 
```{r echo=FALSE, message=FALSE, warning=FALSE}
climate<-read.table("climate.csv", header=T, sep=";", dec=",", na.strings = "NA")
```

We create separate dataframe for each experiment
```{r echo=FALSE, message=FALSE, warning=FALSE}
mature<-subset(climate, type=="mature")
climate_mature_2020<-subset(mature, year=="2020")
head(climate_mature_2020, 3)%>%
  kbl(caption = "climate mature trees 2020") %>%
  kable_classic(full_width = F, html_font = "Calibri")

climate_mature_2018<-subset(mature, year=="2018")
head(climate_mature_2018, 3)%>%
  kbl(caption = "Climate mature trees 2020") %>%
  kable_classic(full_width = F, html_font = "Calibri")

sapling<-subset(climate, type=="sapling")
climate_sapling_2017<-subset(sapling,year=="2017")
head(climate_sapling_2017, 3)%>%
  kbl(caption = "Climate sapling 2017") %>%
  kable_classic(full_width = F, html_font = "Calibri")

climate_sapling_2018_belgium<-subset(sapling, country=="belgium" & year=="2018")
head(climate_sapling_2018_belgium, 3)%>%
  kbl(caption = "Climate sapling 2018 Belgium") %>%
  kable_classic(full_width = F, html_font = "Calibri")
climate_sapling_2018_spain<-subset(sapling, country=="spain" & year=="2018")
head(climate_sapling_2018_spain, 3)%>%
  kbl(caption = "Climate sapling 2018 Spain") %>%
  kable_classic(full_width = F, html_font = "Calibri")
climate_sapling_2018_norway<-subset(sapling, country=="norway" & year=="2018")
head(climate_sapling_2018_norway, 3)%>%
  kbl(caption = "Climate sapling 2018 Norway") %>%
  kable_classic(full_width = F, html_font = "Calibri")
```

# Analyse Fig. 3

## A

output of the GAM model

__Stem thickening saplings 2017__
```{r echo=FALSE, message=FALSE, warning=FALSE}
m_stem_sapling_2017<-gam(percent_TH/100 ~ s(DOY_cumul, by= species, k=3)+species,  data=stem_sapling_2017, family="quasibinomial", method="REML")
summary(m_stem_sapling_2017)
par(mfrow = c(2, 2))
gam.check(m_stem_sapling_2017)
est_stem_sapling_2017<-ggpredict(m_stem_sapling_2017, terms=c("DOY_cumul", "species"))
```

__Senescence saplings 2017__
```{r echo=FALSE, message=FALSE, warning=FALSE}
senescence_sapling_2017$species<-as.factor(senescence_sapling_2017$species)
m_senescence_sapling_2017<-gam(senescence/100 ~ s(DOY,by = species, k=3),  data=senescence_sapling_2017, family="quasibinomial", method="REML")
summary(m_senescence_sapling_2017)
est_senescence_sapling_2017<-ggpredict(m_senescence_sapling_2017, terms=c("DOY", "species"))
```

__Bud development saplings 2017__
```{r echo=FALSE, message=FALSE, warning=FALSE}
bud_sapling_2017$species<-as.factor(bud_sapling_2017$species)
m_bud_sapling_2017<-gam(score ~ s(DOY_cumul,by = species, k=3)+species,  data=bud_sapling_2017, family="quasibinomial", method="REML")
summary(m_bud_sapling_2017)
est_bud_sapling_2017<-ggpredict(m_bud_sapling_2017, terms=c("DOY_cumul", "species"))
```

__Senescence mature trees 2020__
```{r echo=FALSE, message=FALSE, warning=FALSE}
senescence_mature_2020$species<-as.factor(senescence_mature_2020$species)
m_senescence_mature_2020<-gam(senescence/100 ~ s(DOY,by = species, k=3),  data=senescence_mature_2020, family="quasibinomial", method="REML")
summary(m_senescence_mature_2020)
est_senescence_mature_2020<-ggpredict(m_senescence_mature_2020, terms=c("DOY", "species"))
```

__Bud development mature trees 2020__
```{r echo=FALSE, message=FALSE, warning=FALSE}
bud_mature_2020$species<-as.factor(bud_mature_2020$species)
m_bud_mature_2020<-gam(score/4 ~ s(DOY_cumul,by = species, k=3)+species,  data=bud_mature_2020, family="quasibinomial", method="REML")
summary(m_bud_mature_2020)
est_bud_mature_2020<-ggpredict(m_bud_mature_2020, terms=c("DOY_cumul","species"))
```

__Stem thickening mature trees 2020__
```{r echo=FALSE, message=FALSE, warning=FALSE}
stem_mature_2020$species<-as.factor(stem_mature_2020$species)
m_stem_mature_2020<-gam(percent_TH/100 ~ s(DOY_cumul,by = species, k=3)+species,  data=stem_mature_2020, family="quasibinomial", method="REML")
summary(m_stem_mature_2020)
par(mfrow = c(2, 2))
gam.check(m_stem_mature_2020)
est_stem_mature_2020<-ggpredict(m_stem_mature_2020, terms=c("DOY_cumul", "species"))
```

__Stem thickening mature trees spring 2021__
```{r echo=FALSE, message=FALSE, warning=FALSE}
stem_mature_2021$species<-as.factor(stem_mature_2021$species)
stem_mature_2021<-subset(stem_mature_2021, DOY_cumul>480)
m_stem_mature_2021<-gam(percent_TH/100 ~ s(DOY_cumul,by = species, k=3),  data=stem_mature_2021, family="quasibinomial", method="REML")
summary(m_stem_mature_2021)
par(mfrow = c(2, 2))
gam.check(m_stem_mature_2021)
est_stem_mature_2021<-ggpredict(m_stem_mature_2021, terms=c("DOY_cumul", "species"))
est_stem_mature_2021<-subset(est_stem_mature_2021, x<520)
```

## B
output of the GAM models

__coarse roots thickening saplings 2017__
```{r echo=FALSE, message=FALSE, warning=FALSE}
m_root_sapling_2017<-gam(percent_TH/100 ~ s(DOY_cumul, by= species, k=3)+species,  data=root_sapling_2017, family="quasibinomial", method="REML")
summary(m_root_sapling_2017)
par(mfrow = c(2, 2))
gam.check(m_root_sapling_2017)
est_root_sapling_2017<-ggpredict(m_root_sapling_2017, terms=c("DOY_cumul", "species"))
```

__coarse roots thickening mature trees 2020__
```{r echo=FALSE, message=FALSE, warning=FALSE}
m_root_mature_2020<-gam(percent_TH/100 ~ s(DOY_cumul,by = species, k=3)+species,  data=root_mature_2020, family="quasibinomial", method="REML")
summary(m_root_mature_2020)
par(mfrow = c(2, 2))
gam.check(m_root_mature_2020)
est_root_mature_2020<-ggpredict(m_root_mature_2020, terms=c("DOY_cumul", "species"))
```


# Pannel Fig.3

```{r echo=TRUE, message=FALSE, warning=FALSE, results='hide', echo=FALSE, message=FALSE, warning=FALSE,fig.width= 13, fig.height= 13 }
Fig3A_caption<-ggplot()+
    geom_line(data=est_stem_sapling_2017,aes(x, predicted*100,colour=group), size =.5)+
    geom_ribbon(data=est_stem_sapling_2017, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.4)+
    geom_line(data=est_senescence_sapling_2017,aes(x, predicted*40,colour=group),size =.5, linetype="dashed")+
    geom_ribbon(data=est_senescence_sapling_2017, aes(x, predicted*40,ymin = conf.low*40, ymax = conf.high*40, fill=group), alpha = 0.4)+
    geom_line(data=est_bud_sapling_2017,aes(x, predicted*40,colour=group), size= .5,linetype="dotted")+
    geom_ribbon(data=est_bud_sapling_2017, aes(x, predicted*40,ymin = conf.low*40, ymax = conf.high*40, fill=group), alpha = 0.4)+
      scale_fill_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
      scale_color_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
  xlim(210, 510)+
  labs(x=" ", y=" ")+
  scale_y_continuous(sec.axis = sec_axis(~.*2.5, name=" "))+
   theme(panel.background = element_rect(fill = "white",colour="grey"),panel.grid  = element_blank(),plot.background = element_rect(fill = "white",colour="black"), axis.text = element_text(size = 10),axis.title = element_blank(),legend.position="none",axis.text.x=element_text(angle=0, hjust=1))


section_split1<-264
section_split2<-355
section_split3<-445

title<- Fig3A_caption + annotate("rect", xmin = c(-Inf, section_split1,section_split2, section_split3 ), ymin = c(-Inf, -Inf, -Inf, -Inf),
             xmax = c(section_split1,section_split2,section_split3,Inf), ymax = c(Inf,Inf,Inf,Inf),
             fill = c("grey", "lightgrey", "grey","lightgrey"))+
  annotate("text", label = c("Summer","Autum","Winter","Spring"), size = 3, y= rep(mean(layer_scales(Fig3A_caption)$y$range$range),4),
           x = c(253,310,400,475))+
  theme_void()+ theme(legend.position = "none")

Fig3A<-ggplot()+
    geom_line(data=est_stem_mature_2020,aes(x, predicted*100,colour=group),size=1)+
    geom_ribbon(data=est_stem_mature_2020, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.4)+
    geom_line(data=est_senescence_mature_2020,aes(x, predicted*40,colour=group), linetype="dashed",size=1)+
    geom_ribbon(data=est_senescence_mature_2020, aes(x, predicted*40,ymin = conf.low*40, ymax = conf.high*40, fill=group), alpha = 0.4)+
    geom_line(data=est_bud_mature_2020,aes(x, predicted*40,colour=group), linetype="dotted",size=1)+
    geom_ribbon(data=est_bud_mature_2020, aes(x, predicted*40,ymin = conf.low*40, ymax = conf.high*40, fill=group), alpha = 0.4)+
  geom_line(data=est_stem_mature_2021,aes(x, predicted*100,colour=group),size=1)+
    geom_ribbon(data=est_stem_mature_2021, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.4)+
    scale_y_continuous(sec.axis = sec_axis(~.*2.5, name="Bud + leaf development (%) \n Leaf coloration + fall (%)"))+
    scale_fill_manual(values=c("#fc8d59", "#67a9cf"))+
    scale_color_manual(values=c("#fc8d59", "#67a9cf"))+
    xlim(210,520)+
     annotation_custom(ggplotGrob(Fig3A_caption), xmin = 308, xmax = 438, 
                       ymin = 3, ymax = 30)+
   theme(axis.text = element_text(size = 12),axis.title = element_text(size = 14), legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" Day of year", y="Proportion of thickening  \n stem ring (%)")

Fig3B_caption<-ggplot()+
      geom_ribbon(data=est_root_sapling_2017, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.2)+
  geom_line(data=est_root_sapling_2017,aes(x, predicted*100,colour=group), size =0.75)+
      scale_fill_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
      scale_color_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
  xlim(210, 510)+
  theme(panel.background = element_rect(fill = "white",colour="grey"),panel.grid  = element_blank(),plot.background = element_rect(fill = "white",colour="black"), axis.text = element_text(size = 10),axis.title = element_blank(),legend.position="none",axis.text.x=element_text(angle=0, hjust=1))

Fig3B<-ggplot()+
    geom_line(data=est_root_mature_2020,aes(x, predicted*100,colour=group), size=1)+
    geom_ribbon(data=est_root_mature_2020, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.4)+
      scale_fill_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
      scale_color_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
  xlim(210,520)+
  ylim(0,40)+
   annotation_custom(ggplotGrob(Fig3B_caption), xmin = 307, xmax = 428, 
                       ymin = 10, ymax = 40)+
   theme(axis.text = element_text(size = 12),axis.title = element_text(size = 14), axis.text.x=element_text(angle=0, hjust=1),legend.position="none",panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" Day of year", y=" Proportion of thickening \n coarse root ring (%)")


Fig3C<-ggplot()+
  geom_line(data=est_fine_root,aes(x, predicted,colour=group), size=1)+
    geom_ribbon(data=est_fine_root, aes(x, predicted,ymin = conf.low, ymax = conf.high, fill=group), alpha = 0.4)+
   scale_fill_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
    scale_color_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
  xlim(210,520)+
   theme(axis.text = element_text(size = 12),axis.title = element_text(size = 14),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" Day of year Belgium", y=" Cumulative fine \n roots surface (mm²)")

break.vec <- c(as.Date("2020-07-18"),as.Date("2020-10-26"),as.Date("2021-02-04"),as.Date("2021-05-15"))
climate_mature_2020$date<-as.Date(climate_mature_2020$date)
Fig3_D<-ggplot(climate_mature_2020, aes(x = date,))+
  geom_ribbon(aes(ymin = min_TA, ymax = max_TA), fill = "#f1b6da")+
  geom_ribbon(aes(ymin = min_TS, ymax = max_TS), fill = "#80cdc1")+
  scale_color_manual(values=c("#d01c8b",'#018571'))+
  geom_line(aes(y = mean_TA, color = "mean ait temperature"),size=1)+
  geom_line(aes(y = mean_TS, color = "mean soil temperature"),size=1 )+
  geom_hline(yintercept=5, linetype="dashed", color="grey")+
  scale_x_date( breaks = break.vec , date_labels = "%j", limits=c(as.Date("2020-07-28"), as.Date("2021-06-04")))+
 theme(axis.text = element_text(size = 12),axis.title = element_text(size = 14),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" DOY", y="Soil and air \n temperature (ºC) ")

header <-title+ theme(plot.margin = unit(c(0, 1.8, 0, 2.3), "cm"))
Fig3A_bis = Fig3A + theme(plot.margin = unit(c(0, 0, -1, 0.5), "cm"))
Fig3B_bis = Fig3B + theme(plot.margin = unit(c(0.5, 1.9, -1,0.5), "cm"))
Fig3C_bis = Fig3C + theme(plot.margin = unit(c(0.5, 1.9, -1,0.5), "cm"))
Fig3_D_bis = Fig3_D + theme(plot.margin = unit(c(0.5, 1.9, 0, 0.2), "cm"))

plot_fig_3<-plot_grid(header,Fig3A_bis,Fig3B_bis,Fig3C_bis,Fig3_D_bis, nrow=5,labels=c(" ",'A', 'B', 'C',"D"),rel_heights=c(0.7,8,8,6,7),  label_size = 20);plot_fig_3

```

# Analyse Fig. 4
## A
output of the GAM models

__Senescence mature trees 2018__
```{r echo=FALSE, message=FALSE, warning=FALSE}
m_senescence_mature_2018<-gam(senescence/100 ~ s(DOY, k=3),  data=senescence_mature_2018, family="quasibinomial", method="REML")
summary(m_senescence_mature_2018)
est_senescence_mature_2018<-ggpredict(m_senescence_mature_2018, terms=c("DOY"))
```

__Bud development mature trees 2018__
```{r echo=FALSE, message=FALSE, warning=FALSE}
m_bud_mature_2018<-gam(score/4 ~ s(DOY_cumul, k=3),  data=bud_mature_2018, family="quasibinomial", method="REML")
summary(m_bud_mature_2018)
est_bud_mature_2018<-ggpredict(m_bud_mature_2018, terms=c("DOY_cumul"))
```

__Stem thickening mature trees 2018__
```{r echo=FALSE, message=FALSE, warning=FALSE}
m_stem_mature_2018<-gam(percent_TH/100 ~ s(DOY_cumul, k=3),  data=stem_mature_2018, family="quasibinomial", method="REML")
summary(m_stem_mature_2018)
par(mfrow = c(2, 2))
gam.check(m_stem_mature_2018)
est_stem_mature_2018<-ggpredict(m_stem_mature_2018, terms=c("DOY_cumul"))
```

__Stem thickening mature trees spring 2019__
```{r echo=FALSE, message=FALSE, warning=FALSE}
m_stem_mature_2019<-gam(percent_TH/100 ~ s(DOY_cumul,by = species, k=3),  data=stem_mature_2019, family="quasibinomial", method="REML")
summary(m_stem_mature_2019)
par(mfrow = c(2, 2))
gam.check(m_stem_mature_2019)
est_stem_mature_2019<-ggpredict(m_stem_mature_2019, terms=c("DOY_cumul", "species"))

```

## B
output of the GAM models

__Coarse root thickening mature trees 2018__
```{r echo=FALSE, message=FALSE, warning=FALSE}
m_root_mature_2018<-gam(percent_TH/100 ~ s(DOY_cumul, k=3),  data=root_mature_2020, family="quasibinomial", method="REML")
summary(m_root_mature_2018)
par(mfrow = c(2, 2))
gam.check(m_root_mature_2018)
est_root_mature_2018<-ggpredict(m_root_mature_2018, terms=c("DOY_cumul"))
```

# Pannel Fig.4

```{r echo=FALSE, message=FALSE, warning=FALSE,fig.width= 13, fig.height= 10 }
fig4A<-ggplot()+
    geom_line(data=est_stem_mature_2018,aes(x, predicted*100),color="#67a9cf", size=1)+
    geom_ribbon(data=est_stem_mature_2018, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100),fill="#67a9cf", alpha = 0.4)+
    geom_line(data=est_senescence_mature_2018,aes(x, predicted*40),color="#67a9cf", linetype="dashed", size=1)+
    geom_ribbon(data=est_senescence_mature_2018, aes(x, predicted*40,ymin = conf.low*40, ymax = conf.high*40), fill="#67a9cf",alpha = 0.1)+
    geom_line(data=est_bud_mature_2018,aes(x, predicted*40), color="#67a9cf",linetype="dotted", size=1)+
    geom_ribbon(data=est_bud_mature_2018, aes(x, predicted*40,ymin = conf.low*40, ymax = conf.high*40), fill="#67a9cf",alpha = 0.1)+
    geom_line(data=est_stem_mature_2019,aes(x, predicted*100),color="#67a9cf", size=1)+
    geom_ribbon(data=est_stem_mature_2019, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100),fill="#67a9cf", alpha = 0.4)+
    scale_y_continuous(sec.axis = sec_axis(~.*2.5, name="Bud + leaf develoment (%) \n Leaf coloration + fall (%)"))+
    scale_fill_manual(values=c("black"))+
    scale_color_manual(values=c("black"))+
    xlim(210,520)+
   theme(axis.text = element_text(size = 12),axis.title = element_text(size = 14),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" Day of year", y="Proportion of thickening  \n stem ring (%)")

fig4B<-ggplot()+
    geom_line(data=est_root_mature_2018,aes(x, predicted*100),color="#67a9cf", size=1)+
    geom_ribbon(data=est_root_mature_2018, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100),fill="#67a9cf", alpha = 0.4)+
  xlim(210,520)+
  ylim(0,40)+
   theme(axis.text = element_text(size = 12),axis.title = element_text(size = 14),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" Day of year Belgium", y="  Proportion of thickening  \n coarse root ring (%)")

break.vec <- c(as.Date("2018-07-19"),as.Date("2018-10-27"),as.Date("2019-02-04"),as.Date("2019-05-15"))

climate_mature_2018$date<-as.Date(climate_mature_2018$date)

fig4C<-ggplot(climate_mature_2018, aes(x = date,))+
  geom_ribbon(aes(ymin = min_TA, ymax = max_TA), fill = "#f1b6da")+
    geom_ribbon(aes(ymin = min_TS, ymax = max_TS), fill = "#80cdc1")+
    scale_color_manual(values=c("#d01c8b",'#018571'))+
  geom_line(aes(y = mean_TA, color = "mean ait temperature"),size=1)+
  geom_line(aes(y = mean_TS, color = "mean soil temperature"),size=1)+
  geom_hline(yintercept=5, linetype="dashed", color="grey")+
   scale_x_date( breaks = break.vec , date_labels = "%j", limits=c(as.Date("2018-07-29"), as.Date("2019-06-04")))+
 theme(axis.text = element_text(size = 12),axis.title = element_text(size = 14),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" DOY", y="Soil and air \n temperature (ºC) ")

section_split1<-264
section_split2<-355
section_split3<-445

title<- fig4B + annotate("rect", xmin = c(-Inf, section_split1,section_split2, section_split3 ), ymin = c(-Inf, -Inf, -Inf, -Inf),
             xmax = c(section_split1,section_split2,section_split3,Inf), ymax = c(Inf,Inf,Inf,Inf),
             fill = c("grey", "lightgrey", "grey","lightgrey"))+
  annotate("text", label = c("Summer","Autumn","Winter","Spring"), size = 5.5, y= rep(mean(layer_scales(fig4B)$y$range$range),4),
           x = c(230,310,400,475))+
  theme_void()+ theme(legend.position = "none")

header <-title+ theme(plot.margin = unit(c(0, 2, 0, 2.3), "cm"))
fig4A_bis = fig4A + theme(plot.margin = unit(c(0, 0, -1, 0.5), "cm"))
fig4B_bis =  fig4B+ theme(plot.margin = unit(c(0.5, 2, -1,0.5), "cm"))
fig4C_bis = fig4C + theme(plot.margin = unit(c(0.5, 2, 0, 0.5), "cm"))

fig4<-plot_grid(header, fig4A_bis,fig4B_bis,fig4C_bis, nrow=4,labels=c(" ",'A', 'B', 'C'), rel_heights=c(0.8,6,6,5), label_size = 20 );fig4
```

# Analyse fig. 5
## A
```{r echo=FALSE, message=FALSE, warning=FALSE}
temp_TH_mature_2020_stem<- subset(temp_TH_mature_2020, organ=="stem")
temp_TH_mature_2020_stem$species<-as.factor(temp_TH_mature_2020_stem$species)
temp_TH_mature_2020_root<- subset(temp_TH_mature_2020, organ=="roots")
temp_TH_mature_2020_root$species<-as.factor(temp_TH_mature_2020_root$species)
```

output og GAM models 

__relation coarse roots thickening and soil temperature mature trees 2020__
```{r echo=FALSE, message=FALSE, warning=FALSE}
m_temp_TH_mature_2020_root<-gam(mean_persent_TH/100 ~ s(X10days_temp,by = species, k=3)+species,  data=temp_TH_mature_2020_root, family="quasibinomial", method="REML")
summary(m_temp_TH_mature_2020_root)
par(mfrow = c(2, 2))
gam.check(m_temp_TH_mature_2020_root)
est_temp_TH_mature_2020_root<-ggpredict(m_temp_TH_mature_2020_root, terms=c("X10days_temp", "species"))
```

__relation stems thickening and air temperature mature trees 2020__
```{r echo=FALSE, message=FALSE, warning=FALSE}
m_temp_TH_mature_2020_stem<-gam(mean_persent_TH/100 ~ s(X10days_temp,by = species, k=3)+species,  data=temp_TH_mature_2020_stem, family="quasibinomial", method="REML")
summary(m_temp_TH_mature_2020_stem)
par(mfrow = c(2, 2))
gam.check(m_temp_TH_mature_2020_stem)
est_temp_TH_mature_2020_stem<-ggpredict(m_temp_TH_mature_2020_stem, terms=c("X10days_temp", "species"))

```

## B
```{r echo=FALSE, message=FALSE, warning=FALSE}
temp_TH_sampling_2017_stem<- subset(temp_TH_sampling_2017, organ=="stem")
temp_TH_sampling_2017_stem$species<-as.factor(temp_TH_sampling_2017_stem$species)
temp_TH_sampling_2017_root<- subset(temp_TH_sampling_2017, organ=="roots")
temp_TH_sampling_2017_root$species<-as.factor(temp_TH_sampling_2017_root$species)
```

__relation coarse roots thickening and soil temperature saplings 2017__
```{r echo=FALSE, message=FALSE, warning=FALSE}
m_temp_TH_sampling_2017_root<-gam(mean_persent_TH/100 ~ s(X10days_temp,by = species, k=3)+species,  data=temp_TH_sampling_2017_root, family="quasibinomial", method="REML")
summary(m_temp_TH_sampling_2017_root)
par(mfrow = c(2, 2))
gam.check(m_temp_TH_sampling_2017_root)
est_temp_TH_sampling_2017_root<-ggpredict(m_temp_TH_sampling_2017_root, terms=c("X10days_temp", "species"))
```

__relation stems thickening and air temperature saplings 2017__
```{r echo=FALSE, message=FALSE, warning=FALSE}
m_temp_TH_sampling_2017_stem<-gam(mean_persent_TH/100 ~ s(X10days_temp,by = species, k=3)+species,  data=temp_TH_sampling_2017_stem, family="quasibinomial", method="REML")
summary(m_temp_TH_sampling_2017_stem)
par(mfrow = c(2, 2))
gam.check(m_temp_TH_sampling_2017_stem)
est_temp_TH_sampling_2017_stem<-ggpredict(m_temp_TH_sampling_2017_stem, terms=c("X10days_temp", "species"))
```

# Pannel 5

```{r echo=FALSE, message=FALSE, warning=FALSE, fig.width= 7, fig.height= 6}
Fig5A<-ggplot()+
    geom_ribbon(data=est_temp_TH_mature_2020_root, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.4)+
   geom_line(data=est_temp_TH_mature_2020_root,aes(x, predicted*100,colour=group), size=1,linetype="dashed")+
    geom_ribbon(data=est_temp_TH_mature_2020_stem, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.4)+
   geom_line(data=est_temp_TH_mature_2020_stem,aes(x, predicted*100,colour=group), size=1)+
      scale_fill_manual(values=c("#fc8d59", "#67a9cf"))+
      scale_color_manual(values=c("#fc8d59", "#67a9cf"))+
  xlim(2,20)+
  ylim(0,40)+
   theme(axis.text = element_text(size = 12),axis.title = element_text(size = 12),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x="    ", y="Proportion of \n thickening ring (%)")


Fig5B<-ggplot()+
    geom_ribbon(data=est_temp_TH_sampling_2017_stem, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.2)+
  geom_line(data=est_temp_TH_sampling_2017_stem,aes(x, predicted*100,colour=group), size=1)+
    geom_ribbon(data=est_temp_TH_sampling_2017_root, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.2)+
    geom_line(data=est_temp_TH_sampling_2017_root,aes(x, predicted*100,colour=group),size=1, linetype="dashed")+
      scale_fill_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
      scale_color_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
      xlim(2,20)+
  ylim(0,40)+
   theme(axis.text = element_text(size = 12),axis.title = element_text(size = 12),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" Mean 10 days temperature prior sampling (°C) ", y="Proportion of \n thickening ring (%)")

Fig5A_bis<-Fig5A + theme(plot.margin = unit(c(0.5, 0.5, 0, 0.5), "cm"))
Fig5B_bis<-Fig5B + theme(plot.margin = unit(c(0, 0.5, 0, 0.5), "cm"))

Fig5<-plot_grid(Fig5A_bis,Fig5B_bis, nrow=2, rel_heights=c( 1,1),labels=c('A', 'B'));Fig5

```

# Analyse Fig. 6 
## A
output of GAM models 

__Senecence saplings Belgium 2018__
```{r echo=FALSE, message=FALSE, warning=FALSE}
senescence_sapling_2018_belgium$species<-as.factor(senescence_sapling_2018_belgium$species)
m_senescence_sapling_2018_belgium<-gam(senescence/100 ~ s(DOY,by = species, k=3),  data=senescence_sapling_2018_belgium, family="quasibinomial", method="REML")
summary(m_senescence_sapling_2018_belgium)
est_senescence_sapling_2018_belgium<-ggpredict(m_senescence_sapling_2018_belgium, terms=c("DOY", "species"))
```

__Stem thickening saplings Belgium 2018__
```{r echo=FALSE, message=FALSE, warning=FALSE}
stem_sapling_2018_belgium$species<-as.factor(stem_sapling_2018_belgium$species)
m_stem_sapling_2018_belgium<-gam(percent_TH/100 ~ s(DOY_cumul,by = species, k=3),  data=stem_sapling_2018_belgium, family="quasibinomial", method="REML")
summary(m_stem_sapling_2018_belgium)
est_stem_sapling_2018_belgium<-ggpredict(m_stem_sapling_2018_belgium, terms=c("DOY_cumul", "species"))
```

## B
output of GAM models 

__Senecence saplings Spain 2018__
```{r echo=FALSE, message=FALSE, warning=FALSE}
senescence_sapling_2018_spain$species<-as.factor(senescence_sapling_2018_spain$species)
m_senescence_sapling_2018_spain<-gam(senescence/100 ~ s(DOY,by = species, k=3),  data=senescence_sapling_2018_spain, family="quasibinomial", method="REML")
summary(m_senescence_sapling_2018_spain)
est_senescence_sapling_2018_spain<-ggpredict(m_senescence_sapling_2018_spain, terms=c("DOY", "species"))
```

__Stem thickening saplings Spain 2018__
```{r echo=FALSE, message=FALSE, warning=FALSE}
stem_sapling_2018_spain$species<-as.factor(stem_sapling_2018_spain$species)
m_stem_sapling_2018_spain<-gam(percent_TH/100 ~ s(DOY_cumul,by = species, k=3),  data=stem_sapling_2018_spain, family="quasibinomial", method="REML")
summary(m_stem_sapling_2018_spain)
est_stem_sapling_2018_spain<-ggpredict(m_stem_sapling_2018_spain, terms=c("DOY_cumul", "species"))
```

## C
output of GAM models 

__Senecence saplings Norway 2018__
```{r echo=FALSE, message=FALSE, warning=FALSE}
senescence_sapling_2018_norway$species<-as.factor(senescence_sapling_2018_norway$species)
m_senescence_sapling_2018_norway<-gam(senescence/100 ~ s(DOY,by = species, k=3),  data=senescence_sapling_2018_norway, family="quasibinomial", method="REML")
summary(m_senescence_sapling_2018_norway)
est_senescence_sapling_2018_norway<-ggpredict(m_senescence_sapling_2018_norway, terms=c("DOY", "species"))
```

__Stem thickening saplings Norway 2018__
```{r echo=FALSE, message=FALSE, warning=FALSE}
stem_sapling_2018_norway$species<-as.factor(stem_sapling_2018_norway$species)
m_stem_sapling_2018_norway<-gam(percent_TH/100 ~ s(DOY_cumul,by = species, k=3),  data=stem_sapling_2018_norway, family="quasibinomial", method="REML")
summary(m_stem_sapling_2018_norway)
est_stem_sapling_2018_norway<-ggpredict(m_stem_sapling_2018_norway, terms=c("DOY_cumul", "species"))
```

## D
__Coarse root thickening saplings Belgium 2018__
```{r echo=FALSE, message=FALSE, warning=FALSE}
m_root_sapling_2018_belgium<-gam(percent_TH/100 ~ s(DOY_cumul, by= species, k=3)+species,  data=root_sapling_2018_belgium, family="quasibinomial", method="REML")
summary(m_root_sapling_2018_belgium)
par(mfrow = c(2, 2))
gam.check(m_root_sapling_2018_belgium)
est_root_sapling_2018_belgium<-ggpredict(m_root_sapling_2018_belgium, terms=c("DOY_cumul", "species"))
```

## E
__Coarse root thickening saplings Spain 2018__
```{r echo=FALSE, message=FALSE, warning=FALSE}
m_root_sapling_2018_spain<-gam(percent_TH/100 ~ s(DOY_cumul, by= species, k=3)+species,  data=root_sapling_2018_spain, family="quasibinomial", method="REML")
summary(m_root_sapling_2018_spain)
par(mfrow = c(2, 2))
gam.check(m_root_sapling_2018_spain)
est_root_sapling_2018_spain<-ggpredict(m_root_sapling_2018_spain, terms=c("DOY_cumul", "species"))
```

## F
__Coarse root thickening saplings Norway 2018__
```{r echo=FALSE, message=FALSE, warning=FALSE}
m_root_sapling_2018_norway<-gam(percent_TH/100 ~ s(DOY_cumul, by= species, k=3)+species,  data=root_sapling_2018_norway, family="quasibinomial", method="REML")
summary(m_root_sapling_2018_norway)
par(mfrow = c(2, 2))
gam.check(m_root_sapling_2018_norway)
est_root_sapling_2018_norway<-ggpredict(m_root_sapling_2018_norway, terms=c("DOY_cumul", "species"))
```

# Pannel 6 

```{r echo=FALSE, message=FALSE, warning=FALSE,fig.width= 13, fig.height= 10}
Fig6A<-ggplot()+
    geom_line(data=est_stem_sapling_2018_belgium,aes(x, predicted*100,colour=group), size=1)+
    geom_ribbon(data=est_stem_sapling_2018_belgium, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.4)+
    geom_line(data=est_senescence_sapling_2018_belgium,aes(x, predicted*80,colour=group), linetype="dashed", size=1)+
    geom_ribbon(data=est_senescence_sapling_2018_belgium, aes(x, predicted*80,ymin = conf.low*80, ymax = conf.high*80, fill=group), alpha = 0.4)+
    scale_y_continuous(sec.axis = sec_axis(~.*1.25, name=" "))+
    scale_fill_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
    scale_color_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
  xlim(220,335)+
   theme(axis.text = element_text(size = 12),axis.title = element_text(size = 12),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" ", y=" Proportion of thickening \n stem ring (%)")

Fig6B<-ggplot()+
    geom_line(data=est_stem_sapling_2018_spain,aes(x, predicted*100,colour=group), size=1)+
    geom_ribbon(data=est_stem_sapling_2018_spain, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.4)+
    geom_line(data=est_senescence_sapling_2018_spain,aes(x, predicted*80,colour=group), linetype="dashed", size=1)+
    geom_ribbon(data=est_senescence_sapling_2018_spain, aes(x, predicted*80,ymin = conf.low*80, ymax = conf.high*80, fill=group), alpha = 0.4)+
    scale_y_continuous(sec.axis = sec_axis(~.*1.25, name=" "))+
    scale_fill_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
    scale_color_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
  xlim(220,335)+
   theme(axis.text = element_text(size = 12),axis.title = element_text(size = 12),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" ", y=" ")

Fig6C<-ggplot()+
    geom_line(data=est_stem_sapling_2018_norway,aes(x, predicted*100,colour=group), size=1)+
    geom_ribbon(data=est_stem_sapling_2018_norway, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.4)+
    geom_line(data=est_senescence_sapling_2018_norway,aes(x, predicted*80,colour=group), linetype="dashed", size=1)+
    geom_ribbon(data=est_senescence_sapling_2018_norway, aes(x, predicted*80,ymin = conf.low*80, ymax = conf.high*80, fill=group), alpha = 0.4)+
    scale_y_continuous(sec.axis = sec_axis(~.*1.25, name=" Leaf coloration + fall "))+
    scale_fill_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
    scale_color_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
  xlim(220,335)+
   theme(axis.text = element_text(size = 12),axis.title = element_text(size = 12),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" ", y=" ")

Fig6D<-ggplot()+
    geom_line(data=est_root_sapling_2018_belgium,aes(x, predicted*100,colour=group),size =1)+
    geom_ribbon(data=est_root_sapling_2018_belgium, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.4)+
      scale_fill_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
      scale_color_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
  xlim(220,335)+
   ylim(0,50)+
   theme(axis.text = element_text(size = 12),axis.title = element_text(size = 12),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" ", y=" Proportion of thickening \n coarse root ring (%)")

Fig6E<-ggplot()+
    geom_line(data=est_root_sapling_2018_spain,aes(x, predicted*100,colour=group),size =1)+
    geom_ribbon(data=est_root_sapling_2018_spain, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.4)+
      scale_fill_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
      scale_color_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
  xlim(220,335)+
   ylim(0,50)+
   theme(axis.text = element_text(size = 12),axis.title = element_text(size = 12),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" ", y="")

Fig6F<-ggplot()+
    geom_line(data=est_root_sapling_2018_norway,aes(x, predicted*100,colour=group),size =1)+
    geom_ribbon(data=est_root_sapling_2018_norway, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.4)+
      scale_fill_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
      scale_color_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
  xlim(220,335)+
   ylim(0,50)+
   theme(axis.text = element_text(size = 12),axis.title = element_text(size = 12),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" ", y="")

climate_sapling_2018_belgium$date<-as.Date(climate_sapling_2018_belgium$date)
break.vec <- c(as.Date("2018-08-28"),as.Date("2018-09-27"),as.Date("2018-10-27"),as.Date("2018-11-26"))
Fig6G<-ggplot()+
  geom_line( aes(x= date, y = mean_TA), data=climate_sapling_2018_belgium,size=1, color="#d01c8b" )+
    scale_x_date( breaks = break.vec , date_labels = "%j", limits=c(as.Date("2018-08-08"), as.Date("2018-12-01")))+
  ylim(-10,25)+
  geom_hline(yintercept=5, linetype="dashed", color="grey")+
 theme(axis.text = element_text(size = 12),axis.title = element_text(size = 12),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" ", y=" Air \n temperature (ºC)")

Fig6H<-ggplot()+
  geom_line( aes(x= DOY, y = mean_TA), data=climate_sapling_2018_spain,size=1, color="#d01c8b" )+
  ylim(-10,25)+
    xlim(220,335)+
  geom_hline(yintercept=5, linetype="dashed", color="grey")+
 theme(axis.text = element_text(size = 12),axis.title = element_text(size = 12),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x="DOY Autumn 2018 ", y=" Air \n temperature (ºC)")

Fig6I<-ggplot()+
  geom_line( aes(x= DOY, y = mean_TA), data=climate_sapling_2018_norway,size=1, color="#d01c8b" )+
  ylim(-10,25)+
    xlim(220,335)+
  geom_hline(yintercept=5, linetype="dashed", color="grey")+
 theme(axis.text = element_text(size = 12),axis.title = element_text(size = 12),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" ", y=" Air \n temperature (ºC)")

Fig6A_bis = Fig6A + theme(plot.margin = unit(c(0.8, -0.4, -1, 0.2), "cm"))
Fig6B_bis = Fig6B + theme(plot.margin = unit(c(0.8,-0.4, -1, -0.3), "cm"))
Fig6C_bis = Fig6C + theme(plot.margin = unit(c(0.8,0, -1,-0.2), "cm"))

Fig6D_bis = Fig6D + theme(plot.margin = unit(c(0.8, 0.8, -1,0.2), "cm"))
Fig6E_bis = Fig6E + theme(plot.margin = unit(c(0.8, 0.8, -1,-0.7), "cm"))
Fig6F_bis = Fig6F + theme(plot.margin = unit(c(0.8, 1.2, -1,-0.6), "cm"))

Fig6G_bis = Fig6G + theme(plot.margin = unit(c(0.8, 0.8, 0,0), "cm"))
Fig6H_bis = Fig6H + theme(plot.margin = unit(c(0.8, 0.8, 0,-0.8), "cm"))
Fig6I_bis = Fig6I + theme(plot.margin = unit(c(0.8, 1.2, 0,-0.7), "cm"))

column1<- plot_grid(Fig6A_bis, Fig6D_bis, Fig6G_bis, nrow=3, labels=c('A','D','G'))
column2<- plot_grid(Fig6B_bis, Fig6E_bis, Fig6H_bis, nrow=3, labels=c('B','E','H'))
column3<- plot_grid(Fig6C_bis, Fig6F_bis, Fig6I_bis, nrow=3, labels=c('C','F','I'))

column1_bis<-column1 + theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))
column2_bis<-column2 + theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))
column3_bis<-column3 + theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))

Fig6<-plot_grid(column1_bis,column2_bis,column3_bis, ncol=3,rel_widths = c(1.1,1,1));Fig6
```

# SEPPLEMENTARY
# Analyse supp data Fig. 4

## A
Output GAM models

__Coarse root enlargement mature trees 2020__
```{r echo=FALSE, message=FALSE, warning=FALSE}
m_root_mature_2020_EC<-gam(percent_EC/100 ~ s(DOY_cumul,by = species, k=3)+species,  data=root_mature_2020, family="quasibinomial", method="REML")
summary(m_root_mature_2020_EC)
par(mfrow = c(2, 2))
gam.check(m_root_mature_2020_EC)
est_root_mature_2020_EC<-ggpredict(m_root_mature_2020_EC, terms=c("DOY_cumul", "species"))
```

## B
__Coarse root enlargement mature trees 2018__
```{r echo=FALSE, message=FALSE, warning=FALSE}
m_root_mature_2018_EC<-gam(percent_EC/100 ~ s(DOY_cumul, k=3),  data=root_mature_2018, family="quasibinomial", method="REML")
summary(m_root_mature_2018_EC)
par(mfrow = c(2, 2))
gam.check(m_root_mature_2018_EC)
est_root_mature_2018_EC<-ggpredict(m_root_mature_2018_EC, terms=c("DOY_cumul"))
```

## C
__Coarse root enlargement saplings 2017__
```{r echo=FALSE, message=FALSE, warning=FALSE}
m_root_sapling_2017_EC<-gam(percent_EC/100 ~ s(DOY_cumul,by = species, k=3)+species,  data=root_sapling_2017, family="quasibinomial", method="REML")
summary(m_root_sapling_2017_EC)
par(mfrow = c(2, 2))
gam.check(m_root_sapling_2017_EC)
est_root_sapling_2017_EC<-ggpredict(m_root_sapling_2017_EC, terms=c("DOY_cumul", "species"))
```

## D
__Coarse root enlargement saplings Belgium 2018__
```{r echo=FALSE, message=FALSE, warning=FALSE}
root_sapling_2018_belgium$EC_poisson<-round(root_sapling_2018_belgium$percent_EC, 0)
m_root_sapling_2018_belgium_EC<-gam(EC_poisson ~ s(DOY_cumul, by= species, k=3),  data=root_sapling_2018_belgium, family="nb", method="REML")
summary(m_root_sapling_2018_belgium_EC)
par(mfrow = c(2, 2))
gam.check(m_root_sapling_2018_belgium_EC)
est_root_sapling_2018_belgium_EC<-ggpredict(m_root_sapling_2018_belgium_EC, terms=c("DOY_cumul", "species"))
```

## E
__Coarse root enlargement saplings Spain 2018__
```{r echo=FALSE, message=FALSE, warning=FALSE}
m_root_sapling_2018_spain_EC<-gam(percent_EC/100 ~ s(DOY_cumul,by = species, k=3)+species,  data=root_sapling_2018_spain, family="quasibinomial", method="REML")
summary(m_root_sapling_2018_spain_EC)
par(mfrow = c(2, 2))
gam.check(m_root_sapling_2018_spain_EC)
est_root_sapling_2018_spain_EC<-ggpredict(m_root_sapling_2018_spain_EC, terms=c("DOY_cumul", "species"))
```

## F
__Coarse root enlargement saplings Norway 2018__
```{r echo=FALSE, message=FALSE, warning=FALSE}
m_root_sapling_2018_norway_EC<-gam(percent_EC/100 ~ s(DOY_cumul,by = species, k=3)+species,  data=root_sapling_2018_norway, family="quasibinomial", method="REML")
summary(m_root_sapling_2018_norway_EC)
par(mfrow = c(2, 2))
gam.check(m_root_sapling_2018_norway_EC)
est_root_sapling_2018_norway_EC<-ggpredict(m_root_sapling_2018_norway_EC, terms=c("DOY_cumul", "species"))
```

# Panel supp data 2

```{r echo=FALSE, message=FALSE, warning=FALSE, ,fig.height= 9, fig.width= 9}

SDF2A<-ggplot()+
    geom_line(data=est_root_mature_2020,aes(x, predicted*100,colour=group), size=1)+
    geom_ribbon(data=est_root_mature_2020, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.4)+
    geom_point(data= root_mature_2020, aes(DOY_cumul,percent_TH, colour=species))+
      scale_fill_manual(values=c("#fc8d59", "#67a9cf"))+
      scale_color_manual(values=c("#fc8d59", "#67a9cf"))+
  xlim(210,520)+
  ylim(0,80)+
   theme(axis.text = element_text(size = 10),axis.title = element_text(size = 10), axis.text.x=element_text(angle=0, hjust=1),legend.position="none",panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" Day of year", y=" Proportion of thickening \n coarse root ring (%)")

SDF2B<-ggplot()+
    geom_line(data=est_root_mature_2020_EC,aes(x, predicted*100,colour=group),size=1)+
    geom_ribbon(data=est_root_mature_2020_EC, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.4)+
    geom_point(data= root_mature_2020, aes(DOY_cumul,percent_EC, colour=species))+
      scale_fill_manual(values=c("#fc8d59", "#67a9cf"))+
      scale_color_manual(values=c("#fc8d59", "#67a9cf"))+
  xlim(210,520)+
  ylim(0,50)+
   theme(axis.text = element_text(size = 10),axis.title = element_text(size = 10),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" Day of year ", y=" Proportion of enlarging \n coarse root ring (%)")

SDF2C<-ggplot()+
    geom_line(data=est_root_mature_2018,aes(x, predicted*100), color="#67a9cf", size=1)+
    geom_ribbon(data=est_root_mature_2018, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100), fill="#67a9cf",alpha = 0.4)+
    geom_point(data=root_mature_2018, aes(DOY_cumul,percent_TH),color="#67a9cf")+
  xlim(210,520)+
  ylim(0,80)+
   theme(axis.text = element_text(size = 10),axis.title = element_text(size = 10),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" Day of year Belgium", y="  Proportion of thickening  \n coarse root ring (%)")

SDF2D<-ggplot()+
    geom_line(data=est_root_mature_2018_EC,aes(x, predicted*100), color="#67a9cf", size=1)+
    geom_ribbon(data=est_root_mature_2018_EC, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100), fill="#67a9cf",alpha = 0.4)+
    geom_point(data=root_mature_2018, aes(DOY_cumul,percent_EC),color="#67a9cf")+
  xlim(210,520)+
  ylim(0,50)+
   theme(axis.text = element_text(size = 10),axis.title = element_text(size = 10),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" Day of year Belgium", y="  Proportion of enlarging  \n coarse root ring (%)")

SDF2E<-ggplot()+
      geom_ribbon(data=est_root_sapling_2017, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.2)+
      geom_line(data=est_root_sapling_2017,aes(x, predicted*100,colour=group), size =0.75)+
      geom_point(data=root_sapling_2017, aes(DOY_cumul, percent_TH,colour=species ))+
      scale_fill_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
      scale_color_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
  xlim(210, 510)+
  ylim(0,80)+
  theme(axis.text = element_text(size = 11),axis.title = element_text(size = 11), axis.text.x=element_text(angle=0, hjust=1),legend.position="none",panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ 
  xlab("Day of year") + 
  ylab("Proportion of thickening \n  coarse root ring (%)")

SDF2F<-ggplot()+
      geom_ribbon(data=est_root_sapling_2017_EC, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.2)+
      geom_line(data=est_root_sapling_2017_EC,aes(x, predicted*100,colour=group), size =0.75)+
      geom_point(data=root_sapling_2017, aes(DOY_cumul, percent_EC,colour=species ))+
      scale_fill_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
      scale_color_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
  xlim(210, 510)+
  ylim(0,50)+
  theme(axis.text = element_text(size = 11),axis.title = element_text(size = 11), axis.text.x=element_text(angle=0, hjust=1),legend.position="none",panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ 
  xlab("Day of year") + 
  ylab("Proportion of enlarging \n  coarse root ring (%)")



supp_data_fig2<-plot_grid(SDF2A,SDF2B,SDF2C,SDF2D,SDF2E,SDF2F, nrow=3, ncol=2, labels=c('A', 'B', 'C','D', 'E', 'F'));supp_data_fig2

```


# Analyse supp data fig. 3

## A
```{r echo=FALSE, message=FALSE, warning=FALSE, ,fig.height= 3, fig.width= 6}
fit<-rpart(percent_EC ~ species+DOY_cumul, data=root_mature_2020)
pfit<-prune(fit, cp= 0.01)
rpart.plot(pfit, fallen.leaves=TRUE, tweak=1, gap=0, space=0.5, type=3,clip.right.labs = FALSE,suffix="%", faclen = 0, clip.facs = TRUE, varlen=3,  extra=0, box.palette=c("white"))
```

## B
```{r echo=FALSE, message=FALSE, warning=FALSE, ,fig.height= 3, fig.width= 6}
fit<-rpart(percent_TH ~ species+DOY_cumul, data=root_mature_2020)
pfit<-prune(fit, cp= 0.01)
rpart.plot(pfit, fallen.leaves=TRUE, tweak=1, gap=0, space=0.5, type=3,clip.right.labs = FALSE,suffix="%", faclen = 0, clip.facs = TRUE, varlen=3,  extra=0, box.palette=c("white"))
```

## C
```{r echo=FALSE, message=FALSE, warning=FALSE, ,fig.height= 3, fig.width= 6}
fit<-rpart(percent_EC ~ species+DOY_cumul, data=root_sapling_2017)
pfit<-prune(fit, cp= 0.01)
rpart.plot(pfit, fallen.leaves=TRUE, tweak=1, gap=0, space=0.5, type=3,clip.right.labs = FALSE,suffix="%", faclen = 0, clip.facs = TRUE, varlen=3,  extra=0, box.palette=c("white"))
```

## D
```{r echo=FALSE, message=FALSE, warning=FALSE, ,fig.height= 3, fig.width= 6}
fit<-rpart(percent_TH ~ species+DOY_cumul, data=root_sapling_2017)
pfit<-prune(fit, cp= 0.01)
rpart.plot(pfit, fallen.leaves=TRUE, tweak=1, gap=0, space=0.5, type=3,clip.right.labs = FALSE,suffix="%", faclen = 0, clip.facs = TRUE, varlen=3,  extra=0, box.palette=c("white"))
```

## E
```{r echo=FALSE, message=FALSE, warning=FALSE, ,fig.height= 3, fig.width= 6}
root_sapling_2018<-subset(exp_2018)
fit<-rpart(percent_EC ~ species+DOY_cumul+country, data=exp_2018)
pfit<-prune(fit, cp= 0.01)
rpart.plot(pfit, fallen.leaves=TRUE, tweak=1, gap=0, space=0.5, type=3,clip.right.labs = FALSE,suffix="%", faclen = 0, clip.facs = TRUE, varlen=3,  extra=0, box.palette=c("white"))
```

## F
```{r echo=FALSE, message=FALSE, warning=FALSE, ,fig.height= 4, fig.width= 6}

fit<-rpart(percent_TH ~ species+DOY_cumul+country, data=exp_2018)
pfit<-prune(fit, cp= 0.01)
rpart.plot(pfit, fallen.leaves=TRUE, tweak=1, gap=0, space=0.5, type=3,clip.right.labs = FALSE,suffix="%", faclen = 0, clip.facs = TRUE, varlen=3,  extra=0, box.palette=c("white"))
```

# Pannel supp data Fig. 4
```{r echo= FALSE, message=FALSE, warning=FALSE, fig.width=7, fig.height= 9}
SDF4A<-ggplot()+
    geom_line(data=est_root_mature_2020_EC,aes(x, predicted*100,colour=group),size=1)+
    geom_ribbon(data=est_root_mature_2020_EC, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.4)+
      scale_fill_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
      scale_color_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
  xlim(210,520)+
  ylim(0,10)+
   theme(axis.text = element_text(size = 12),axis.title = element_text(size = 12),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" ", y=" Proportion of enlarging \n coarse root ring (%)")

SDF4B<-ggplot()+
    geom_line(data=est_root_mature_2018_EC,aes(x, predicted*100),size=1,color="#67a9cf")+
    geom_ribbon(data=est_root_mature_2018_EC, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100),fill="#67a9cf", alpha = 0.4)+
  xlim(210,520)+
  ylim(0,10)+
   theme(axis.text = element_text(size = 12),axis.title = element_text(size = 12),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" ", y=" Proportion enlarging \n coarse root ring (%)")


SDF4C<-ggplot()+
      geom_ribbon(data=est_root_sapling_2017_EC, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.2)+
  geom_line(data=est_root_sapling_2017_EC,aes(x, predicted*100,colour=group),size = 1)+
      scale_fill_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
      scale_color_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
  xlim(210, 490)+
  ylim(0,10)+
   theme(axis.text = element_text(size = 10),axis.title = element_text(size = 10),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" ", y=" Proportion of enlarging \n coarse root ring (%)")

SDF4D<-ggplot()+
    geom_line(data=est_root_sapling_2018_belgium_EC,aes(x, predicted,colour=group),size =1)+
    geom_ribbon(data=est_root_sapling_2018_belgium_EC, aes(x, predicted,ymin = conf.low, ymax = conf.high, fill=group), alpha = 0.4)+    scale_fill_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
      scale_color_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
  xlim(220,335)+
  ylim(0,20)+
   theme(axis.text = element_text(size = 10),axis.title = element_text(size = 10),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" Day of year", y=" Proportion of enlarging \n coarse root ring (%)")

SDF4E<-ggplot()+
    geom_line(data=est_root_sapling_2018_spain_EC,aes(x, predicted*100,colour=group),size =1)+
    geom_ribbon(data=est_root_sapling_2018_spain_EC, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.4)+      scale_fill_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
      scale_color_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
  xlim(220,335)+
  ylim(0,20)+
   theme(axis.text = element_text(size = 10),axis.title = element_text(size = 10),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" Day of year", y=" ")

SDF4F<-ggplot()+
    geom_line(data=est_root_sapling_2018_norway_EC,aes(x, predicted*100,colour=group),size =1)+
    geom_ribbon(data=est_root_sapling_2018_norway_EC, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.4)+      scale_fill_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
      scale_color_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
  xlim(220,335)+
  ylim(0,20)+
   theme(axis.text = element_text(size = 10),axis.title = element_text(size = 10),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" Day of year", y=" ")

SDF4D_bis<-SDF4D + theme(plot.margin = unit(c(0.5, 0.2, 0, 0), "cm"))
SDF4E_bis<- SDF4E+ theme(plot.margin = unit(c(0.5, 0.2, 0, 0), "cm"))
SDF4F_bis <- SDF4F + theme(plot.margin = unit(c(0.5, 0.2, 0, 0), "cm"))
SDF4A_bis<-SDF4A + theme(plot.margin = unit(c(0.5, 0.2, 0, 1), "cm"))
SDF4B_bis<- SDF4B+ theme(plot.margin = unit(c(0.5, 0.2, 0, 1), "cm"))
SDF4C_bis <- SDF4C + theme(plot.margin = unit(c(0.5, 0.2, 0, 1), "cm"))

line3<-plot_grid(SDF4D_bis,SDF4E_bis,SDF4F_bis, nrow=1, labels=c('D', 'E', 'F'))
part1<-plot_grid(SDF4A_bis,SDF4B_bis,SDF4C_bis, nrow=3,rel_heights=c( 1, 1,1), labels=c('A', 'B', 'C'))
SDF4<-plot_grid(part1, line3,nrow=2, rel_heights=c(8, 3));SDF4

```

# Analyse supp data Fig. 5

## A
```{r echo=FALSE, message=FALSE, warning=FALSE}
temp_TH_mature_2018_stem<- subset(temp_TH_mature_2018, organ=="stem")
temp_TH_mature_2018_stem$species<-as.factor(temp_TH_mature_2018_stem$species)
temp_TH_mature_2018_root<- subset(temp_TH_mature_2018, organ=="roots")
temp_TH_mature_2018_root$species<-as.factor(temp_TH_mature_2018_root$species)
```

Output GAM models

__relation coarse roots thickening and soil temperature mature trees 2018__
```{r echo=FALSE, message=FALSE, warning=FALSE}
m_temp_TH_mature_2018_root<-gam(mean_persent_TH/100 ~ s(X10days_temp, k=3),  data=temp_TH_mature_2018_root, family="quasibinomial", method="REML")
summary(m_temp_TH_mature_2018_root)
par(mfrow = c(2, 2))
gam.check(m_temp_TH_mature_2018_root)
est_temp_TH_mature_2018_root<-ggpredict(m_temp_TH_mature_2018_root, terms=c("X10days_temp"))
```

__relation stem thickening and air temperature mature trees 2018__
```{r echo=FALSE, message=FALSE, warning=FALSE}
m_temp_TH_mature_2018_stem<-gam(mean_persent_TH/100 ~ s(X10days_temp, k=3),  data=temp_TH_mature_2018_stem, family="quasibinomial", method="REML")
summary(m_temp_TH_mature_2018_stem)
par(mfrow = c(2, 2))
gam.check(m_temp_TH_mature_2018_stem)
est_temp_TH_mature_2018_stem<-ggpredict(m_temp_TH_mature_2018_stem, terms=c("X10days_temp"))

```

## B
```{r echo=FALSE, message=FALSE, warning=FALSE}
temp_TH_sapling_2018_stem<- subset(temp_TH_sampling_2018_bel, organ=="stem")
temp_TH_sapling_2018_stem$species<-as.factor(temp_TH_sapling_2018_stem$species)
temp_TH_sapling_2018_root<- subset(temp_TH_sampling_2018_bel, organ=="roots")
temp_TH_sapling_2018_root$species<-as.factor(temp_TH_sapling_2018_root$species)
```

__relation coarse roots thickening and soil temperature saplings Belgium 2018__
```{r echo=FALSE, message=FALSE, warning=FALSE}
m_temp_TH_sapling_2018_root<-gam(mean_persent_TH/100 ~ s(X10days_temp,by = species, k=3),  data=temp_TH_sapling_2018_root, family="quasibinomial", method="REML")
summary(m_temp_TH_sapling_2018_root)
par(mfrow = c(2, 2))
gam.check(m_temp_TH_sapling_2018_root)
est_temp_TH_sapling_2018_root<-ggpredict(m_temp_TH_sapling_2018_root, terms=c("X10days_temp","species"))
```

__relation stem thickening and air temperature saplings Belgium 2018__
```{r echo=FALSE, message=FALSE, warning=FALSE}
m_temp_TH_sapling_2018_stem<-gam(mean_persent_TH/100 ~ s(X10days_temp,by = species, k=3), data=temp_TH_sapling_2018_stem, family="quasibinomial", method="REML")
summary(m_temp_TH_sapling_2018_stem)
par(mfrow = c(2, 2))
gam.check(m_temp_TH_sapling_2018_stem)
est_temp_TH_sapling_2018_stem<-ggpredict(m_temp_TH_sapling_2018_stem, terms=c("X10days_temp","species"))
```

# Pannel supp data Fig.5

```{r echo=FALSE, message=FALSE, warning=FALSE, results='hide',fig.width= 13, fig.height= 5}

SDF5A<-ggplot()+
    geom_line(data=est_temp_TH_mature_2018_root,aes(x, predicted*100,colour=group),size=1,linetype="dashed")+
    geom_ribbon(data=est_temp_TH_mature_2018_root, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.4)+
  geom_line(data=est_temp_TH_mature_2018_stem,aes(x, predicted*100,colour=group),size=1)+
    geom_ribbon(data=est_temp_TH_mature_2018_stem, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.4)+
      scale_fill_manual(values=c("#67a9cf"))+
      scale_color_manual(values=c("#67a9cf"))+
  xlim(2,20)+
  ylim(0,75)+
   theme(axis.text = element_text(size = 14),axis.title = element_text(size = 14),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" Mean 10 days temperature prior to sampling ", y="Proportion of \n thickening ring (%)")

SDF5B<-ggplot()+
  geom_line(data=est_temp_TH_sapling_2018_root,aes(x, predicted*100,colour=group),size=1, linetype="dashed")+
    geom_ribbon(data=est_temp_TH_sapling_2018_root, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.4)+
  geom_line(data=est_temp_TH_sapling_2018_stem,aes(x, predicted*100,colour=group),size=1)+
    geom_ribbon(data=est_temp_TH_sapling_2018_stem, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.4)+
      scale_fill_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
      scale_color_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
  xlim(2,20)+
  ylim(0,75)+
   theme(axis.text = element_text(size = 14),axis.title = element_text(size = 14),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x="Mean 10 days temperature prior to sampling", y="")

SDF5A_bis<-SDF5A + theme(plot.margin = unit(c(0, 1, 0, 0), "cm"))
SDF5B_bis<-SDF5B + theme(plot.margin = unit(c(0, 1, 0, 0), "cm"))

SDF5<-plot_grid(SDF5A_bis,SDF5B_bis, nrow=1, labels=c('A', 'B')); SDF5 
```



# Pannel supp data Fig. 6

```{r echo=FALSE, message=FALSE, warning=FALSE, results='hide',fig.width= 10, fig.height= 10}
SDF6A<-ggplot()+
    geom_ribbon(data=est_temp_TH_mature_2020_root, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.4)+
   geom_line(data=est_temp_TH_mature_2020_root,aes(x, predicted*100,colour=group), size=1,linetype="dashed")+
      geom_point(data= temp_TH_mature_2020_root, aes(X10days_temp,mean_persent_TH, colour=species))+
      scale_fill_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
      scale_color_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
  xlim(2,20)+
  ylim(0,50)+
   theme(axis.text = element_text(size = 12),axis.title = element_text(size = 12),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" Mean 10 days temperature prior sampling (°C) ", y="Proportion of \n thickening ring (%)")

SDF6B<-ggplot()+
   geom_line(data=est_temp_TH_mature_2018_root,aes(x, predicted*100),size=1, linetype="dashed",colour="#67a9cf")+
    geom_ribbon(data=est_temp_TH_mature_2018_root, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100), alpha = 0.4,fill="#67a9cf")+
    geom_point(data=temp_TH_mature_2018_root, aes(X10days_temp,mean_persent_TH), colour="#67a9cf")+
      scale_fill_manual(values=c("#67a9cf"))+
      scale_color_manual(values=c("#67a9cf"))+
  xlim(2,20)+
  ylim(0,50)+
   theme(axis.text = element_text(size = 12),axis.title = element_text(size = 12),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" Mean 10 days temperature prior to sampling ", y="Proportion of \n thickening ring (%)")

SDF6C<-ggplot()+
    geom_ribbon(data=est_temp_TH_sampling_2017_root, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.4)+
  geom_point(data=temp_TH_sampling_2017_root , aes(X10days_temp,mean_persent_TH, colour=species))+
   geom_line(data=est_temp_TH_sampling_2017_root,aes(x, predicted*100,colour=group),size=1)+
      scale_fill_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
      scale_color_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
  xlim(2,20)+
  ylim(0,50)+
   theme(axis.text = element_text(size = 12),axis.title = element_text(size = 12),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" Mean 10 days temperature prior sampling (°C) ", y="Proportion of \n thickening ring (%)")

SDF6D<-ggplot()+
  geom_line(data=est_temp_TH_sapling_2018_root,aes(x, predicted*100,colour=group),size=1)+
    geom_ribbon(data=est_temp_TH_sapling_2018_root, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.4)+
  geom_point(data=temp_TH_sapling_2018_root , aes(X10days_temp,mean_persent_TH, colour=species))+
      scale_fill_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
      scale_color_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
  xlim(2,20)+
  ylim(0,50)+
   theme(axis.text = element_text(size = 12),axis.title = element_text(size = 12),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" Mean 10 days temperature prior sampling (°C) ", y="Proportion of \n thickening ring (%)")


SDF6<-plot_grid(SDF6A,SDF6B,SDF6C,SDF6D, nrow=2, labels=c('A', 'B',"C", 'D'));SDF6 
```


# Pannel supp data Fig. 7
```{r echo=FALSE, message=FALSE, warning=FALSE, results='hide',fig.width= 8, fig.height= 10}
SDF7A<-ggplot()+
    geom_ribbon(data=est_root_sapling_2018_belgium, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.4)+
   geom_line(data=est_root_sapling_2018_belgium,aes(x, predicted*100,colour=group),size =1)+
    geom_point(data=root_sapling_2018_belgium, aes(DOY_cumul, percent_TH,colour=species))+
      scale_fill_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
      scale_color_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
  xlim(220,335)+
   ylim(0,80)+
   theme(axis.text = element_text(size = 10),axis.title = element_text(size = 10),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" Day of year", y=" Proportion of thickening \n coarse root ring (%)")

SDF7B<-ggplot()+
    geom_ribbon(data=est_root_sapling_2018_spain, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.4)+
      geom_line(data=est_root_sapling_2018_spain,aes(x, predicted*100,colour=group),size =1)+
    geom_point(data=root_sapling_2018_spain, aes(DOY_cumul,percent_TH, colour=species))+
      scale_fill_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
      scale_color_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
  xlim(220,335)+
   ylim(0,80)+
   theme(axis.text = element_text(size = 10),axis.title = element_text(size = 10),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" Day of year ", y= "Proportion of thickening \n  coarse root ring (%)")

SDF7C<-ggplot()+
    geom_ribbon(data=est_root_sapling_2018_norway, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.4)+
      geom_line(data=est_root_sapling_2018_norway,aes(x, predicted*100,colour=group),size =1)+
    geom_point(data= root_sapling_2018_norway, aes(DOY_cumul,percent_TH, colour=species))+
      scale_fill_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
      scale_color_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
  xlim(220,335)+
  ylim(0,80)+
   theme(axis.text = element_text(size = 10),axis.title = element_text(size = 10),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" Day of year ", y="Proportion of thickening \n  coarse root ring (%)")

SDF7D<-ggplot()+
    geom_line(data=est_root_sapling_2018_belgium_EC,aes(x, predicted,colour=group),size =1)+
    geom_ribbon(data=est_root_sapling_2018_belgium_EC, aes(x, predicted,ymin = conf.low, ymax = conf.high, fill=group), alpha = 0.4)+ 
    geom_point(data=root_sapling_2018_belgium, aes(DOY_cumul,EC_poisson,colour=species))+
  scale_fill_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
      scale_color_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
  xlim(220,335)+
  ylim(0,50)+
   theme(axis.text = element_text(size = 10),axis.title = element_text(size = 10),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" Day of year", y=" Proportion of enlarging \n coarse root ring (%)")

SDF7E<-ggplot()+
    geom_line(data=est_root_sapling_2018_spain_EC,aes(x, predicted*100,colour=group),size =1)+
    geom_ribbon(data=est_root_sapling_2018_spain_EC, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.4)+
    geom_point(data=root_sapling_2018_spain, aes(DOY_cumul,percent_EC,colour=species))+
  scale_fill_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
      scale_color_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
  xlim(220,335)+
  ylim(0,50)+
   theme(axis.text = element_text(size = 10),axis.title = element_text(size = 10),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" Day of year", y=" Proportion of enlarging \n coarse root ring (%)")

SDF7F<-ggplot()+
    geom_line(data=est_root_sapling_2018_norway_EC,aes(x, predicted*100,colour=group),size =1)+
    geom_ribbon(data=est_root_sapling_2018_norway_EC, aes(x, predicted*100,ymin = conf.low*100, ymax = conf.high*100, fill=group), alpha = 0.4)+
    geom_point(data=root_sapling_2018_spain, aes(DOY_cumul,percent_EC, colour=species))+
  scale_fill_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
      scale_color_manual(values=c("#fc8d59", "#67a9cf","#d7191c", "#5e3c99"))+
  xlim(220,335)+
    ylim(0,50)+
   theme(axis.text = element_text(size = 10),axis.title = element_text(size = 10),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" Day of year", y=" Proportion of enlarging \n coarse root ring (%)")


SDF7<-plot_grid(SDF7A,SDF7B,SDF7C,SDF7D,SDF7E,SDF7F, nrow=3, ncol=2, labels=c('A', 'B', 'C','D', 'E', 'F'));SDF7

```


# Pannel supp data Fig. 9
```{r echo=FALSE, message=FALSE, warning=FALSE, results='hide',fig.width= 8, fig.height= 3}

climate_sapling_2017$date<-as.Date(climate_sapling_2017$date)
break.vec <- c(as.Date("2017-09-07"),as.Date("2017-10-27"),as.Date("2017-12-16"),as.Date("2018-02-05"),as.Date("2018-03-26"), as.Date("2018-05-15"))

SDF9A<-ggplot()+
      geom_ribbon( aes(x = date, ymin = min_TS, ymax = max_TS), data=climate_sapling_2017 ,fill = "#80cdc1")+
  geom_line( aes(x= date, y = mean_TS), data=climate_sapling_2017,size=1, color='#018571' )+
   geom_hline(yintercept=5, linetype="dashed", color="grey")+
  ylim(-10,25)+
    scale_x_date( breaks = break.vec , date_labels = "%j", limits=c(as.Date("2017-09-10"), as.Date("2018-05-05")))+
 theme(axis.text = element_text(size = 12),axis.title = element_text(12),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" Day of year", y="Within pot \n temperature (ºC) ")

break.vec <- c(as.Date("2018-08-28"),as.Date("2018-09-27"),as.Date("2018-10-27"),as.Date("2018-11-26"))

SDF9B<-ggplot()+
   geom_ribbon( aes(date, ymin = min_TS, ymax = max_TS), data=climate_sapling_2018_belgium ,fill = "#80cdc1")+
  geom_line( aes(x= date, y = mean_TS), data=climate_sapling_2018_belgium,size=1, color='#018571' )+
    scale_x_date( breaks = break.vec , date_labels = "%j", limits=c(as.Date("2018-08-08"), as.Date("2018-12-01")))+
  ylim(-10,25)+
     geom_hline(yintercept=5, linetype="dashed", color="grey")+
 theme(axis.text = element_text(size = 12),axis.title = element_text(size = 12),legend.position="none",axis.text.x=element_text(angle=0, hjust=1),panel.background = element_rect(fill = "white",colour="grey"),panel.grid.major = element_line(colour = "white"),legend.key = element_rect(colour = NA, fill = NA))+ labs(x=" Day of year ", y=" ")

SDF9<-plot_grid(SDF9A, SDF9B, nrow=1, labels=c('A','B'),rel_widths = c(2.5, 1.5)); SDF9

```

