Below Is the Supporting Information

## Below is the Supporting information

(i) Supplementary figures

Fig. S1Results of Partial Least Squares regression between grassland aboveground productivity and daily unfilted and 15-day running mean filtered soil temperature at 15 cm depth.

Fig. S2Inter-annual variation in air temperature of the 97 national weather stations across the whole Qinghai-Tibetan Plateau during 1997-2011. Red dots mean that there was a significant warming trend during the past 15 years, while the size indicates the magnitude of the increase in air temperature. Grey dots represent stations for which no significant trends for air temperature were found.

Fig.S3Seasonal (a) and monthly (b) differences between soil temperatures at different depths and air temperatures during 1997-2011.

Fig.S4Mean daily and monthly standard deviation of air and soil temperatures at different depths during 1997-2011.

Fig.S5Results of Partial Least Squares regression between grassland aboveground productivity and daily air and soil temperature (at 15 cm depth) from previous September to August. See caption of Fig. 3 for a full explanation.

Fig.S6Monthly mean air and soil temperature at different depths as well as monthly precipitation during the growing season; inter-annual variation in precipitation during the growing season in the past 15 years (a); and daily mean precipitation, evaporation, wind speed and frozen soil depth in our study area (b). The preseason (March-April) period was marked in red dashed lines in this figure.

(ii) Codesfor reproducing the analyses

Code 1:Codes for reproducing the PLS analysis (Fig. 3) in the present study in R language.

#START

windowsFonts(

A=windowsFont("Times New Roman"),

B=windowsFont("Bookman Old Style"),

C=windowsFont("Comic Sans MS"),

D=windowsFont("Symbol")

)

require(chillR)

require(pls)

out_path<-"G:/"

runn_mean <- 15

split_month<-8

weather_data_frame<-read.csv("G:/weather.csv")

bio_data_frame<-read.csv("G:/aug.csv")

PLS_results_path<-paste(out_path,"PLS_chill_force",sep="")

use_Tmean=FALSE

crossvalidate = TRUE

leg1<-182:365

leg<-182:365

lc<-length(leg1)

weather_file <- weather_data_frame

weather_file$Tmax <- suppressWarnings(as.numeric(as.character(weather_file$Tmax)))

weather_file$Tmin <- suppressWarnings(as.numeric(as.character(weather_file$Tmin)))

suppressWarnings(weather_file[which(weather_file$Tmin > weather_file$Tmax),

c("Tmin", "Tmax")] <- NA)

bio_data <- bio_data_frame

runn_mean<-15

Tmin_gaps <- interpolate_gaps(weather_file$Tmin)

weather_file$Tmin <- Tmin_gaps[[1]]

weather_file[, "Tmin_interpolated"] <- 0

suppressWarnings(weather_file[Tmin_gaps[[2]], "Tmin_interpolated"] <- 1)

Tmax_gaps <- interpolate_gaps(weather_file$Tmax)

weather_file$Tmax <- Tmax_gaps[[1]]

weather_file[, "Tmax_interpolated"] <- 0

suppressWarnings(weather_file[Tmax_gaps[[2]], "Tmax_interpolated"] <- 1)

if (use_Tmean) {

Tmean_gaps <- interpolate_gaps(weather_file$Tmean)

weather_file$Tmean <- Tmean_gaps[[1]]

weather_file[, "Tmean_interpolated"] <- 0

weather_file[Tmean_gaps[[2]], "Tmean_interpolated"] <- 1

} else {

weather_file[, "Tmean"] <- (weather_file$Tmax + weather_file$Tmin)/2

}

weather_file[which(weather_file$Month <= split_month), "Season"] <- weather_file[which(weather_file$Month <=

split_month), "Year"]

weather_file[which(weather_file$Month > split_month), "Season"] <- weather_file[which(weather_file$Month >

split_month), "Year"] + 1

weather_file[, "Date"] <- weather_file$Month * 100 + weather_file$Day

weather_file[, "JDay"] <- strptime(paste(weather_file$Month,

"/", weather_file$Day, "/", weather_file$Year, sep = ""),

"%m/%d/%Y")$yday + 1

ww <- weather_file[, "Tmean"]

rr <- weather_file[, "Tmean"]

for (dd in 1:length(ww)) {

if (dd < ceiling(runn_mean/2)) {

rr[dd] <- mean(ww[1:(dd + floor(runn_mean/2))])

}

if ((dd >= ceiling(runn_mean/2)) & (dd <= length(ww) -

ceiling(runn_mean/2))) {

rr[dd] <- mean(ww[(dd - floor(runn_mean/2)):(dd +

floor(runn_mean/2))])

}

if (dd > (length(ww) - ceiling(runn_mean/2))) {

rr[dd] <- mean(ww[(dd - floor(runn_mean/2)):length(ww)])

}

}

weather_file[, "runn"] <- rr

