Here briefly summarized the process to figure visualization based on ggplot2 packages for the manuscript (DOI coming soon). Note the descriptions may be updated at any time, and there is concern that they may contain incomplete and redundant code.The data used are not published online, but may be shared by contacting the corresponding author.

Figure 2

(a)Scatter plots by PC1 and PC2. The blue dots indicate the western provenance regions and the red dots indicate the eastern provenance regions. (b)Permutations of PC1 and PC2 loadings consisting of monthly meteorological factors. In this analysis, April and May were marked as spring, June through August as summer, September and October as fall, and November through March as winter. Abbreviations for meteorological factors was summarized in Supplemental Table 1.

#package preparation
library(ggplot2)
library(ggfortify)
library(ggrepel)
library(patchwork)
library(reshape2)

#Data adjustment
m<-read.csv("/Users/sugai/Desktop/tgg_meteodata.csv")
#dim(m) #200 70
mp<-prcomp(m[,7:70],scale=TRUE)
#autoplot(mp) #PC1:48.8,PC2:27.6
mpx<-as.data.frame(mp$x)[,1:2]
mpx$region<-m$region
mpx$row<-m$row
mpx$col<-m$col
mpx<-transform(mpx,region=factor(region,levels=c("OTI","URK","ASR","SRM","HST","ASB","KTB","KCA")))

mpl<-as.data.frame(mp$rotation)[,1:2]
mpl$factor<-as.factor(rownames(mpl))
write.csv(mpl,"/Users/sugai/Desktop/tgg_meteodata_loading.csv")
#Season factors were manually added (11~3:Winter,4,5:Spring,6~8:Summer,9,10:Autumn)
mpl2<-read.csv("/Users/sugai/Desktop/tgg_meteodata_loading2.csv")
mpl2<-transform(mpl2,season=factor(season,levels=c("Spring","Summer","Autumn","Winter")))

t<-read.csv("/Users/sugai/Desktop/tgg_toporogy.csv")
t2<-merge(mpx,t,by=c("region","row","col"))
M<-cor(t2[,c(4:5,7:9)])
M2<-t(M[1:2,3:5])
M3<-melt(M2)
M3<-transform(M3,Var1=factor(Var1,levels=c("lat","lon","aveAl")))

#Visualization
a<-
ggplot(data=mpx,
       mapping=aes(y=PC1,x=PC2,colour=region,fill=region,label=region))+
    stat_ellipse(geom = "polygon",alpha=.7)+
  geom_point(data=mpx,
             mapping=aes(y=PC1,x=PC2,fill=region),
             shape=21,colour="black",size=1,alpha=.6)+
  scale_fill_brewer(palette = "RdBu")+
  scale_color_brewer(palette = "RdBu")+
  theme_classic(base_family = "Helvetica")+
  labs(y="PC1 (48.8%)",x="PC2 (27.6%)")+
  theme(legend.position="bottom",
        legend.direction = "horizontal",
        legend.text = element_text(size=rel(1),colour="black"),
        legend.title = element_text(size=rel(1),colour="black"),
        plot.title = element_text(size=rel(2),colour="black"),
        axis.text = element_text(size=rel(1),colour="black"),
        axis.title = element_text(size=rel(2),colour="black"))+
  guides(fill = guide_legend(title="Provenance region",nrow = 1, byrow = TRUE))+
  guides(colour = guide_legend(title="Provenance region",nrow = 1, byrow = TRUE))

b<-
ggplot(data=mpl2,aes(x=reorder(x=factor,X=PC1,fun=mean),y=PC1,fill=season))+
  geom_col(alpha=.75,color="black")+
    labs(y="PC1 loading",x="")+
theme_classic(base_family = "Helvetica")+
    scale_fill_manual(values=c("#ffcccc","#ff0066","#99ccff","#0066cc"))+
    theme(legend.position="top",legend.justification = "center",
          legend.direction = "horizontal",
          legend.text = element_text(size=rel(.7),colour="black"),
          legend.title = element_text(size=rel(.7),colour="black"),
          axis.text.y = element_text(size=rel(.5),colour="black"),
          axis.text.x = element_text(size=rel(1),colour="black"),
          axis.title = element_text(size=rel(2),colour="black"))+
    coord_flip()

c<-
ggplot(data=mpl2,aes(x=reorder(x=factor,X=PC2,fun=mean),y=PC2,fill=season))+
  geom_col(alpha=.75,color="black")+
    labs(y="PC2 loading",x="")+
theme_classic(base_family = "Helvetica")+
    scale_fill_manual(values=c("#ffcccc","#ff0066","#99ccff","#0066cc"))+
    theme(legend.position="none",
          plot.background = element_blank(),
          axis.text.y = element_text(size=rel(.5),colour="black"),
          axis.text.x = element_text(size=rel(1),colour="black"),
          axis.title = element_text(size=rel(2),colour="black"))+
    coord_flip()

d<-ggplot(data=M3,aes(x=Var1,y=value,fill=Var1))+
    geom_col(alpha=.7)+
    theme_classic()+
    labs(y="Correlation coefficient",x="")+
    scale_x_discrete(labels=c("lat"="Latitude","lon"="Longitude","aveAl"="Altitude"))+
    scale_fill_manual(values=c("#ff3399","#c2059c","#2952c2"))+
    theme(legend.position="none",
          plot.title = element_text(size=rel(2),colour="black"),
          plot.background = element_blank(),
          axis.text = element_text(size=rel(1),colour="black"),
          axis.title = element_text(size=rel(1.3),colour="black"))+
    facet_wrap(.~M3$Var2,scales="free_y")+
    coord_flip()

