Long-Term Trends in CO2 Concentrations & Carbon Isotope Composition

Dominick Cifelli

May 3, 2023


library(isocalcR)
library(tidyr)
library(dplyr)

Attaching package: 'dplyr'
The following objects are masked from 'package:stats':

    filter, lag
The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union
library(cowplot)
library(ggplot2)
library(ggthemes)

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

    theme_map
library(gt)
library(Metrics)
library(car)
Loading required package: carData

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

    recode

Introduction

The Earth’s climate is affected by the rampant release of greenhouse gases, especially carbon dioxide (CO2), into the atmosphere. This surge in CO2 levels carries profound consequences for the world’s ecosystems and poses formidable challenges for maintaining resources sustainably and conserving our environment. Grasping the shifts in atmospheric CO2 concentrations (Ca) and its isotopic composition (δ13C) across time is pivotal for gauging the scale and implications of human-triggered climate alterations. Over the last century, the combustion of fossil fuels, deforestation, and other human endeavors have driven a substantial uptick in atmospheric CO2 levels. The industrial revolution marked the genesis of a staggering CO2 rise, an ascent that has gained pace in recent decades. The current CO2 concentration eclipses pre-industrial levels by approximately 50%, soaring to around 420 parts per million (ppm). This unparalleled surge in CO2 holds momentous implications for the Earth’s climate system and the functioning of terrestrial carbon absorption.

Carbon isotopes unveil precious insights into the mechanisms propelling changes in atmospheric CO2 concentrations. Carbon exists in two stable forms, 12C and 13C, differing in their neutron count. The proportion of these isotopes in atmospheric CO2 (δ13C.atm) is influenced by myriad factors, encompassing fossil fuel emissions, the dynamics of the carbon cycle, and the way plants respond to shifting CO2 levels. Unveiling the temporal patterns in δ13C.atm is crucial for disentangling the contributions from different CO2 sources and sinks, and for evaluating the repercussions of human activities on global carbon circulation. Beyond being a tracer of atmospheric CO2, δ13C composition furnishes crucial insights into the internal processes of plants. During photosynthesis, plants favor the lighter isotope, 12C, resulting in reduced δ13C values as atmospheric CO2 levels climb. Variations in δ13C composition in plant tissues, like wood, can mirror changes in carbon assimilation paths, water efficiency, and nutrient availability. Consequently, delving into the nexus between δ13C values of woody organic matter and environmental elements can unveil plant adaptations to shifting CO2 levels and climatic shifts. Despite the growing reservoir of research on atmospheric CO2 levels and carbon isotope composition, there persists a need for all-encompassing inquiries that synthesize long-term data and employ advanced statistical models. Such investigations can offer a more nuanced understanding of the temporal dynamics and ramifications of these variables. Hence, the ambition of this study is to scrutinize the temporal trends in atmospheric CO2 concentrations and carbon isotope composition over the last 80 years, leveraging data from Belmecheri and Lavergne (2020) and the isocalcR toolkit developed by Mathias and Hudiburg. Furthermore, I aspire to probe the interplay between δ13C composition of woody organic matter and environmental components, such as water-use efficiency and climatic parameters.

The isocalcR toolkit bridges the gap for standardized tools to analyze stable isotope data from wood and leaf organic matter. Equipped with an array of functions and recommended reference data, isocalcR facilitates the computation of common isotope-based physiological indices. It calculates indices like leaf carbon isotope discrimination, the intercellular-to-atmospheric CO2 ratio, and intrinsic water use efficiency—each fraction corresponding to key phenological processes such as carbon assimilation (iWUE Mesophyll) and light-dependent processes (iWUE Photorespiration). Additionally, isocalcR incorporates endorsed atmospheric CO2 and δ13CO2 data, assuring precision and currentness in calculations.

Therefore, my aim is to compute fractions of water use efficiency using the formulas outlined in Mathias & Hudiburg, subsequently comparing each rendition of WUE. I will also identify trends in atmospheric CO2 concentrations, δ13CO2 composition of woody organic matter, and intrinsic water use efficiency. Subsequently, I will explore the interrelation between δ13CO2 composition of woody organic matter and environmental factors like mean growth temperature (MGT_C) and elevation (Elevation_m). Overall, my endeavor is to unveil trends and correlations across various environmental factors and plant physiology. To accomplish this, i will use the following packages:

Package Description
isocalcR isocalcR Provides functions for calculating isotopic distributions and mass-to-charge ratios for chemical compounds
dplyr Package for data manipulation in R, providing a concise and intuitive syntax for filtering, transforming, and summarizing data
ggplot2 Data visualization tool that allows users to create complex and customized plots with ease
ggthemes Extends the functionality of ggplot2 by providing additional themes and scales for creating visually appealing and professional-looking plots
cowplot Facilitates the creation of complex and publication-ready plots by allowing users to arrange multiple plots into grids and customize their appearance
gt Tool for creating elegant and customizable tables from data frames or other tabular data sources, offering a range of formatting options and support for various output formats
Metrics Provides a collection of commonly used evaluation metrics for machine learning and statistical models, allowing users to assess the performance and accuracy of models
car Offers a wide range of functions for regression modeling and diagnostic testing, including tools for assessing model assumptions, conducting hypothesis tests, and performing regression diagnostics