VIP <- function(object) {

if (object$method != "oscorespls")

stop("Only implemented for orthogonal scores algorithm. Refit with 'method = \"oscorespls\"'")

if (nrow(object$Yloadings) > 1)

stop("Only implemented for single-response models")

SS <- c(object$Yloadings)^2 * colSums(object$scores^2)

Wnorm2 <- colSums(object$loading.weights^2)

SSW <- sweep(object$loading.weights^2, 2, SS/Wnorm2,"*")

sqrt(nrow(SSW) * apply(SSW, 1, cumsum)/cumsum(SS))

}

pls_ncomp <- function(indep, dep, threshold = 30) {

dat <- data.frame(dep)

dat$runn <- indep

if (length(dep) > 15) {

suppressWarnings(pls_out <- plsr(dep ~ runn, data = dat,

ncomp = 10, validation = "CV"))

suppressWarnings(pls_cv <- crossval(pls_out, segments = 10))

res <- data.frame(ncomp = c(1:10), explvar = explvar(pls_out),

cumul = NA)

res$cumul[1] <- res$explvar[1]

for (i in 2:nrow(res)) {

res$cumul[i] <- res$cumul[i - 1] + res$explvar[i]

}

ncomp <- which(res$cumul > threshold)[1]

if (is.na(ncomp))

ncomp <- 10

} else ncomp <- 2

return(ncomp)

}

seasons <- unique(weather_file$Season)

for (yy in seasons) {

yearweather <- weather_file[which(weather_file$Season ==

yy), ]

weathervector <- c(yearweather$runn[1:365])

if (yy == seasons[1])

year_res <- weathervector

else year_res <- rbind(year_res, weathervector)

if (nrow(yearweather) == 365) {

labdates <- yearweather$Date

labJdates <- yearweather$JDay

}

}

colnames(year_res) <- c(paste("runn_", 1:365, sep = ""))

year_res <- cbind(Season = seasons, year_res)

data_set <- year_res

full_seasons <- which(!is.na(rowMeans(data_set)))

data_set <- data_set[full_seasons, ]

newseasons <- data_set[, "Season"]

suppressWarnings(bio_data <- bio_data[which(bio_data[, "Year"] %in%

newseasons), ])

suppressWarnings(bio_data <- bio_data[which(!is.na(as.numeric(as.character(bio_data$pheno)))),

])

suppressWarnings(bio <- as.numeric(as.character(bio_data$pheno)))

indep <- as.matrix(data_set[which(data_set[, "Season"] %in%

bio_data$Year), ])

indep <- indep[, 2:ncol(indep)]

if (crossvalidate) {

ncomp <- pls_ncomp(indep = indep, dep = bio)

}else ncomp = crossvalidate

sdindep <- apply(indep, 2, sd)

sdindep[which(sdindep == 0)] <- 1

PLS_output <- plsr(bio ~ indep[,leg1], ncomp, validation = "none",

method = "oscorespls", scale = sdindep[leg1])

out <- data.frame()

out[1:lc, "Date"] <- labdates[leg1]

out[1:lc, "JDay"] <- labJdates[leg1]

out[1:lc, "Coeff"] <- coef(PLS_output)

out[1:lc, "VIP"] <- VIP(PLS_output)[ncomp, ]

colnames(out) <- c("Date", "JDay", "Coefficient", "VIP")

if (is.na(out[1, "VIP"])) {

out[, "VIP"] <- 0

out[, "coef"] <- 0

nothing <- TRUE

}else {

nothing <- FALSE

}

color_bar_maker <- function(columns_yn, columns_quant, threshold,

col1, col2, col3) {

color_bars <- c(rep(NA, length(columns_yn)))

color_bars[which(columns_yn >= threshold & columns_quant < 0)] <- col1

color_bars[which(columns_yn >= threshold & columns_quant >=0)] <- col2

color_bars[which(!columns_yn >= threshold)] <- col3

return(color_bars)

}

strptime(paste("01/", split_month + 1, "/2001", sep = ""),

format = "%d/%m/%Y") - 86400

leg <- strptime(strptime(paste("01/", split_month + 1, "/2001",

sep = ""), format = "%d/%m/%Y") - 86400 + leg * 86400, "%Y-%m-%d")

tick_marks <- leg[sapply(strsplit(as.character(leg), "-"), "[", 3) == "01"]

tick_labels <- as.Date(tick_marks, origin = "1999-1-2")

tick_labels <- as.POSIXlt(tick_labels)$mon

tick_labels <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun","Jul", "Aug", "Sep", "Oct", "Nov", "Dec")[tick_labels + 1]

mean_clim <- colMeans(indep[,leg1], na.rm = TRUE)

dev_clim <- apply(indep[,leg1], 2, sd, na.rm = TRUE)

tick_label_pos <- leg[sapply(strsplit(as.character(leg), "-"), "[", 3) == "15"]

