These codes produce the results of figure 5 and figure S4. The files used are located in the folder named “Correlations” and contain the data (charcoal measurements, fire histories, pollen percentages, PAR, and mean July temperature) based on the outputs created to produce the figures and results for each site (Rosalia_corr.txt, Rosalia_corr_charcoal.txt, and Pikku_corr.txt, Pikku_charcoal.txt).

#install.packages("ggplot2") #if not already installed
library("ggplot2")
## Warning: package 'ggplot2' was built under R version 4.0.5
windowsFonts(Times=windowsFont("Times New Roman"))
#install.packages("ggridges") #if not already installed
library(ggridges)
#install.packages("reshape2") #if not already installed
library(reshape2)
#install.packages("viridis") #if not already installed
library(viridis)
## Warning: package 'viridis' was built under R version 4.0.5
## Loading required package: viridisLite
## Warning: package 'viridisLite' was built under R version 4.0.5
#install.packages("gridExtra") #if not already installed
library(gridExtra)
dirs<-c('Correlations/')
Files<-list.files(dirs)

Pikku_charcoal<-read.table("Correlations/Pikku_corr_charcoal.txt",head=T,sep=" ")
Rosalia_charcoal<-read.table("Correlations/Rosalia_corr_charcoal.txt",head=T,sep=" ")
Pikku_Fire<-read.table("Correlations/Pikku_corr.txt",head=T,sep=" ")
Rosalia_Fire<-read.table("Correlations/Rosalia_corr.txt",head=T,sep=" ")

##Rosalia

BB<-data.frame()
FF<-data.frame()
FS<-data.frame()
TJ<-data.frame()
for (u in 1:999){  #number of iterations
COL_BB<-c(6:18)
COL_FF<-c(6:18)
COL_FS<-c(6:18)
COL_TJ<-c(6:18)
x<-Rosalia_charcoal[sample(nrow(Rosalia_charcoal), 110), ]
for (i in 1:13){
BB_x<-cor.test(x[,COL_BB[i]],x[,2], method = c("pearson"))$estimate[[1]]
FF_x<-cor.test(x[,COL_FF[i]],x[,3], method = c("pearson"))$estimate[[1]]
FS_x<-cor.test(x[,COL_FS[i]],x[,4], method = c("pearson"))$estimate[[1]]
TJ_x<-cor.test(x[,COL_TJ[i]],x[,5], method = c("pearson"))$estimate[[1]]

BB[u,i]<-BB_x
FF[u,i]<-FF_x
FS[u,i]<-FS_x
TJ[u,i]<-TJ_x
}
}

colnames(BB)<-c("Pinus","Picea","Betula","Alnus","Juniperus","Salix","Ericaceae","Poaceae","Artemisia","Cyperaceae","Sphagnum","PAR","ROC")
colnames(FF)<-c("Pinus","Picea","Betula","Alnus","Juniperus","Salix","Ericaceae","Poaceae","Artemisia","Cyperaceae","Sphagnum","PAR","ROC")
colnames(FS)<-c("Pinus","Picea","Betula","Alnus","Juniperus","Salix","Ericaceae","Poaceae","Artemisia","Cyperaceae","Sphagnum","PAR","ROC")  
colnames(TJ)<-c("Pinus","Picea","Betula","Alnus","Juniperus","Salix","Ericaceae","Poaceae","Artemisia","Cyperaceae","Sphagnum","PAR","ROC")
BB[is.na(BB)] <- 0
FF[is.na(FF)] <- 0
FS[is.na(FS)] <- 0
TJ[is.na(TJ)] <- 0

LW<-data.frame()
for (u in 1:999){  #number of iterations
COL_LW<-c(2:5,8:20)
x<-Rosalia_Fire[sample(nrow(Rosalia_Fire), 84), ]
for (i in 1:17){
LW_x<-cor.test(x[,COL_LW[i]],x[,6], method = c("pearson"))$estimate[[1]]
LW[u,i]<-LW_x
}
}
  
colnames(LW)<-c("Biomass burned","Fire frequency","Fire size","July temperature","Pinus","Picea","Betula","Alnus","Juniperus","Salix","Ericaceae","Poaceae","Artemisia","Cyperaceae","Sphagnum","PAR","ROC")

LW[is.na(LW)] <- 0