Methods

Data Exploration & Manipulation

First I read in the data and explore summary statistics


data("CO2data")
data("piru13C")
summary(CO2data)
       yr               Ca           d13C.atm     
 Min.   :   0.0   Min.   :276.2   Min.   :-8.656  
 1st Qu.: 505.5   1st Qu.:277.8   1st Qu.:-6.430  
 Median :1011.0   Median :279.0   Median :-6.420  
 Mean   :1011.0   Mean   :282.7   Mean   :-6.478  
 3rd Qu.:1516.5   3rd Qu.:281.2   3rd Qu.:-6.400  
 Max.   :2022.0   Max.   :415.0   Max.   :-6.280  
summary(piru13C)
      Year          Site             wood.d13C          MGT_C      
 Min.   :1940   Length:223         Min.   :-25.69   Min.   :16.06  
 1st Qu.:1958   Class :character   1st Qu.:-23.44   1st Qu.:17.44  
 Median :1977   Mode  :character   Median :-22.90   Median :17.78  
 Mean   :1977                      Mean   :-22.87   Mean   :17.83  
 3rd Qu.:1995                      3rd Qu.:-22.13   3rd Qu.:18.33  
 Max.   :2014                      Max.   :-20.99   Max.   :19.22  
  Elevation_m        frac  
 Min.   :1033   Min.   :2  
 1st Qu.:1033   1st Qu.:2  
 Median :1060   Median :2  
 Mean   :1099   Mean   :2  
 3rd Qu.:1206   3rd Qu.:2  
 Max.   :1206   Max.   :2

CO2data: yr = Years (0-2022), Ca = Atmospheric CO2 concentration (ppm), d13C.atm = Atmospheric carbion isotoppe (δ13C) ratio

Here I merge the piru13C and CO2data datasets based on the “Year” variable in piru13C and the “yr” variable in CO2data. The all = TRUE argument ensures that all observations from both datasets are included in the merged dataset.


combined_data <- merge(piru13C, CO2data, by.x = "Year", by.y = "yr", all = TRUE)

This creates a new varaible in our dataset that contrains converted wood δ13C composition to d13C values

combined_data$d13C_wood <- combined_data$wood.d13C / 1000

Now I can calculate water use effeciency fractions using the formulas detailed in Mathias & Hudiburg


combined_data$iWUE_simple <- mapply(d13C.to.iWUE,
                                    d13C.plant = combined_data$wood.d13C,
                                    year = combined_data$Year,
                                    elevation = combined_data$Elevation_m,
                                    temp = combined_data$MGT_C,
                                    method = "simple",
                                    tissue = "wood")

combined_data$iWUE_photorespiration <- mapply(d13C.to.iWUE,
                                              d13C.plant = combined_data$wood.d13C,
                                              year = combined_data$Year,
                                              elevation = combined_data$Elevation_m,
                                              temp = combined_data$MGT_C,
                                              method = "photorespiration",
                                              frac = combined_data$frac)

combined_data$iWUE_mesophyll <- mapply(d13C.to.iWUE,
                                       d13C.plant = combined_data$wood.d13C,
                                       year = combined_data$Year,
                                       elevation = combined_data$Elevation_m,
                                       temp = combined_data$MGT_C,
                                       method = "mesophyll",
                                       frac = combined_data$frac)

Here I prepare the data for visualization. Specifically,the data_long dataset will have the columns Year, Ca, d13C.atm, Variable, and Value. The Variable column will contain the names of the original selected columns (d13C_wood, iWUE_simple, iWUE_photorespiration, iWUE_mesophyll), and the Value column will contain the corresponding values.


data_long <- combined_data %>%
  select(Year, d13C_wood, iWUE_simple, iWUE_photorespiration, iWUE_mesophyll, Ca, d13C.atm) %>%
  pivot_longer(cols = c(d13C_wood, iWUE_simple, iWUE_photorespiration, iWUE_mesophyll),
               names_to = "Variable",
               values_to = "Value")

Now I create a new dataset called data_subset by performing the following steps: Selecting columns from the combined_data dataset & remove missing values.


data_subset <- combined_data %>% 
  select(Year, wood.d13C, MGT_C, Elevation_m, Ca, iWUE_mesophyll, iWUE_photorespiration, iWUE_simple) %>%
  na.omit()

Linear regression models

Here I create the four linear regression models.


m1 <- lm(d13C.atm ~ Ca, data = CO2data) # Atmospheric CO2 concentration (ppm) & 13C sequestration
m1sum <- summary(m1)  

m2 <- lm(wood.d13C ~ iWUE_mesophyll + iWUE_photorespiration + iWUE_simple, data = data_subset) # 13C concentration & Intrinsic Water Use Effeciency
m2sum <- summary(m2)

m3 <- lm(Elevation_m ~ wood.d13C + iWUE_mesophyll + iWUE_photorespiration + iWUE_simple, data = data_subset) # Elevation & 13C concentration + iWUE fractions
m3sum <- summary(m3)