bmp(paste(PLS_results_path, ".bmp", sep = ""), width = 1000, height = 1000, pointsize = 10)

par(mfcol = c(2, 1))

par(mar = c(6.1, 9.8, 6.1, 2.1))

par(mgp = c(4, 1.5, 0))

color_bars <- color_bar_maker(out[1:lc, "VIP"], out[1:lc, "VIP"], threshold = 0.8,

col1 = "DARK BLUE", col2 = "DARK BLUE",

col3 = "DARK GREY")

plot(leg, out[1:lc, "VIP"], main = "VIP", xlab = NA, ylab = NA,

xaxs = "i", xaxt = "n", yaxs = "i", cex.lab = 4, cex.axis = 4,

cex.main = 5, type = "h", col = color_bars, lwd = 6,family="A")

tick_marks <- grconvertX(tick_marks, from = "user", to = "user")

X_coords <- grconvertX(leg, from = "user", to = "user")

ticks_labels <- grconvertX(tick_label_pos, from = "user", to = "user")

axis(1, lwd.ticks = 3, at = tick_marks, labels = FALSE, cex.axis = 4, padj = 1)

axis(2, lwd.ticks = 3, labels = FALSE)

box(which = "plot", lwd = 3)

axis(1, lwd.ticks = 0, at = ticks_labels, labels = tick_labels, cex.axis = 4, padj = 1,family="A")

mtext(side = 2, text = "VIP", line = 6, cex = 4,family="A")

if (nothing) {

text(x = (max(X_coords) - min(X_coords))/2 + min(X_coords),

y = 0.5, adj = 0.5, labels = ("No yield produced in any model run"),

cex = 2)

}

tick_marks2 <- leg[sapply(strsplit(as.character(leg), "-"), "[", 3) == "01"]

tick_marks2<-tick_marks2[2:6]

ticks_labels2 <- grconvertX(tick_marks2, from = "user",to = "user")

abline(v=ticks_labels2,lwd=2,lty="dashed",col="DARK GREY")

color_bars <- color_bar_maker(out[1:lc, "VIP"], out[1:lc, "Coefficient"],

threshold = 0.8, col1 = "RED", col2 = "DARK GREEN", col3 = "DARK GREY")

plot(leg, out[1:lc, "Coefficient"], main = "Model coefficients",

ylab = NA, xlab = NA, xaxs = "i", xaxt = "n", yaxs = "i",

cex.lab = 4, cex.axis = 4, cex.main = 5, type = "h",

col = color_bars, lwd = 6,family="A")

axis(1, lwd.ticks = 3, at = tick_marks, labels = FALSE, cex.axis = 4, padj = 1)

axis(2, lwd.ticks = 3, labels = FALSE)

box(which = "plot", lwd = 3)

axis(1, lwd.ticks = 0, at = ticks_labels, labels = tick_labels, cex.axis = 4, padj = 1,family="A")

mtext(side = 2, text = "Model coefficients", line = 6, cex = 4,family="A")

if (nothing) {

text(x = (max(X_coords) - min(X_coords))/2 + min(X_coords),

y = 0.5, adj = 0.5, labels = ("No yield produced in any model run"),

cex = 2)

}

tick_marks2 <- leg[sapply(strsplit(as.character(leg), "-"),"[", 3) == "01"]

tick_marks2<-tick_marks2[2:6]

ticks_labels2 <- grconvertX(tick_marks2, from = "user",to = "user")

abline(v=ticks_labels2,lwd=2,lty="dashed",col="DARK GREY")

dev.off()

p<-predict(PLS_output)

ll<-length(bio)

df<-data.frame(obs=bio,pred=p[(ll*(ncomp-1)+1):(ll*ncomp)])

rmse <- function(obs, pred) sqrt(mean((obs-pred)^2))

rmse<-rmse(df$obs,df$pred)

write.csv(rmse,"G:/rmse.csv",row.names = FALSE)

write.csv(out, paste("G:/", "VIP_MC.csv", sep = ""), row.names = FALSE)

# END

Code 2: Codes for the analysis and drawing of daily and monthly temperature in Fig. 2.

#Taking Soil layer 20cm as an example

#START

windowsFonts(

A=windowsFont("Times New Roman"),

B=windowsFont("Bookman Old Style"),

C=windowsFont("Comic Sans MS"),

D=windowsFont("Symbol")

)

require(chillR)

out_path<-"G:/"

runn_mean <- 15

split_month<-8

weather_data_frame<-read.csv("G:/weather.csv")

bio_data_frame<-read.csv("G:/aug.csv")

PLS_results_path<-paste(out_path,"climate_index",sep="")

use_Tmean = FALSE

main_lable<-"Soil layer 20cm"

soil_air<-read.csv("G:/soil_air.csv")

#soil_air<-soil_air[,1:4]

#soil_air<-soil_air[,c(1:3,5)]

soil_air<-soil_air[,c(1:3,6)] #Soil layer 20cm