BB_plot<-t(BB)
BB_plot<-data.frame(rownames(BB_plot),BB_plot)
BB_plot<-melt(BB_plot,by="rownames.BB_plot.")
## Using rownames.BB_plot. as id variables
colnames(BB_plot)[1]<-'Species'
BB_plot$Species<-factor(BB_plot$Species,levels=rev(colnames(BB)))

FF_plot<-t(FF)
FF_plot<-data.frame(rownames(FF_plot),FF_plot)
FF_plot<-melt(FF_plot,by="rownames.FF_plot.")
## Using rownames.FF_plot. as id variables
colnames(FF_plot)[1]<-'Species'
FF_plot$Species<-factor(FF_plot$Species,levels=rev(colnames(FF)))

FS_plot<-t(FS)
FS_plot<-data.frame(rownames(FS_plot),FS_plot)
FS_plot<-melt(FS_plot,by="rownames.FS_plot.")
## Using rownames.FS_plot. as id variables
colnames(FS_plot)[1]<-'Species'
FS_plot$Species<-factor(FS_plot$Species,levels=rev(colnames(FS)))

LW_plot<-t(LW)
LW_plot<-data.frame(rownames(LW_plot),LW_plot)
LW_plot<-melt(LW_plot,by="rownames.LW_plot.")
## Using rownames.LW_plot. as id variables
colnames(LW_plot)[1]<-'Species'
LW_plot$Species<-factor(LW_plot$Species,levels=rev(colnames(LW)))

TJ_plot<-t(TJ)
TJ_plot<-data.frame(rownames(TJ_plot),TJ_plot)
TJ_plot<-melt(TJ_plot,by="rownames.TJ_plot.")
## Using rownames.TJ_plot. as id variables
colnames(TJ_plot)[1]<-'Species'
TJ_plot$Species<-factor(TJ_plot$Species,levels=rev(colnames(TJ)))

RosBB_plot<-BB_plot
RosFF_plot<-FF_plot
RosFS_plot<-FS_plot
RosLW_plot<-LW_plot
RosTJ_plot<-TJ_plot
#Pikku         
BB<-data.frame()
FF<-data.frame()
FS<-data.frame()
TJ<-data.frame()
for (u in 1:999){  #number of iterations
COL_BB<-c(6:18)
COL_FF<-c(6:18)
COL_FS<-c(6:18)
COL_TJ<-c(6:18)
x<-Pikku_charcoal[sample(nrow(Pikku_charcoal), 108), ]
for (i in 1:13){
BB_x<-cor.test(x[,COL_BB[i]],x[,2], method = c("pearson"))$estimate[[1]]
FF_x<-cor.test(x[,COL_FF[i]],x[,3], method = c("pearson"))$estimate[[1]]
FS_x<-cor.test(x[,COL_FS[i]],x[,4], method = c("pearson"))$estimate[[1]]
TJ_x<-cor.test(x[,COL_TJ[i]],x[,5], method = c("pearson"))$estimate[[1]]

BB[u,i]<-BB_x
FF[u,i]<-FF_x
FS[u,i]<-FS_x
TJ[u,i]<-TJ_x
}
}

colnames(BB)<-c("Pinus","Picea","Betula","Alnus","Juniperus","Salix","Ericaceae","Poaceae","Artemisia","Cyperaceae","Sphagnum","PAR","ROC")
colnames(FF)<-c("Pinus","Picea","Betula","Alnus","Juniperus","Salix","Ericaceae","Poaceae","Artemisia","Cyperaceae","Sphagnum","PAR","ROC")
colnames(FS)<-c("Pinus","Picea","Betula","Alnus","Juniperus","Salix","Ericaceae","Poaceae","Artemisia","Cyperaceae","Sphagnum","PAR","ROC")  
colnames(TJ)<-c("Pinus","Picea","Betula","Alnus","Juniperus","Salix","Ericaceae","Poaceae","Artemisia","Cyperaceae","Sphagnum","PAR","ROC")
BB[is.na(BB)] <- 0
FF[is.na(FF)] <- 0
FS[is.na(FS)] <- 0
TJ[is.na(TJ)] <- 0