m4 <- lm(MGT_C ~ wood.d13C + iWUE_mesophyll + iWUE_photorespiration + iWUE_simple, data = data_subset) # Avg growth temperature & 13C concentration + iWUE fractions
m4sum <- summary(m4)

Statistical Analysis

Pearson Correlation Assessment: statistical measure that quantifies the linear relationship between two continuous variables (δ13CO2 / iWUE and environmental factors). It assesses the strength and direction of the association between the variables. Note that this test is sensitive to linear relationships but may not capture nonlinear relationships as it assumes that the relationship between the variables is linear and that the data follow a bivariate normal distribution.


corCO2data <- cor(CO2data[,2:3], method = "pearson")
cor_piruCO2 <- cor(data_subset[,2:8], method = "pearson")

Outlier Tests: Using the outlierTest function from the car package I detect outliers in the regression models. The Output is based on the studentized residuals of the model.


outm1 <- outlierTest(m1)
outm2 <- outlierTest(m2)
outm3 <- outlierTest(m3)
outm4 <- outlierTest(m4)

ANOVA: Now I perform an analysis of variance (ANOVA) to compare the means of two or more groups to determine if there are significant differences between them.


fit1 <- aov(m1)
fit2 <- aov(m2)
fit3 <- aov(m3)
fit4 <- aov(m4)

Non-constant variance Tests: Here I perform a NCV test on the regression models to to assess whether the assumption of constant variance (homoscedasticity) holds in the residuals of a regression model


ncv_m1 <- ncvTest(m1)

ncv_m2 <- ncvTest(m2)

ncv_m3 <- ncvTest(m3)

ncv_m4 <- ncvTest(m4)

Kruskal-Wallis test: The Kruskal-Wallis test is similar to an ANOVA as it is used to comare the medians of two or more independent groups.


kru_m1 <- kruskal.test(d13C.atm ~ Ca, data = CO2data)

kru_m2 <- kruskal.test(wood.d13C ~ iWUE_simple, data = data_subset)

kru_m3 <- kruskal.test(Elevation_m ~ iWUE_simple, data = data_subset)

kru_m4 <- kruskal.test(MGT_C ~ iWUE_simple, data = data_subset)

Results & Discussion

CR Plots

The CR plots provided here serve as illustrative representations of diagnostic plots utilized to evaluate the linearity assumption within linear regression models. These plots showcase the alignment of the fitted values (component) along the x-axis and the standardized residuals (residual) plotted along the y-axis.


crPlots(m1)


crPlots(m3)


crPlots(m4)


ggplot(data_long, aes(x = Year, y = Value, color = Variable)) +
  geom_point(alpha = 0.5) +
  geom_smooth(aes(group = Variable), color = "#000000") +
  theme_bw() +
  labs(x = "Year", y = "iWUE") +
  scale_x_continuous(breaks = seq(1940, 2020, 10), limits = c(1940, 2020)) +
  scale_color_tableau() +
  facet_wrap(~Variable, scales = "free_y") +
  theme(legend.position = "bottom")