colnames(soil_air)<-c("Year","Month","Day","Td")

require(pls)

weather_file <- weather_data_frame

weather_file$Tmax <- suppressWarnings(as.numeric(as.character(weather_file$Tmax)))

weather_file$Tmin <- suppressWarnings(as.numeric(as.character(weather_file$Tmin)))

suppressWarnings(weather_file[which(weather_file$Tminweather_file$Tmax),

c("Tmin", "Tmax")] <- NA)

bio_data <- bio_data_frame

runn_mean<-15

Tmin_gaps <- interpolate_gaps(weather_file$Tmin)

weather_file$Tmin <- Tmin_gaps[[1]]

weather_file[, "Tmin_interpolated"] <- 0

suppressWarnings(weather_file[Tmin_gaps[[2]], "Tmin_interpolated"] <- 1)

Tmax_gaps <- interpolate_gaps(weather_file$Tmax)

weather_file$Tmax <- Tmax_gaps[[1]]

weather_file[, "Tmax_interpolated"] <- 0

suppressWarnings(weather_file[Tmax_gaps[[2]], "Tmax_interpolated"] <- 1)

if (use_Tmean) {

Tmean_gaps <- interpolate_gaps(weather_file$Tmean)

weather_file$Tmean <- Tmean_gaps[[1]]

weather_file[, "Tmean_interpolated"] <- 0

weather_file[Tmean_gaps[[2]], "Tmean_interpolated"] <- 1

} else {

weather_file[, "Tmean"] <- (weather_file$Tmax + weather_file$Tmin)/2

}

weather_file[which(weather_file$Month <= split_month), "Season"] <- weather_file[which(weather_file$Month <=

split_month), "Year"]

weather_file[which(weather_file$Monthsplit_month), "Season"] <- weather_file[which(weather_file$Month

split_month), "Year"] + 1

weather_file[, "Date"] <- weather_file$Month * 100 + weather_file$Day

weather_file[, "JDay"] <- strptime(paste(weather_file$Month,

"/", weather_file$Day, "/",weather_file$Year, sep = ""), "%m/%d/%Y")$yday + 1

ww <- weather_file[, "Tmean"]

rr <- weather_file[, "Tmean"]

for (dd in 1:length(ww)) {

if (dd < ceiling(runn_mean/2)) {

rr[dd] <- mean(ww[1:(dd + floor(runn_mean/2))])

}

if ((dd >= ceiling(runn_mean/2)) & (dd <= length(ww) -

ceiling(runn_mean/2))) {

rr[dd] <- mean(ww[(dd - floor(runn_mean/2)):(dd +

floor(runn_mean/2))])

}

if (dd > (length(ww) - ceiling(runn_mean/2))) {

rr[dd] <- mean(ww[(dd - floor(runn_mean/2)):length(ww)])

}

}

weather_file[, "runn"] <- rr

seasons <- unique(weather_file$Season)

for (yy in seasons) {

yearweather <- weather_file[which(weather_file$Season ==

yy), ]

weathervector <- c(yearweather$runn[1:365])

if (yy == seasons[1])

year_res <- weathervector

elseyear_res <- rbind(year_res, weathervector)

if (nrow(yearweather) == 365) {

labdates <- yearweather$Date

labJdates <- yearweather$JDay

}

}

colnames(year_res) <- c(paste("runn_", 1:365, sep = ""))

year_res <- cbind(Season = seasons, year_res)

data_set <- year_res

full_seasons <- which(!is.na(rowMeans(data_set)))

data_set <- data_set[full_seasons, ]

newseasons <- data_set[, "Season"]

suppressWarnings(bio_data <- bio_data[which(bio_data[, "Year"] %in%

newseasons), ])

suppressWarnings(bio_data <- bio_data[which(!is.na(as.numeric(as.character(bio_data$pheno)))),

])

suppressWarnings(bio <- as.numeric(as.character(bio_data$pheno)))

indep <- as.matrix(data_set[which(data_set[, "Season"] %in% bio_data$Year), ])

indep <- indep[, 2:ncol(indep)]

sdindep <- apply(indep, 2, sd)

sdindep[which(sdindep == 0)] <- 1

lc <- 365

leg <- 1:365

strptime(paste("01/", split_month + 1, "/2001", sep = ""), format = "%d/%m/%Y") - 86400

leg <- strptime(strptime(paste("01/", split_month + 1, "/2001",

sep = ""), format = "%d/%m/%Y") - 86400 + leg * 86400, "%Y-%m-%d")

tick_marks <- leg[sapply(strsplit(as.character(leg), "-"), "[", 3) == "01"]

tick_labels <- as.Date(tick_marks, origin = "1999-1-2")

tick_labels <- as.POSIXlt(tick_labels)$mon

tick_labels <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun",

"Jul", "Aug", "Sep", "Oct", "Nov", "Dec")[tick_labels + 1]