LW<-data.frame()
for (u in 1:999){  #number of iterations
COL_LW<-c(2:5,8:20)
x<-Pikku_Fire[sample(nrow(Pikku_Fire), 92), ]
for (i in 1:17){
LW_x<-cor.test(x[,COL_LW[i]],x[,6], method = c("pearson"))$estimate[[1]]
LW[u,i]<-LW_x
}
}
  
colnames(LW)<-c("Biomass burned","Fire frequency","Fire size","July temperature","Pinus","Picea","Betula","Alnus","Juniperus","Salix","Ericaceae","Poaceae","Artemisia","Cyperaceae","Sphagnum","PAR","ROC")

LW[is.na(LW)] <- 0

BB_plot<-t(BB)
BB_plot<-data.frame(rownames(BB_plot),BB_plot)
BB_plot<-melt(BB_plot,by="rownames.BB_plot.")
## Using rownames.BB_plot. as id variables
colnames(BB_plot)[1]<-'Species'
BB_plot$Species<-factor(BB_plot$Species,levels=rev(colnames(BB)))

FF_plot<-t(FF)
FF_plot<-data.frame(rownames(FF_plot),FF_plot)
FF_plot<-melt(FF_plot,by="rownames.FF_plot.")
## Using rownames.FF_plot. as id variables
colnames(FF_plot)[1]<-'Species'
FF_plot$Species<-factor(FF_plot$Species,levels=rev(colnames(FF)))

FS_plot<-t(FS)
FS_plot<-data.frame(rownames(FS_plot),FS_plot)
FS_plot<-melt(FS_plot,by="rownames.FS_plot.")
## Using rownames.FS_plot. as id variables
colnames(FS_plot)[1]<-'Species'
FS_plot$Species<-factor(FS_plot$Species,levels=rev(colnames(FS)))

LW_plot<-t(LW)
LW_plot<-data.frame(rownames(LW_plot),LW_plot)
LW_plot<-melt(LW_plot,by="rownames.LW_plot.")
## Using rownames.LW_plot. as id variables
colnames(LW_plot)[1]<-'Species'
LW_plot$Species<-factor(LW_plot$Species,levels=rev(colnames(LW)))

TJ_plot<-t(TJ)
TJ_plot<-data.frame(rownames(TJ_plot),TJ_plot)
TJ_plot<-melt(TJ_plot,by="rownames.TJ_plot.")
## Using rownames.TJ_plot. as id variables
colnames(TJ_plot)[1]<-'Species'
TJ_plot$Species<-factor(TJ_plot$Species,levels=rev(colnames(TJ)))

PikBB_plot<-BB_plot
PikFF_plot<-FF_plot
PikFS_plot<-FS_plot
PikLW_plot<-LW_plot
PikTJ_plot<-TJ_plot

##Graphs

RosBB_plot<-cbind(rep(1,nrow(RosBB_plot)),RosBB_plot)
colnames(RosBB_plot)[1]<-"Site"
RosFF_plot<-cbind(rep(1,nrow(RosFF_plot)),RosFF_plot)
colnames(RosFF_plot)[1]<-"Site"
RosFS_plot<-cbind(rep(1,nrow(RosFS_plot)),RosFS_plot)
colnames(RosFS_plot)[1]<-"Site"
RosLW_plot<-cbind(rep(1,nrow(RosLW_plot)),RosLW_plot)
colnames(RosLW_plot)[1]<-"Site"
RosTJ_plot<-cbind(rep(1,nrow(RosTJ_plot)),RosTJ_plot)
colnames(RosTJ_plot)[1]<-"Site"

PikBB_plot<-cbind(rep(2,nrow(PikBB_plot)),PikBB_plot)
colnames(PikBB_plot)[1]<-"Site"
PikFF_plot<-cbind(rep(2,nrow(PikFF_plot)),PikFF_plot)
colnames(PikFF_plot)[1]<-"Site"
PikFS_plot<-cbind(rep(2,nrow(PikFS_plot)),PikFS_plot)
colnames(PikFS_plot)[1]<-"Site"
PikLW_plot<-cbind(rep(2,nrow(PikLW_plot)),PikLW_plot)
colnames(PikLW_plot)[1]<-"Site"
PikTJ_plot<-cbind(rep(2,nrow(PikTJ_plot)),PikTJ_plot)
colnames(PikTJ_plot)[1]<-"Site"