`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'


print(m1sum)

Call:
lm(formula = d13C.atm ~ Ca, data = CO2data)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.165004 -0.033674 -0.001794  0.029147  0.215545 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.419e+00  1.707e-02  -83.13   <2e-16 ***
Ca          -1.789e-02  6.031e-05 -296.68   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.04166 on 2021 degrees of freedom
Multiple R-squared:  0.9776,    Adjusted R-squared:  0.9775 
F-statistic: 8.802e+04 on 1 and 2021 DF,  p-value: < 2.2e-16
print(m2sum)

Call:
lm(formula = wood.d13C ~ iWUE_mesophyll + iWUE_photorespiration + 
    iWUE_simple, data = data_subset)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.16693 -0.04550 -0.01113  0.03768  0.29708 

Coefficients:
                       Estimate Std. Error  t value Pr(>|t|)    
(Intercept)           -17.54772    0.15069 -116.450  < 2e-16 ***
iWUE_mesophyll         -3.12192    0.02338 -133.502  < 2e-16 ***
iWUE_photorespiration   1.59197    0.02226   71.517  < 2e-16 ***
iWUE_simple             0.09703    0.01453    6.677 1.98e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.07539 on 219 degrees of freedom
Multiple R-squared:  0.9933,    Adjusted R-squared:  0.9932 
F-statistic: 1.077e+04 on 3 and 219 DF,  p-value: < 2.2e-16
print(m3sum)

Call:
lm(formula = Elevation_m ~ wood.d13C + iWUE_mesophyll + iWUE_photorespiration + 
    iWUE_simple, data = data_subset)

Residuals:
    Min      1Q  Median      3Q     Max 
-96.087 -42.177  -8.935  37.646 130.304 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)           6994.3276   847.4852   8.253 1.47e-14 ***
wood.d13C              333.3020    47.9107   6.957 4.01e-11 ***
iWUE_mesophyll         884.0872   150.4895   5.875 1.57e-08 ***
iWUE_photorespiration -481.3544    77.8884  -6.180 3.11e-09 ***
iWUE_simple              0.3205    11.3045   0.028    0.977    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 53.45 on 218 degrees of freedom
Multiple R-squared:  0.5148,    Adjusted R-squared:  0.5059 
F-statistic: 57.82 on 4 and 218 DF,  p-value: < 2.2e-16
print(m4sum)

Call:
lm(formula = MGT_C ~ wood.d13C + iWUE_mesophyll + iWUE_photorespiration + 
    iWUE_simple, data = data_subset)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.059675 -0.003322  0.005101  0.008265  0.013859 

Coefficients:
                       Estimate Std. Error  t value Pr(>|t|)    
(Intercept)           -0.690341   0.208528   -3.311  0.00109 ** 
wood.d13C             -0.004443   0.011789   -0.377  0.70665    
iWUE_mesophyll         1.406583   0.037029   37.986  < 2e-16 ***
iWUE_photorespiration -2.722274   0.019165 -142.045  < 2e-16 ***
iWUE_simple            1.797626   0.002782  646.272  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.01315 on 218 degrees of freedom
Multiple R-squared:  0.9996,    Adjusted R-squared:  0.9996 
F-statistic: 1.356e+05 on 4 and 218 DF,  p-value: < 2.2e-16

Summary of Model 1: The analysis of Model 1 has revealed a highly significant correlation between atmospheric δ13C and atmospheric CO2 concentration (Ca) (p < 2.2e-16). The model effectively explains the variations in atmospheric δ13C, as evidenced by the adjusted R-squared value of 0.9775, indicating that approximately 97.75% of the fluctuation in atmospheric δ13C can be attributed to Ca concentration.

The estimated coefficient for Ca is -0.01789, with a standard error of 6.031e-05. This negative coefficient signifies that an elevation in Ca concentration corresponds to a decrease in atmospheric δ13C. Although the effect size is relatively modest, it remains statistically significant. The intercept term (-1.419) denotes the projected atmospheric δ13C when Ca concentration is zero.

The residuals of the model adhere to a normal distribution, as suggested by their descriptive statistics. The residual standard error of 0.04166 denotes the average magnitude of the discrepancies between observed and projected atmospheric δ13C values.

These outcomes are in line with previous research (Mathias et al., 2018; Adams et al., 2020; Guerrieria et al., 2020; Lavergne et al., 2022) that have corroborated a negative correlation between atmospheric CO2 concentration and carbon isotopes. The decline in atmospheric δ13C as Ca concentration increases can be attributed to plants’ propensity to favor 12C during photosynthesis. This selectivity leads to a reduction in the relative abundance of 13C within atmospheric CO2.

Summary of Model 2: The outcomes of the multiple linear regression analysis for wood δ13C data showcase significant results (p < 2.2e-16), implying a robust connection between wood δ13C and predictor variables: iWUE_mesophyll, iWUE_photorespiration, and iWUE_simple. The model boasts an adjusted R-squared value of 0.9932, indicating that approximately 99.32% of the variance in wood δ13C can be illuminated through the linear amalgamation of these variables.

The coefficient estimates for the predictor variables are as follows: iWUE_mesophyll (-3.12192), iWUE_photorespiration (1.59197), and iWUE_simple (0.09703). All coefficients are statistically significant (p < 0.001), signifying the profound impact of these variables on wood δ13C.

The intercept term (-17.54772) depicts the estimated wood δ13C when all predictor variables are zero. The negative coefficient for iWUE_mesophyll suggests that heightened mesophyll intrinsic water use efficiency corresponds to lower wood δ13C. Conversely, the positive coefficient for iWUE_photorespiration signifies that enhanced photorespiration intrinsic water use efficiency is linked to higher wood δ13C. Similarly, the positive coefficient for iWUE_simple indicates that greater simple intrinsic water use efficiency corresponds to higher wood δ13C.

The residual distribution adheres to normality, as demonstrated by descriptive statistics. The residual standard error of 0.07539 denotes the average magnitude of disparities between observed and projected wood δ13C values.

The high R-squared value underscores the model’s efficacy in explaining data variability, while the significant coefficients underscore the pivotal role of water use efficiency in shaping wood carbon isotopes. The negative correlation between iWUE_mesophyll and wood δ13C suggests that plants with heightened water use efficiency exhibit lower wood δ13C values, reflecting their 13C discrimination during photosynthesis. Conversely, the positive correlation between iWUE_photorespiration and wood δ13C implies that plants with enhanced photorespiration intrinsic water use efficiency tend to have higher wood δ13C values, potentially due to increased carbon loss during photorespiration.

Summary of Model 3: The multiple linear regression analysis applied to Elevation_m data unveils statistically significant results (p < 2.2e-16), underscoring the substantial connection between elevation and predictor variables: wood.d13C, iWUE_mesophyll, and iWUE_photorespiration. The adjusted R-squared value of 0.5059 indicates that around 50.59% of elevation variance can be attributed to these variables’ linear combination.

The coefficient estimates for the predictor variables are: wood.d13C (333.3020), iWUE_mesophyll (884.0872), and iWUE_photorespiration (-481.3544). All coefficients, excluding iWUE_simple, are statistically significant (p < 0.05), signifying their influence on elevation. The positive coefficient for wood.d13C implies that increased wood δ13C corresponds to higher elevation, potentially due to elevation-related environmental factors affecting wood carbon isotopes. Moreover, the positive coefficient for iWUE_mesophyll indicates that higher mesophyll intrinsic water use efficiency correlates with elevated elevations. Conversely, the negative coefficient for iWUE_photorespiration suggests that heightened photorespiration intrinsic water use efficiency is tied to lower elevations.

The residuals conform to normal distribution, as indicated by descriptive statistics. The residual standard error of 53.45 quantifies the average magnitude of discrepancies between observed and projected elevation values.

The moderate R-squared value highlights the model’s capability to elucidate a considerable portion of elevation variation, with space for additional unanalyzed factors. The significant coefficients underscore wood δ13C, mesophyll intrinsic water use efficiency, and photorespiration intrinsic water use efficiency as contributors to elevation variance. These findings intimate that physiological processes and carbon isotopic discrimination may intertwine with elevational gradients.

Summary of Model 4: The multiple linear regression analysis conducted on MGT_C data yields highly significant results (p < 2.2e-16), indicating a robust relationship between MGT_C and predictor variables: wood.d13C, iWUE_mesophyll, iWUE_photorespiration, and iWUE_simple.

The coefficient estimates for the predictor variables are as follows: wood.d13C (-0.004443), iWUE_mesophyll (1.406583), iWUE_photorespiration (-2.722274), and iWUE_simple (1.797626). All coefficients, excluding wood.d13C, are statistically significant (p < 0.05), emphasizing MGT_C’s impact on iWUE_mesophyll, iWUE_photorespiration, and iWUE_simple.

The positive coefficients for iWUE_mesophyll and iWUE_simple underscore their correlation with higher MGT_C values. Conversely, the negative coefficient for iWUE_photorespiration suggests an inverse relationship between photorespiration intrinsic water use efficiency and MGT_C values.

The residuals adhere to normality, as indicated by descriptive statistics. The residual standard error of 0.01315 captures the average magnitude of differences between observed and projected MGT_C values.

The high R-squared value signals the model’s adeptness at explaining a substantial portion of MGT_C variation, highlighting its relationship with predictor variables. The significant coefficients emphasize MGT_C’s role in shaping physiological plant processes.

These insights collectively underscore the intricate connections between atmospheric δ13C, wood δ13C, elevation, mean growth temperature (MGT_C), and water use efficiency variables. The models offer a comprehensive perspective on these interrelationships and their ecological implications. Further investigation is warranted to delve into underlying mechanisms and additional contributing factors.


print(corCO2data)
                 Ca   d13C.atm
Ca        1.0000000 -0.9887132
d13C.atm -0.9887132  1.0000000
print(cor_piruCO2)
                        wood.d13C         MGT_C   Elevation_m           Ca
wood.d13C              1.00000000 -1.920377e-01  6.499000e-01  0.022421975
MGT_C                 -0.19203768  1.000000e+00 -2.160729e-05  0.141279486
Elevation_m            0.64990003 -2.160729e-05  1.000000e+00 -0.005140573
Ca                     0.02242198  1.412795e-01 -5.140573e-03  1.000000000
iWUE_mesophyll         0.61323133 -2.532179e-02  3.666682e-01  0.801618096
iWUE_photorespiration  0.63546030 -3.305377e-02  3.808482e-01  0.784163081
iWUE_simple            0.63983670 -1.340447e-02  3.863895e-01  0.780950896
                      iWUE_mesophyll iWUE_photorespiration iWUE_simple
wood.d13C                 0.61323133            0.63546030  0.63983670
MGT_C                    -0.02532179           -0.03305377 -0.01340447
Elevation_m               0.36666816            0.38084821  0.38638952
Ca                        0.80161810            0.78416308  0.78095090
iWUE_mesophyll            1.00000000            0.99958960  0.99920154
iWUE_photorespiration     0.99958960            1.00000000  0.99975112
iWUE_simple               0.99920154            0.99975112  1.00000000

CO2 Data Correlation Summary: The correlation matrix for the CO2 data variables underscores a robust negative correlation (-0.9887) between atmospheric carbon dioxide concentration (Ca) and carbon isotope ratio (d13C.atm). This observation implies that as Ca levels increase, the d13C.isotope ratio diminishes. The correlation coefficient of -0.9887 underscores the powerful linear connection between these variables, signifying that fluctuations in atmospheric CO2 are intimately linked to changes in carbon isotopes. This alignment conforms to the established understanding that the combustion of fossil fuels, releasing carbon dioxide into the atmosphere, results in a lowered d13C.isotope ratio due to the isotopic constitution of fossil fuels.

piruCO2 Data Correlation Summary: The correlation matrix for the piruCO2 variables uncovers the subsequent correlations:

  1. A positive correlation between wood.d13C and Elevation_m (0.6499) implies that as elevation escalates, the wood carbon isotope ratio tends to ascend.
  2. A positive correlation between wood.d13C and Ca (0.0224) suggests a feeble connection between the wood carbon isotope ratio and atmospheric carbon dioxide concentration.
  3. Positive correlations between wood.d13C and intrinsic water-use efficiency measures, iWUE_mesophyll, iWUE_photorespiration, and iWUE_simple (0.6132, 0.6355, and 0.6398 respectively) convey that as the wood carbon isotope ratio increases, the intrinsic water-use efficiency measures also exhibit a tendency to rise.
  4. A negative correlation between MGT_C and wood.d13C (-0.192) indicates a subtle inverse relationship between the carbon isotope ratio and the mean growth temperature of the trees.

The remaining correlations among the variables are notably strong and positive, underlining a high degree of linkage. These correlations include MGT_C and Ca (0.1413), Elevation_m and intrinsic water-use efficiency measures, iWUE_mesophyll, iWUE_photorespiration, and iWUE_simple (0.3667, 0.3808, and 0.3864 respectively), and Ca and intrinsic water-use efficiency measures, iWUE_mesophyll, iWUE_photorespiration, and iWUE_simple (0.8016, 0.7842, and 0.7810 respectively). These observations imply that these variables are positively associated and have a tendency to increase in tandem.

In essence, the correlation matrix affords insights into the interconnections among the variables within the piruCO2 dataset. It reveals the strength and direction of the relationships between them, furnishing a comprehensive perspective on their associations.


print(outm1)
     rstudent unadjusted p-value Bonferroni p
2022 5.305261         1.2484e-07   0.00025255
2023 4.627304         3.9397e-06   0.00796990
2021 4.494333         7.3738e-06   0.01491700
print(outm2)
     rstudent unadjusted p-value Bonferroni p
2162  4.27554         2.8517e-05    0.0063593
print(outm3)
No Studentized residuals with Bonferroni p < 0.05
Largest |rstudent|:
     rstudent unadjusted p-value Bonferroni p
2160 2.504012           0.013015           NA
print(outm4)
      rstudent unadjusted p-value Bonferroni p
2024 -4.875862         2.0925e-06   0.00046663
2023 -4.811001         2.8095e-06   0.00062651
2022 -4.732518         3.9973e-06   0.00089141

Model 1 Outlier Test Interpretation: From the provided output, it’s evident that the studentized residuals for the years 2022, 2023, and 2021 stand out significantly from zero, indicating the potential presence of outliers. The unadjusted p-values and Bonferroni-adjusted p-values are available for each year, allowing for further analysis and decision-making on how to handle these potential outliers in the analysis.

Model 2 Outlier Test Interpretation: In Model 2, a specific observation with index number 2162 exhibits a studentized residual that significantly differs from zero, suggesting a potential outlier. The unadjusted p-value and Bonferroni-adjusted p-value are provided for this particular observation, offering insights into how to approach this potential outlier in the analysis.

Model 3 Outlier Test Interpretation: For Model 3, the unadjusted p-value associated with a residual is 0.013015, indicating a relatively moderate level of statistical significance. The Bonferroni-adjusted p-value is marked as NA, possibly due to the lack of other prominent outliers in the model.

Overall Interpretation of Outlier Tests for Model 3: The absence of significant outliers in Model 3 implies that the data points conform reasonably well to the assumptions of the linear regression model.

Model 4 Outlier Test Interpretation: In Model 4, certain observations, specifically those from the years 2024, 2023, and 2022, display substantial studentized residuals coupled with low unadjusted p-values. This suggests that these observations deviate significantly from the anticipated values as predicted by the regression model. The Bonferroni-adjusted p-values for these observations also being very small indicate strong statistical evidence against the null hypothesis of no outlier effect.

Overall Interpretation of Outlier Tests for Model 4: The outcomes imply that observations from 2024, 2023, and 2022 might be influential outliers, wielding notable influence over the predictive capabilities of the regression model for MGT_C.

In summary, these outlier test results furnish valuable information for evaluating the impact of potential outliers on the regression models. They aid in making informed decisions about whether to retain, adjust, or exclude these outliers in the subsequent analyses.


print(fit1)
Call:
   aov(formula = m1)

Terms:
                       Ca Residuals
Sum of Squares  152.75146   3.50743
Deg. of Freedom         1      2021

Residual standard error: 0.04165926
Estimated effects may be unbalanced
print(fit2)
Call:
   aov(formula = m2)

Terms:
                iWUE_mesophyll iWUE_photorespiration iWUE_simple Residuals
Sum of Squares        69.49156             113.80254     0.25335   1.24464
Deg. of Freedom              1                     1           1       219

Residual standard error: 0.07538766
Estimated effects may be unbalanced
print(fit3)
Call:
   aov(formula = m3)

Terms:
                wood.d13C iWUE_mesophyll iWUE_photorespiration iWUE_simple
Sum of Squares   542150.7         2089.6              116524.0         2.3
Deg. of Freedom         1              1                     1           1
                Residuals
Sum of Squares   622825.1
Deg. of Freedom       218

Residual standard error: 53.45087
Estimated effects may be unbalanced
print(fit4)
Call:
   aov(formula = m4)

Terms:
                wood.d13C iWUE_mesophyll iWUE_photorespiration iWUE_simple
Sum of Squares    3.46150        1.28552              16.83339    72.24427
Deg. of Freedom         1              1                     1           1
                Residuals
Sum of Squares    0.03771
Deg. of Freedom       218

Residual standard error: 0.01315185
Estimated effects may be unbalanced

Model 1 ANOVA Results Interpretation: The ANOVA table provides essential information about the distribution of sum of squares and degrees of freedom among the terms in Model 1. Specifically, the term “Ca” (atmospheric CO2 concentration) has a sum of squares of 152.75146 with one degree of freedom. The remaining sum of squares (3.50743) and degrees of freedom (2021) are allocated to the residuals, which signify the unexplained variability in the model. The residual standard error, offering an estimate of the standard deviation of residuals, is also presented. The note “Estimated effects may be unbalanced” indicates that the distribution of the atmospheric CO2 variable may not be even or uniform. This observation applies to the other models as well.

Model 2 ANOVA Results Interpretation: In Model 2, the ANOVA results show the sum of squares for the three Water Use Efficiency fractions: 69.49156, 113.80254, and 0.25335, each with one degree of freedom. The remaining sum of squares (1.24464) and degrees of freedom (219) correspond to the residuals. These figures provide insights into the proportion of variability attributed to each predictor variable and the residuals.

Model 3 ANOVA Results Interpretation: For Model 3, the ANOVA outcomes exhibit the sum of squares for the four Water Use Efficiency fractions: 542150.7, 2089.6, 116524.0, and 2.3, each with one degree of freedom. The remaining sum of squares (622825.1) and degrees of freedom (218) are attributed to the residuals. These statistics provide a breakdown of the variability apportioned to each predictor term and the residuals in the model.

Model 4 ANOVA Results Interpretation: Model 4’s ANOVA results showcase the sum of squares for the four Water Use Efficiency fractions: 3.46150, 1.28552, 16.83339, and 72.24427, respectively. Each term is associated with one degree of freedom. The residual sum of squares (0.03771) and degrees of freedom (218) account for the variability that remains unexplained by the model. These values grant insight into the relative contributions of each predictor variable and the residuals in the model’s overall variability.

In summary, the ANOVA results help dissect the distribution of variability within each regression model by breaking it down into the effects of predictor terms and residuals. These statistics facilitate understanding the significance and impact of each variable within the model.

print(ncv_m1)
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 777.4753, Df = 1, p = < 2.22e-16
print(ncv_m2)
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 45.17054, Df = 1, p = 1.806e-11
print(ncv_m3)
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 1.562031, Df = 1, p = 0.21137
print(ncv_m4)
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 48.59042, Df = 1, p = 3.1541e-12

Model 1 NCV Results Interpretation: For Model 1, the Non-Constant Variance (NCV) test indicates that the chi-square value is 777.4753 with 1 degree of freedom. The p-value, which is < 2.22e-16, suggests strong evidence against the assumption of constant variance. This implies that the variability of residuals is not consistent across all levels of the predictor variable. The non-constant variance can impact the reliability of the regression model’s predictions.

Model 2 NCV Results Interpretation: In Model 2, the NCV test yields a chi-square value of 45.17054 with 1 degree of freedom. The p-value of 1.806e-11 indicates strong evidence against the assumption of constant variance. Like in Model 1, this suggests that the variability of residuals is not uniform across different levels of the predictor variable.

Model 3 NCV Results Interpretation: For Model 3, the NCV test produces a chi-square value of 1.562031 with 1 degree of freedom. The p-value of 0.21137 indicates weak evidence against the assumption of constant variance. This suggests that the assumption of constant variance may hold reasonably well for Model 3, as the evidence against it is not strong.

Model 4 NCV Results Interpretation: In Model 4, the NCV test results in a chi-square value of 48.59042 with 1 degree of freedom. The p-value of 3.1541e-12 suggests strong evidence against the assumption of constant variance. This is similar to the findings in Models 1 and 2, indicating that the variability of residuals is not constant across different levels of the predictor variable.

Across all models, the non-constant variance test does not support the assumption of constant variance. This suggests that the variability of residuals changes as the predictor variables change. This violation of the assumption can lead to biased or inefficient parameter estimates and unreliable hypothesis tests. To account for the heteroscedasticity (non-constant variance) in the residuals, it may be advisable to explore alternative modeling approaches that can handle such conditions. Additionally, enhancing the models by incorporating different calculations and arguments could potentially yield improved results, considering the complexity of the relationships being investigated.


print(kru_m1)

    Kruskal-Wallis rank sum test

data:  d13C.atm by Ca
Kruskal-Wallis chi-squared = 1476.9, df = 754, p-value < 2.2e-16
print(kru_m2)

    Kruskal-Wallis rank sum test

data:  wood.d13C by iWUE_simple
Kruskal-Wallis chi-squared = 222, df = 222, p-value = 0.4874
print(kru_m3)

    Kruskal-Wallis rank sum test

data:  Elevation_m by iWUE_simple
Kruskal-Wallis chi-squared = 222, df = 222, p-value = 0.4874
print(kru_m4)

    Kruskal-Wallis rank sum test

data:  MGT_C by iWUE_simple
Kruskal-Wallis chi-squared = 222, df = 222, p-value = 0.4874

Model 1 Kruskal Test Results Interpretation: In Model 1, the Kruskal-Wallis test statistic is 1476.9, with 754 degrees of freedom. The p-value reported as “< 2.2e-16” indicates strong evidence against the null hypothesis of no difference in the medians. This suggests that there are significant differences in the medians of δ13C among the different levels of total atmospheric CO2 concentrations (Ca). Therefore, it can be concluded that the variable Ca has a statistically significant effect on the atmospheric carbon isotope composition (δ13C). Further investigations, such as post-hoc tests or pairwise comparisons, could help identify which specific groups have different median δ13C values.

Model 2-4 Kruskal Test Results Interpretation: For the remaining models (Models 2 to 4), the Kruskal-Wallis test results indicate that there is not sufficient evidence to conclude that there are significant differences in the medians of δ13C for the predictor variables being analyzed. This suggests that the variables related to iWUE fractions and environmental factors (such as Temperature and Elevation) may not have a statistically significant effect on δ13C values.

Overall, the Kruskal-Wallis test is a non-parametric method used to compare medians across different groups or levels of a categorical variable. In Model 1, the test results provide evidence that Ca significantly influences δ13C, while for Models 2 to 4, there isn’t enough evidence to support significant differences in δ13C related to the examined predictor variables. This insight can help guide further analysis and understanding of the relationships between these variables and δ13C in the context of the specific research or study.

Discussion & Conclusion

Throughout this study I sought to investigate the temporal trends in atmospheric carbon dioxide (CO2) concentrations (Ca) and carbon isotope composition (δ13C.atm) over the past 80 years, as well as the relationship between δ13C composition of woody organic matter and environmental factors. The analysis utilized data from Belmecheri and Lavergne (2020) and Mathias & Thomas (2021) and employed statistical modeling techniques to explore the dynamics and implications of these variables.

Our analysis revealed a consistent and significant increase in atmospheric CO2 concentrations over time, in line with previous studies. Since the industrial revolution, the average annual CO2 concentration has been steadily rising, reaching approximately 420 parts per million (ppm) in recent years. This upward trend in CO2 levels has important implications for δ13C values, as the preferential uptake of the lighter carbon isotope, 12C, by plants during photosynthesis results in lower δ13C values in response to increased CO2 concentrations.

The analysis of δ13C values over the study period confirmed the expected changes, demonstrating an increased proportion of carbon derived from atmospheric CO2 in plant tissues. This finding has potential implications for plant growth, nutrient uptake, and water-use efficiency, as it suggests altered carbon assimilation pathways and isotope discrimination patterns under elevated CO2 levels.

Moreover, our investigation into the relationship between δ13C composition of woody organic matter and environmental factors yielded valuable insights. The analysis revealed a highly significant negative relationship between atmospheric δ13C and Ca concentration. The negative coefficient for Ca suggests that as atmospheric CO2 concentrations increase, there is a corresponding decrease in δ13C values. This finding aligns with the principle of isotopic discrimination, as plants exhibit a higher affinity for 12C when CO2 availability is higher, leading to a decline in the relative abundance of 13C in atmospheric CO2.

Furthermore, the multiple linear regression analysis identified significant relationships between wood δ13C and water use efficiency (WUE) fractions, namely iWUE_mesophyll, iWUE_photorespiration, and iWUE_simple. The negative coefficient for iWUE_mesophyll indicates that higher mesophyll intrinsic water use efficiency is associated with lower wood δ13C values, reflecting the ability of plants to discriminate against 13C during photosynthesis. Conversely, the positive coefficient for iWUE_photorespiration suggests that higher photorespiration intrinsic water use efficiency is associated with higher wood δ13C values, likely due to increased carbon loss during photorespiration. The positive coefficient for iWUE_simple also indicates an increase in wood δ13C with increased simple intrinsic water use efficiency.

In summary, this study provides compelling evidence of the temporal trends in atmospheric CO2 concentrations and carbon isotope composition over the past 80 years. The findings underscore the impact of human activities on the increasing atmospheric CO2 levels and its consequent influence on the carbon isotope ratios in plant tissues. The negative relationship between atmospheric δ13C and CO2 concentrations highlights the role of isotopic discrimination in carbon assimilation processes. Additionally, the relationships between wood δ13C and WUE fractions highlight the intricate connections between plant physiological processes and carbon isotopic composition.

These findings contribute to our understanding of the changing dynamics of atmospheric CO2 and carbon isotope composition and provide valuable insights into the interactions between environmental factors and plant carbon assimilation. This knowledge is crucial for accurately predicting the responses of ecosystems to ongoing climate change and improving our ability to mitigate its impacts. Future research should continue to explore the complex relationships between atmospheric CO2, δ13C values, and environmental factors to further refine our understanding of these intricate processes and their implications for global carbon cycling.