mean_clim <- colMeans(indep, na.rm = TRUE)

dev_clim <- apply(indep, 2, sd, na.rm = TRUE)

write.csv(mean_clim,"G:/mean_clim.csv",row.names=FALSE)

mean_clim2<-read.csv("G:/mean_clim.csv")

mean_clim_for_Jday<-data.frame("mean"=rep(NA,365),"Date"=rep(NA,365),"Date2"=rep(NA,365))

mean_clim_for_Jday[,"mean"]<-mean_clim2[1:365,1]

mean_clim_for_Jday[,"Date"]<-leg[1:365]$yday+1

mean_clim_for_Jday[1:122,"Date2"]<-mean_clim_for_Jday[1:122,"Date"]-243

mean_clim_for_Jday[123:365,"Date2"]<-mean_clim_for_Jday[123:365,"Date"]+122

mean_clim_for_Jday<-mean_clim_for_Jday[which(mean_clim_for_Jday$Date2 %in% c(182:212)),]

located_date<-mean_clim_for_Jday[min(which(mean_clim_for_Jday$mean>0)),"Date"]-59

bmp(paste(PLS_results_path, ".bmp", sep = ""), width = 2000, height = 2000, pointsize = 20)

par(mfcol = c(3, 1))

par(mar = c(6.1, 9.8, 6.1, 2.1))

par(mgp = c(4, 1.5, 0))

tick_label_pos <- leg[sapply(strsplit(as.character(leg), "-"), "[", 3) == "15"]

plot(leg, mean_clim[1:lc], main = main_lable, ylab = NA,

xlab = NA, xaxs = "i", yaxs = "i", xaxt = "n", cex.lab = 4,

cex.axis = 4, cex.main = 5, type = "l", lwd = 3, col = "BLACK",

ylim = c(–-15,20),family="A")

tick_marks <- grconvertX(tick_marks, from = "user", to = "user")

X_coords <- grconvertX(leg, from = "user", to = "user")

ticks_labels <- grconvertX(tick_label_pos, from = "user", to = "user")

arrows(X_coords, mean_clim[1:lc] + dev_clim[1:lc], X_coords,

mean_clim[1:lc] - dev_clim[1:lc], angle = 90, code = 3,

lwd = 6, length = 0, col = "DARK GREY")

lines(leg, mean_clim[1:lc], lwd = 3)

axis(1, lwd.ticks = 3, at = tick_marks, labels = FALSE, cex.axis = 4,padj = 1,family="A")

axis(2, lwd.ticks = 3, labels = FALSE)

box(which = "plot", lwd = 3)

axis(1, lwd.ticks = 0, at = ticks_labels, labels = tick_labels, cex.axis = 4, padj = 1,family="A")

mtext(side = 2, text = expression("Mean temperature ("^"o" * "C)"), line = 5, cex = 3,family="A")

abline(h=0,lwd=4,lty="dashed",col="DARK GREY")

located_date2<-as.character(located_date)

tick_marks2 <- leg[sapply(strsplit(as.character(leg), "-"), "[", 3) == located_date2]

tick_marks3<-tick_marks2[7]

ticks_labels3 <- grconvertX(tick_marks3, from = "user",to = "user")

abline(v=ticks_labels3,lwd=4,lty="dashed",col="DARK GREY")

mtext(paste(located_date," Mar",sep=""), line=-5, side=1,at =(ticks_labels3),

cex = 2.3,family="A")

monthly_data<-weather_file[which(weather_file$Season %in% bio_data$Year),]

month_mean<-data.frame()

tempo<-data.frame()

for (i in 1:12){

monthly<-monthly_data[which(monthly_data$Month==i),]

year<-unique(monthly$Year)

q=1

for (y in year){

month_yr<-monthly[which(monthly$Year==y),]

tempo[q,1]<-mean(month_yr$Tmean)

q=q+1}

month_mean[i,1]<-i

month_mean[i,2]<-mean(tempo[,1])

month_mean[i,3]<-sd(tempo[,1])

}

colnames(month_mean)<-c("Month","Values","Sd")

up<-month_mean$Values[c(9:12,1:8)]+month_mean$Sd[c(9:12,1:8)]

down<-month_mean$Values[c(9:12,1:8)]-month_mean$Sd[c(9:12,1:8)]

lines(x=rep(tick_label_pos[1],2),y=c(down[1],up[1]),lwd=8,col="DARK GREEN")

lines(x=rep(tick_label_pos[2],2),y=c(down[2],up[2]),lwd=8,col="DARK GREEN")

lines(x=rep(tick_label_pos[3],2),y=c(down[3],up[3]),lwd=8,col="DARK GREEN")

lines(x=rep(tick_label_pos[4],2),y=c(down[4],up[4]),lwd=8,col="DARK GREEN")

lines(x=rep(tick_label_pos[5],2),y=c(down[5],up[5]),lwd=8,col="DARK GREEN")