ComBB<-rbind(RosBB_plot,PikBB_plot)
ComBB$Site<-factor(ComBB$Site,labels=c("Rosalia","Pikku"))
ComFF<-rbind(RosFF_plot,PikFF_plot)
ComFF$Site<-factor(ComFF$Site,labels=c("Rosalia","Pikku"))
ComFS<-rbind(RosFS_plot,PikFS_plot)
ComFS$Site<-factor(ComFS$Site,labels=c("Rosalia","Pikku"))
ComLW<-rbind(RosLW_plot,PikLW_plot)
ComLW$Site<-factor(ComLW$Site,labels=c("Rosalia","Pikku"))
ComTJ<-rbind(RosTJ_plot,PikTJ_plot)
ComTJ$Site<-factor(ComTJ$Site,labels=c("Rosalia","Pikku"))

cols <-  alpha(c("Rosalia" = "gold3", "Pikku" = "lightseagreen"),0.4)

BB_p<-    
ggplot(ComBB,aes(x=value,y=Species,fill=Site))+
   theme_minimal()+
   geom_density_ridges_gradient(scale = 1, size = 0.3, rel_min_height = 0.01)+
   scale_fill_manual(values=cols)+
   scale_x_continuous(limits=c(-1,1))+
   geom_vline(xintercept = 0,linetype="dotted")+
   labs(title="Biomass burned")+
   theme(text=element_text(family="Times"),
   axis.title.y=element_blank(),
   axis.title.x=element_blank(),
   axis.text.y=element_text(face="italic"),
   axis.text.x=element_text())
BB_p
## Picking joint bandwidth of 0.0161

FF_p<-    
ggplot(ComFF,aes(x=value,y=Species,fill=Site))+
   theme_minimal()+
   geom_density_ridges_gradient(scale = 1, size = 0.3, rel_min_height = 0.01)+
   scale_fill_manual(values=cols)+
   scale_x_continuous(limits=c(-1,1))+
   geom_vline(xintercept = 0,linetype="dotted")+
   labs(title="Fire frequency")+
   theme(text=element_text(family="Times"),
   axis.title.y=element_blank(),
   axis.title.x=element_blank(),
   axis.text.y=element_text(face="italic"),
   axis.text.x=element_text())
FF_p
## Picking joint bandwidth of 0.0143

FS_p<-    
ggplot(ComFS,aes(x=value,y=Species,fill=Site))+
   theme_minimal()+
   geom_density_ridges_gradient(scale = 1, size = 0.3, rel_min_height = 0.01)+
   scale_fill_manual(values=cols)+
   scale_x_continuous(limits=c(-1,1))+
   geom_vline(xintercept = 0,linetype="dotted")+
   labs(title="Fire size")+
   theme(text=element_text(family="Times"),
   axis.title.y=element_blank(),
   axis.title.x=element_blank(),
   axis.text.y=element_text(face="italic"),
   axis.text.x=element_text())
FS_p
## Picking joint bandwidth of 0.0124

TJ_p<-    
ggplot(ComTJ,aes(x=value,y=Species,fill=Site))+
   theme_minimal()+
   geom_density_ridges_gradient(scale = 1, size = 0.3, rel_min_height = 0.01)+
   scale_fill_manual(values=cols)+
   scale_x_continuous(limits=c(-1,1))+
   geom_vline(xintercept = 0,linetype="dotted")+
   labs(title="July temperature")+
   theme(text=element_text(family="Times"),
   axis.title.y=element_blank(),
   axis.title.x=element_blank(),
   axis.text.y=element_text(face="italic"),
   axis.text.x=element_text())
TJ_p
## Picking joint bandwidth of 0.0149

LW_p<-    
ggplot(ComLW,aes(x=value,y=Species,fill=Site))+
   theme_minimal()+
   geom_density_ridges_gradient(scale = 1, size = 0.3, rel_min_height = 0.01)+
   scale_fill_manual(values=cols)+
   scale_x_continuous(limits=c(-1,1))+
   geom_vline(xintercept = 0,linetype="dotted")+
   labs(title="Portion of charcoal with wood aspect")+
   theme(text=element_text(family="Times"),
   axis.title.y=element_blank(),
   axis.title.x=element_blank(),
   axis.text.y=element_text(face="italic"),
   axis.text.x=element_text())
LW_p
## Picking joint bandwidth of 0.0169