(a|{(d/(b|c)) + plot_layout(heights = c(.5, 3))}) + plot_annotation(tag_levels = 'a', tag_prefix = '(',
                                                                   tag_sep = '.', tag_suffix = ')')

Figure 3

Scatter plots of PCA in the selected 8 candidate traits for snow damage resistance. The blue dots indicate the western provenance regions and the red dots indicate the eastern provenance regions. Green arrows with annotated text denote the PC loading of traits as abbreviations shown in Table 1.

#package preparation
library(ggplot2)
library(ggfortify)
library(ggrepel)
library(patchwork)
library(reshape2)

#Data adjustment
st<-read.csv("/Users/sugai/Desktop/tgg_trait_candi.csv")
#dim(st)
#72 11

sp<-prcomp(st[,4:11],scale=TRUE)
#autoplot(sp)
#PC1:56.9,PC2:14.0

spx<-as.data.frame(sp$x)[,1:2]
spx$region<-st$region
spx<-transform(spx,region=factor(region,levels=c("OTI","URK","ASR","SRM","HST","ASB","KTB","KCA")))
spl<-as.data.frame(sp$rotation)[,1:2]
spl$label<-as.factor(rownames(spl))
spl$angle<-c(-35,-35,-15,-35,30,30,30,-15)

#Visualization
ggplot()+
  geom_point(data=spx,
             mapping=aes(x=PC1,y=PC2,fill=region),
             shape=21,colour="black",size=4,alpha=.6)+
  scale_fill_brewer(palette = "RdBu")+
  theme_classic(base_family = "Helvetica")+
  labs(x="PC1 (56.9%)",y="PC2 (14.0%)")+
  geom_segment(data=spl,aes(x=0,y=0,xend=PC1*6,yend=PC2*6),
               colour="#008700",arrow=arrow(length = unit(0.05, "inches"),type = "closed"),
               linetype=1,size=.3,alpha=.7)+
  geom_text(data=spl,aes(x=(PC1*6.8),y=(PC2*6.5), label=label,angle=angle),
                   size=4,colour="black",fontface="bold",show.legend = FALSE)+
  theme(legend.position="bottom",
        legend.direction = "horizontal",
        legend.text = element_text(size=rel(1),colour="black"),
        legend.title = element_text(size=rel(1),colour="black"),
        plot.title = element_text(size=rel(2),colour="black"),
        axis.text = element_text(size=rel(1),colour="black"),
        axis.title = element_text(size=rel(2),colour="black"))+
  guides(fill = guide_legend(title="Provenance region",nrow = 1, byrow = TRUE))

Figure 4

Correlations between PC1 of the candidate traits for snow damage resistance and mean daylight hour from January to March in provenance environments.The blue dots indicate the western provenance regions and the red dots indicate the eastern provenance regions.Black dotted line denote the linear regression.

#package preparation
library(ggplot2)
library(ggfortify)
library(ggrepel)
library(patchwork)
library(reshape2)
library(ggpmisc)

#function define
data_summary <- function(data, varname, groupnames){
  require(plyr)
  summary_func <- function(x, col){
    c(mean = mean(x[[col]], na.rm=TRUE),
      sd = sd(x[[col]], na.rm=TRUE))
  }
  data_sum<-ddply(data, groupnames, .fun=summary_func,
                  varname)
  data_sum <- rename(data_sum, c("mean" = varname))
 return(data_sum)
}

#Data adjustment
st<-read.csv("/Users/sugai/Desktop/tgg_trait_candi.csv")
#dim(st)
#72 11
sp<-prcomp(st[,4:11],scale=TRUE)
spx<-as.data.frame(sp$x)[,1:2]
spx$region<-st$region


m<-read.csv("/Users/sugai/Desktop/tgg_meteodata.csv")
#dim(m)
#200 70
res<-data_summary(data=m,varname=colnames(m)[2],groupnames = "region")
for(i in 7:70){
rest<-data_summary(data=m,varname=colnames(m)[i],groupnames = "region")
rest[,3]<-NULL
res<-merge(res,rest,by="region")
}

spm<-merge(spx,res,by="region")
spm<-transform(spm,region=factor(region,levels=c("OTI","URK","ASR","SRM","HST","ASB","KTB","KCA")))

#Visualization
ggplot()+
    geom_point(data=spm,mapping=aes(x=TDH3*0.1,y=PC1,fill=region),
               shape=21,colour="black",size=3,alpha=.7)+
    scale_fill_brewer(palette = "RdBu")+
    theme_classic(base_family = "Helvetica")+
    labs(y="PC1 of candidate traits",x="Total monthly daylight hour in March")+
    geom_smooth(method="lm",data=spm,aes(x=TDH3*0.1,y=PC1),size=.5,linetype=2,colour="black")+
    stat_poly_eq(formula = y ~ x,data=spm,aes(x=TDH3*0.1,y=PC1,
                                              label = paste(stat(eq.label),
                                                            stat(p.value.label),
                                                            sep = "~~~")),parse = TRUE,size=9)+
    theme(legend.position="bottom",
          legend.direction = "horizontal",
          legend.text = element_text(size=rel(1),colour="black"),
          legend.title = element_text(size=rel(1),colour="black"),
          plot.title = element_text(size=rel(2),colour="black"),
          axis.text = element_text(size=rel(1),colour="black"),
          axis.title = element_text(size=rel(2),colour="black"))+
    guides(fill=guide_legend(title="Provenance region",nrow=1,byrow=TRUE))