lines(x=rep(tick_label_pos[6],2),y=c(down[6],up[6]),lwd=8,col="DARK GREEN")

lines(x=rep(tick_label_pos[7],2),y=c(down[7],up[7]),lwd=8,col="DARK GREEN")

lines(x=rep(tick_label_pos[8],2),y=c(down[8],up[8]),lwd=8,col="DARK GREEN")

lines(x=rep(tick_label_pos[9],2),y=c(down[9],up[9]),lwd=8,col="DARK GREEN")

lines(x=rep(tick_label_pos[10],2),y=c(down[10],up[10]),lwd=8,col="DARK GREEN")

lines(x=rep(tick_label_pos[11],2),y=c(down[11],up[11]),lwd=8,col="DARK GREEN")

lines(x=rep(tick_label_pos[12],2),y=c(down[12],up[12]),lwd=8,col="DARK GREEN")

points(tick_label_pos, month_mean$Values[c(9:12,1:8)],col = "RED",cex=2.5,pch=19)

soil_air[which(soil_air$Month <= split_month), "Season"] <- soil_air[which(soil_air$Month <=

split_month), "Year"]

soil_air[which(soil_air$Monthsplit_month), "Season"] <- soil_air[which(soil_air$Month

split_month), "Year"] + 1

soil_air<-soil_air[which(soil_air$Season %in% bio_data$Year),]

soil_air_mean<-data.frame()

tempo2<-data.frame()

for (i in 1:12){

soil<-soil_air[which(soil_air$Month==i),]

year<-unique(soil$Year)

q=1

for (y in year){

soil_yr<-soil[which(soil$Year==y),]

tempo2[q,1]<-mean(soil_yr$Td)

q=q+1}

soil_air_mean[i,1]<-i

soil_air_mean[i,2]<-mean(tempo2[,1])

soil_air_mean[i,3]<-sd(tempo2[,1])

}

colnames(soil_air_mean)<-c("Month","Values","Sd")

tick_label_pos3 <- leg[sapply(strsplit(as.character(leg), "-"), "[", 3) == "15"]

up2<-soil_air_mean$Values[c(9:12,1:8)]+soil_air_mean$Sd[c(9:12,1:8)]

down2<-soil_air_mean$Values[c(9:12,1:8)]-soil_air_mean$Sd[c(9:12,1:8)]

lines(x=rep(tick_label_pos3[1],2),y=c(down2[1],up2[1]),lwd=8,col="DARK GREY")

lines(x=rep(tick_label_pos3[2],2),y=c(down2[2],up2[2]),lwd=8,col="DARK GREY")

lines(x=rep(tick_label_pos3[3],2),y=c(down2[3],up2[3]),lwd=8,col="DARK GREY")

lines(x=rep(tick_label_pos3[4],2),y=c(down2[4],up2[4]),lwd=8,col="DARK GREY")

lines(x=rep(tick_label_pos3[5],2),y=c(down2[5],up2[5]),lwd=8,col="DARK GREY")

lines(x=rep(tick_label_pos3[6],2),y=c(down2[6],up2[6]),lwd=8,col="DARK GREY")

lines(x=rep(tick_label_pos3[7],2),y=c(down2[7],up2[7]),lwd=8,col="DARK GREY")

lines(x=rep(tick_label_pos3[8],2),y=c(down2[8],up2[8]),lwd=8,col="DARK GREY")

lines(x=rep(tick_label_pos3[9],2),y=c(down2[9],up2[9]),lwd=8,col="DARK GREY")

lines(x=rep(tick_label_pos3[10],2),y=c(down2[10],up2[10]),lwd=8,col="DARK GREY")

lines(x=rep(tick_label_pos3[11],2),y=c(down2[11],up2[11]),lwd=8,col="DARK GREY")

lines(x=rep(tick_label_pos3[12],2),y=c(down2[12],up2[12]),lwd=8,col="DARK GREY")

points(tick_label_pos3, soil_air_mean$Values[c(9:12,1:8)],col = rgb(0,0.5,1),cex=2.5,pch=19)

gl<-round(mean(soil_air_mean$Values),digits=2)

tick_marks4 <- leg[sapply(strsplit(as.character(leg), "-"), "[", 3) == "01"]

tick_marks4<-tick_marks4[11]

ticks_labels4 <- grconvertX(tick_marks4, from = "user",to = "user")

mtext(paste("Annual mean (Soil-Air)= ",gl,sep=""), line=-19, side=1,at =(ticks_labels4),

cex = 2.3,family="A")

dev.off()

write.csv(month_mean,"G:/month_mean.csv",row.names=FALSE)

write.csv(soil_air_mean,"G:/soil_air_mean.csv",row.names=FALSE)

# END

Code 3: Codes for the Kriging interpolation in Fig. 4.

#Taking Soil layer 20cm as an example

#START

windowsFonts(

A=windowsFont("Times New Roman"),

B=windowsFont("Bookman Old Style"),

C=windowsFont("Comic Sans MS"),

D=windowsFont("Symbol")

)

library(chillR)

library(fields)

weather<-read.csv("G:/weather.csv")

weather<-fix_weather(weather)

file_path<-"G:/"

weather_data_frame = weather$weather

split_month = 8

pheno = read.csv("G:/aug.csv")

use_Tmean = FALSE

Start_JDay_chill = 65

End_JDay_chill = 118

Start_JDay_heat = 123

End_JDay_heat = 243

outpath = file_path

file_name = "pheno_trend_plot"

plot_title = "20cm soil temperature"

ylabel = NA

xlabel = NA

legend_label = NA

image_type = "png"

colorscheme = "normal"

fonttype="serif"

weather_file <- weather_data_frame

weather_file$Tmax <- suppressWarnings(as.numeric(as.character(weather_file$Tmax)))

weather_file$Tmin <- suppressWarnings(as.numeric(as.character(weather_file$Tmin)))

suppressWarnings(weather_file[which(weather_file$Tminweather_file$Tmax),

c("Tmin", "Tmax")] <- NA)

Tmin_gaps <- interpolate_gaps(weather_file$Tmin)

weather_file$Tmin <- Tmin_gaps[[1]]

weather_file[, "Tmin_interpolated"] <- 0

suppressWarnings(weather_file[Tmin_gaps[[2]], "Tmin_interpolated"] <- 1)

Tmax_gaps <- interpolate_gaps(weather_file$Tmax)

weather_file$Tmax <- Tmax_gaps[[1]]

weather_file[, "Tmax_interpolated"] <- 0

suppressWarnings(weather_file[Tmax_gaps[[2]], "Tmax_interpolated"] <- 1)

if (use_Tmean) {

Tmean_gaps <- interpolate_gaps(weather_file$Tmean)

weather_file$Tmean <- Tmean_gaps[[1]]

weather_file[, "Tmean_interpolated"] <- 0

weather_file[Tmean_gaps[[2]], "Tmean_interpolated"] <- 1

} else {

weather_file[, "Tmean"] <- (weather_file$Tmax + weather_file$Tmin)/2

}

weather_file[which(weather_file$Month <= split_month), "Season"] <- weather_file[which(weather_file$Month <=

split_month), "Year"]

weather_file[which(weather_file$Monthsplit_month), "Season"] <- weather_file[which(weather_file$Month

split_month), "Year"] + 1

weather_file[, "Date"] <- weather_file$Month * 100 + weather_file$Day

weather_file[, "JDay"] <- strptime(paste(weather_file$Month,

"/", weather_file$Day, "/", weather_file$Year, sep = ""),

"%m/%d/%Y")$yday + 1

ww <- weather_file[, "Tmean"]

rr <- weather_file[, "Tmean"]

sea <- unique(weather_file$Season)

res <- data.frame(Season = sea, Chill_Tmean = NA, Heat_Tmean = NA, Year_Tmean = NA)

if (Start_JDay_chillEnd_JDay_chill)

chill_days <- c(Start_JDay_chill:366, 1:End_JDay_chill) else {

chill_days <- Start_JDay_chill:End_JDay_chill}

if (Start_JDay_heatEnd_JDay_heat)

heat_days <- c(Start_JDay_heat:366, 1:End_JDay_heat) else {

heat_days <- Start_JDay_heat:End_JDay_heat}

for (s in sea) {

res[which(res$Season == s), "Chill_Tmean"] <- mean(weather_file[which(weather_file$Season ==

s & weather_file$JDay %in% chill_days), "Tmean"])

res[which(res$Season == s), "Heat_Tmean"] <- mean(weather_file[which(weather_file$Season ==

s & weather_file$JDay %in% heat_days), "Tmean"])

res[which(res$Season == s), "Year_Tmean"] <- mean(weather_file[which(weather_file$Season ==

s), "Tmean"])

}

for (i in 1:nrow(pheno)) {

if (pheno[i, 1] %in% res$Season) {

pheno[i, "Chill_Tmean"] <- res[which(res$Season ==

pheno[i, 1]), "Chill_Tmean"]

pheno[i, "Heat_Tmean"] <- res[which(res$Season ==

pheno[i, 1]), "Heat_Tmean"]

pheno[i, "Year_Tmean"] <- res[which(res$Season ==

pheno[i, 1]), "Year_Tmean"]

}

}

Make_date <- function(Jday) {

leg <- strptime(strptime(paste(Jday, "/2001", sep = ""),

format = "%j/%Y"), "%Y-%m-%d")

comps <- strsplit(as.character(leg), "-")

JM <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul",

"Aug", "Sep", "Oct", "Nov", "Dec")[as.numeric(comps[[1]][2])]

return(paste(comps[[1]][3], JM))

}

if (is.na(xlabel))

xlabel <- paste("Mean temperature during ",

Make_date(Start_JDay_chill), " - ", Make_date(End_JDay_chill), " (°C)", sep = "")

if (is.na(ylabel))

ylabel <- paste("Mean temperature during ",

Make_date(Start_JDay_heat), " - ", Make_date(End_JDay_heat), " (°C)", sep = "")

pheno[, 2] <- as.numeric(as.character(pheno[, 2]))

pheno <- pheno[which(!is.na(pheno[, "pheno"])), ]

pheno <- pheno[which(!is.na(pheno[, "Chill_Tmean"])), ]

pheno <- pheno[which(!is.na(pheno[, "Heat_Tmean"])), ]

pheno <- pheno[which(!is.na(pheno[, "Year_Tmean"])), ]

k <- Krig(x = as.matrix(pheno[, c("Chill_Tmean", "Heat_Tmean")]), Y = pheno$pheno)

legend_label<-expression("ANPP (g m"^"-2" * " yr"^"-1"*")")

if (image_type == "tiff")

tiff(paste(outpath, file_name, ".tif", sep = ""), width = 1400,

height = 1500) else {

png(paste(outpath, file_name, ".png", sep = ""), width = 1400,height = 1500)}

par(family = fonttype)

par(mai = c(2.5, 2.5, 1.5, 1),lwd=4)

if (colorscheme == "bw")

colors <- colorRampPalette(c("#F9F9F9", "#343434"))(255) else {

colors <- tim.colors(255)}

surface(k, add.legend=FALSE,col = colors, type = "C", xaxt = "n", yaxt = "n",

xlab=NA,ylab=NA, cex.lab =5, cex.axis = 5, labcex = 5, asp = 1,

legend.mar=10,axis.args = list(cex.axis = 4,lwd.ticks=4,family="A"))

mtext(text = plot_title, side = 3, cex = 6, line = 1.7)

points(pheno[, c("Chill_Tmean", "Heat_Tmean")], pch = 16,cex=2.5,col="DARK BLUE")

axis(1,lwd.ticks=4,cex.axis=5,padj=0.8)

axis(2,lwd.ticks=4,cex.axis=5,padj=-0.3)

mtext(text = xlabel, side = 1, cex = 5.5, line = 7.7)

mtext(text = ylabel, side = 2, cex = 5.5, line = 7)

dev.off()

pheno <- pheno[, which(!is.na(pheno[1, ]))]

write.csv(pheno, paste(outpath, file_name, ".csv", sep = ""), row.names = FALSE)

#END

Code 4: Codes for the linear regression analysis in Fig. 5.

#Taking Soil layer 20cm as an example

#START

windowsFonts(

A=windowsFont("Times New Roman"),

B=windowsFont("Bookman Old Style"),

C=windowsFont("Comic Sans MS"),

D=windowsFont("Symbol")

)

par(mfrow = c(2, 1))

data<-read.csv("C:/Users/Administrator/Desktop/data.csv")

library(ggplot2)

library(gcookbook)

hwp1<-ggplot(data[,c("Chill_Tmean","pheno")],aes(x=Chill_Tmean,y=pheno))+geom_point(colour="RED")

hwp1+xlab("Mean temp (6 Mar-28 Apr)")+ylab(expression("ANPP (g m"^"-2" * " yr"^"-1"*")"))+xlim(0.8,3.1)+

theme(axis.title.x=element_text(size=14,family="A",colour="WHITE"),axis.title.y=element_text(size=14,family="A"),

axis.text.x=element_text(size=12,family="A"),axis.text.y=element_text(size=12,family="A"))+guides(colour=FALSE)+

stat_smooth(method=lm,colour="RED")

ggsave("C:/Users/Administrator/Desktop/hwp1.png",width=5.4,height=5.4,unit="cm",dpi=1000)

hwp2<-ggplot(data[,c("Heat_Tmean","pheno")],aes(x=Heat_Tmean,y=pheno))+geom_point(colour="BLUE")

hwp2+xlab("Mean temp (3 May-31 Aug)")+ylab(expression("ANPP (g m"^"-2" * " yr"^"-1"*")"))+xlim(12.8,15.5)+

theme(axis.title.x=element_text(size=14,family="A",colour="WHITE"),axis.title.y=element_text(size=14,family="A"),

axis.text.x=element_text(size=12,family="A"),axis.text.y=element_text(size=12,family="A"))+guides(colour=FALSE)+

stat_smooth(method=lm,colour="BLUE")

ggsave("C:/Users/Administrator/Desktop/hwp2.png",width=5.4,height=5.4,unit="cm",dpi=1000)

#END

#Note: We only listed parts of all the codes used in the present study. If you like them #and want to know more, please contact us (E-mail: ), #and we are quite pleased to share what we have with you. Thanks again for your #reading of our long codes in R language, and best wishes to you.

1