Microbial Response to Rainfall – Gas Chromatography Analysis

---

R Notebook: Provides reproducible analysis for gas chromatography (GC) data in the following manuscript:

Citation: Romanowicz, KJ, Crump, BC, Kling, GW. (2021) Rainfall alters permafrost soil redox conditions, but meta-omics show divergent microbial community responses by tundra type in the arctic. Soil Systems 5(1): 17. https://doi.org/10.3390/soilsystems5010017

GitHub Repository: https://github.com/kromanowicz/2021-Romanowicz-SoilSystems

NCBI BioProject: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA666429

Accepted for Publication: Soil Systems 10 March 2021

Experiment

This R Notebook provides complete reproducibility of the data analysis in “Rainfall alters permafrost soil redox conditions, but meta-omics show divergent microbial community responses by tundra type in the Arctic” by Romanowicz, Crump, and Kling. In this experiment, mesocosms containing soil from the active layer of two dominant tundra types were subjected to simulated rainfall to alter redox conditions. The microbial functional potential (metagenomics) and gene expression (metatranscriptomics) patterns were measured during saturated anoxic redox conditions prior to rainfall and at multiple time points following the simulated rainfall event. Other measurements include soil properties as well as microbial respiration (CO₂) and methane (CH₄) production from soil subsamples collected at each sampling time point. The purpose was to determine if rainfall, as a form of soil oxidation, is sufficient to alter the anoxic redox conditions in arctic tundra and enhance the microbial degradation of organic carbon and CH₄ to CO₂.

Conceptual Figure. A total of 12 tundra mesocosms (3 replicates x 2 tundra types x 2 sets of response cores) were acclimated initially under anoxic redox conditions to mimic field conditions (T0). Dissolved oxygen was supplied to soils through the downward flow of oxygenated water during a simulated rainfall event. Dissolved oxygen will likely change the redox gradient directly following rainfall (T4) as a short-term effect. Anoxic conditions will likely be re-established after 24 hours (T24) as the pulse of oxygen is consumed through abiotic and biotic soil processes. Under anoxic redox conditions (T0), microorganisms likely degrade organic carbon through anaerobic and fermentation pathways, producing CH₄ and reducing Fe(III) to Fe(II). Rainfall-induced soil oxidation (T4) should stimulate heterotrophic microorganisms that degrade organic carbon and CH₄ through aerobic metabolic pathways, releasing CO₂. Soil oxidation should also stimulate aerobic autotrophic iron oxidizing bacteria that oxidize Fe(II) to Fe(III) and convert CO₂ into microbial biomass. The long-term response (T24) will likely be a combination of aerobic and anaerobic metabolism as well as a combination of reduction and oxidation iron reactions as dissolved oxygen is consumed. The predicted redox conditions and predicted redox reactions for coupled Fe(II)/Fe(III) cycling, as well as the microbial-induced release of CO₂ or CH₄ at each time point are based on the predicted availability of dissolved oxygen entering tundra soils through simulated rainfall.

Soil Sampling for Microbial Gene Expression

An initial soil sampling event for microbial activity was conducted at the end of the anoxic acclimation period (4-7 days) in all mesocosm replicates, representing sampling time point T0. Mesocosms were then flushed to simulate a rainfall event. Additional soil sampling events were conducted at T4 (4-hrs) and T24 (24-hrs) following the rainfall event to determine the temporal extent of microbial gene expression. Soil cores (2.54 cm diameter, 30 cm length) were extracted in duplicate from each mesocosm replicate at each sampling time point and homogenized by depth in 10-cm increments. The 10-20 cm soil increment, composed of organic soil in all mesocosm replicates, was chosen for microbial gene expression analysis and preserved in RNAlater Stabilization Reagent in sterile tubes at 4°C for 18 hours and then stored at -80°C until extraction.

Field Experiment. Tundra soil cores were collected from field sites in August 2017 (top left) and placed in buckets to establish the mesocosm experiment (bottom left). Tussock tundra cores were composed of an organic soil layer overlying a mineral soil layer (top middle) while wet sedge tundra cores were composed entirely of organic soil (bottom middle). Soil subsampling for microbial activity was taken from the 10-20 cm depth of duplicate soil cores in Tussock (top right) and Wet Sedge (bottom right).

# Make a vector of required packages
required.packages <- c("ape","cowplot","data.table","devtools","dplyr","DT","ggplot2","ggpubr","grid","gridExtra","kableExtra","knitr","pheatmap","png","RColorBrewer","reshape","rstatix","statmod","stringr","tibble","tidyr","tidyverse","vegan","yaml")

# Load required packages
lapply(required.packages, library, character.only = TRUE)

Soil Incubation

A subset of each soil sample collected from the 10-20 cm depth during the mesocosm experiment were sealed in jars for an incubation experiment to measure microbial respiration and methane production. Jars containing soil from the wet sedge tundra ecosystem were purged with N₂ to reproduce anaerobic field conditions. Each jar was sealed with an airtight lid that contained a septum for gas sampling. Jars were incubated for 5 days in the dark at field temperature (4\(^\circ\)C). Following the 5-day incubation, we removed 20 mL of headspace gas in an evacuated syringe and immediately analyzed CO₂ and CH₄ concentrations via gas chromotography.

Microbial Respiration

Plot the CO₂ data for microbial respiration.

# Place categories in preferred order for plotting
co2_gas$Veg<-factor(co2_gas$Veg, levels=c("Tussock","Wet Sedge"))

co2_gas$Sample <- factor(co2_gas$Sample,levels=c("Tuss-T0","WS-T0","Tuss-T4","WS-T4","Tuss-T24","WS-T24"))

co2_gas_barplot<-ggplot(co2_gas, aes(x = Veg, y = Mean, fill = Sample)) + geom_bar(stat = "identity", position = "dodge", color="black") + geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = 0.2, position = position_dodge(0.9)) + ylab(expression(atop("Microbial Respiration", paste(CO[2]~(ng~C~g^{"-1"}~soil~hr^{"-1"}))))) + theme_classic() + theme(axis.text=element_text(size=10), axis.title=element_text(size=12), axis.title.x=element_blank()) + theme(legend.position = "bottom", legend.title=element_blank(), legend.text=element_text(size=8), panel.background = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank()) + scale_size(guide=FALSE) + scale_fill_manual(values=c("palegreen","lightskyblue1", "green3", "dodgerblue1", "darkgreen", "midnightblue")) + scale_y_continuous(expand = c(0, 0), limits = c(0, 120)) + annotate(geom="text", x=1.0, y=85, label="italic(N.S.)", parse=TRUE) + annotate(geom="text", x=2.0, y=23, label="italic(N.S.)", parse=TRUE) + geom_segment(aes(x = 1.0, y = 90, xend = 1.0, yend = 100)) + geom_segment(aes(x = 1.0, y = 100, xend = 2.0, yend = 100)) + geom_segment(aes(x = 2.0, y = 100, xend = 2.0, yend = 32)) + annotate("text", x = 1.5, y = 110, label = "Mean Difference (4.5x)") + annotate("text", x = 1.5, y = 103.5, size=3, label = "Paired~t-test~(italic(p) < 0.001)", parse = TRUE)

co2_gas_barplot

Run statistics for differences between sampling time points within each tundra ecosystem.

# Load CO2 data
co2_anova<-read.csv("Stats.Data/co2.anova.csv")

# Run ANOVA for Tuss differences
tuss_co2_anova<-aov(Tuss_CO2~Time, data=co2_anova)
summary.aov(tuss_co2_anova)

            Df Sum Sq Mean Sq F value Pr(>F)
Time         2   21.4    10.7   0.033  0.968
Residuals    6 1941.3   323.6               

TukeyHSD(tuss_co2_anova)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = Tuss_CO2 ~ Time, data = co2_anova)

$Time
            diff       lwr      upr     p adj
T24-T0  3.687107 -41.37615 48.75037 0.9660478
T4-T0   2.553929 -42.50933 47.61719 0.9835136
T4-T24 -1.133178 -46.19644 43.93008 0.9967251

# Run ANOVA for WS differences
ws_co2_anova<-aov(WS_CO2~Time, data=co2_anova)
summary.aov(ws_co2_anova)

            Df Sum Sq Mean Sq F value Pr(>F)
Time         2  12.61   6.303   1.512  0.294
Residuals    6  25.02   4.170               

TukeyHSD(ws_co2_anova)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = WS_CO2 ~ Time, data = co2_anova)

$Time
             diff       lwr      upr     p adj
T24-T0  2.0835744 -3.032353 7.199502 0.4701724
T4-T0  -0.7038709 -5.819798 4.412056 0.9078917
T4-T24 -2.7874453 -7.903373 2.328482 0.2898464

Run t-test for differences in CO2 production between tundra ecosystems

t.test(co2.t.test$Tuss,co2.t.test$WS,paired=TRUE)

    Paired t-test

data:  co2.t.test$Tuss and co2.t.test$WS
t = 47.676, df = 2, p-value = 0.0004397
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 40.79162 48.88477
sample estimates:
mean of the differences 
               44.83819 

Methane Production

Plot the CH₄ data for methane production.

# Place categories in preferred order for plotting
ch4_gas$Veg<-factor(ch4_gas$Veg, levels=c("Tussock","Wet Sedge"))

ch4_gas$Sample <- factor(ch4_gas$Sample,levels=c("Tuss-T0","WS-T0","Tuss-T4","WS-T4","Tuss-T24","WS-T24"))

ch4_gas_barplot<-ggplot(ch4_gas, aes(x = Veg, y = Mean, fill = Sample)) + geom_bar(stat = "identity", position = "dodge", color="black") + geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = 0.2, position = position_dodge(0.9)) + ylab(expression(atop("Methane Production", paste(CH[4]~(pg~C~g^{"-1"}~soil~hr^{"-1"}))))) + theme_classic() + theme(axis.text=element_text(size=10), axis.title=element_text(size=12), axis.title.x=element_blank()) + theme(legend.position = "bottom", legend.title=element_blank(), legend.text=element_text(size=8), panel.background = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank()) + scale_size(guide=FALSE) + scale_fill_manual(values=c("palegreen", "lightskyblue1", "green3", "dodgerblue1", "darkgreen", "midnightblue")) + scale_y_continuous(expand = c(0, 0), limits = c(-10, 40)) + geom_hline(aes(yintercept = 0), color="black") + annotate(geom="text", x=0.7, y=25, label="a", parse=TRUE) + annotate(geom="text", x=1.0, y=-6, label="b", parse=TRUE) + annotate(geom="text", x=1.3, y=-6, label="b", parse=TRUE) + annotate(geom="text", x=1.7, y=13, label="a", parse=TRUE) + annotate(geom="text", x=2.0, y=-8, label="b", parse=TRUE) + annotate(geom="text", x=2.3, y=-8, label="b", parse=TRUE)

ch4_gas_barplot

Run statistics for differences between sampling time points within each tundra ecosystem.

# Run ANOVA for Tuss differences
tuss_ch4_anova<-aov(Tuss_CH4~Time, data=ch4_anova)
tuss_ch4_anova

Call:
   aov(formula = Tuss_CH4 ~ Time, data = ch4_anova)

Terms:
                    Time Residuals
Sum of Squares  758.3840   48.1623
Deg. of Freedom        2         6

Residual standard error: 2.833204
Estimated effects may be unbalanced

TukeyHSD(tuss_ch4_anova)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = Tuss_CH4 ~ Time, data = ch4_anova)

$Time
              diff        lwr        upr     p adj
T24-T0 -19.6798033 -26.777647 -12.581960 0.0003546
T4-T0  -19.2590856 -26.356929 -12.161242 0.0003999
T4-T24   0.4207178  -6.677126   7.518561 0.9819852

# Run ANOVA for WS differences
ws_ch4_anova<-aov(WS_CH4~Time, data=ch4_anova)
ws_ch4_anova

Call:
   aov(formula = WS_CH4 ~ Time, data = ch4_anova)

Terms:
                    Time Residuals
Sum of Squares  327.3446   17.3485
Deg. of Freedom        2         6

Residual standard error: 1.700418
Estimated effects may be unbalanced

TukeyHSD(ws_ch4_anova)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = WS_CH4 ~ Time, data = ch4_anova)

$Time
            diff        lwr       upr     p adj
T24-T0 -13.29290 -17.552844 -9.032951 0.0001826
T4-T0  -12.22735 -16.487295 -7.967402 0.0002923
T4-T24   1.06555  -3.194397  5.325496 0.7349851

Reproducibility

The session information is provided for full reproducibility.

devtools::session_info()

─ Session info ──────────────────────────────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.1.2 (2021-11-01)
 os       macOS Big Sur 10.16
 system   x86_64, darwin17.0
 ui       RStudio
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       America/Los_Angeles
 date     2022-01-24
 rstudio  1.4.1106 Tiger Daylily (desktop)
 pandoc   2.11.4 @ /Applications/RStudio.app/Contents/MacOS/pandoc/ (via rmarkdown)

─ Packages ──────────────────────────────────────────────────────────────────────────────────────────
 package      * version date (UTC) lib source
 abind          1.4-5   2016-07-21 [1] CRAN (R 4.1.0)
 ape          * 5.6-1   2022-01-07 [1] CRAN (R 4.1.2)
 assertthat     0.2.1   2019-03-21 [1] CRAN (R 4.1.0)
 backports      1.4.1   2021-12-13 [1] CRAN (R 4.1.0)
 brio           1.1.3   2021-11-30 [1] CRAN (R 4.1.0)
 broom          0.7.11  2022-01-03 [1] CRAN (R 4.1.2)
 cachem         1.0.6   2021-08-19 [1] CRAN (R 4.1.0)
 callr          3.7.0   2021-04-20 [1] CRAN (R 4.1.0)
 car            3.0-12  2021-11-06 [1] CRAN (R 4.1.0)
 carData        3.0-5   2022-01-06 [1] CRAN (R 4.1.2)
 cellranger     1.1.0   2016-07-27 [1] CRAN (R 4.1.0)
 cli            3.1.1   2022-01-20 [1] CRAN (R 4.1.2)
 cluster        2.1.2   2021-04-17 [1] CRAN (R 4.1.2)
 colorspace     2.0-2   2021-06-24 [1] CRAN (R 4.1.0)
 cowplot      * 1.1.1   2020-12-30 [1] CRAN (R 4.1.0)
 crayon         1.4.2   2021-10-29 [1] CRAN (R 4.1.0)
 data.table   * 1.14.2  2021-09-27 [1] CRAN (R 4.1.0)
 DBI            1.1.2   2021-12-20 [1] CRAN (R 4.1.0)
 dbplyr         2.1.1   2021-04-06 [1] CRAN (R 4.1.0)
 desc           1.4.0   2021-09-28 [1] CRAN (R 4.1.0)
 devtools     * 2.4.3   2021-11-30 [1] CRAN (R 4.1.0)
 digest         0.6.29  2021-12-01 [1] CRAN (R 4.1.0)
 dplyr        * 1.0.7   2021-06-18 [1] CRAN (R 4.1.0)
 DT           * 0.20    2021-11-15 [1] CRAN (R 4.1.0)
 ellipsis       0.3.2   2021-04-29 [1] CRAN (R 4.1.0)
 evaluate       0.14    2019-05-28 [1] CRAN (R 4.1.0)
 fansi          1.0.2   2022-01-14 [1] CRAN (R 4.1.2)
 farver         2.1.0   2021-02-28 [1] CRAN (R 4.1.0)
 fastmap        1.1.0   2021-01-25 [1] CRAN (R 4.1.0)
 forcats      * 0.5.1   2021-01-27 [1] CRAN (R 4.1.0)
 fs             1.5.2   2021-12-08 [1] CRAN (R 4.1.0)
 generics       0.1.1   2021-10-25 [1] CRAN (R 4.1.0)
 ggplot2      * 3.3.5   2021-06-25 [1] CRAN (R 4.1.0)
 ggpubr       * 0.4.0   2020-06-27 [1] CRAN (R 4.1.0)
 ggsignif       0.6.3   2021-09-09 [1] CRAN (R 4.1.0)
 glue           1.6.0   2021-12-17 [1] CRAN (R 4.1.0)
 gridExtra    * 2.3     2017-09-09 [1] CRAN (R 4.1.0)
 gtable         0.3.0   2019-03-25 [1] CRAN (R 4.1.0)
 haven          2.4.3   2021-08-04 [1] CRAN (R 4.1.0)
 highr          0.9     2021-04-16 [1] CRAN (R 4.1.0)
 hms            1.1.1   2021-09-26 [1] CRAN (R 4.1.0)
 htmltools      0.5.2   2021-08-25 [1] CRAN (R 4.1.0)
 htmlwidgets    1.5.4   2021-09-08 [1] CRAN (R 4.1.0)
 httr           1.4.2   2020-07-20 [1] CRAN (R 4.1.0)
 jquerylib      0.1.4   2021-04-26 [1] CRAN (R 4.1.0)
 jsonlite       1.7.3   2022-01-17 [1] CRAN (R 4.1.2)
 kableExtra   * 1.3.4   2021-02-20 [1] CRAN (R 4.1.0)
 knitr        * 1.37    2021-12-16 [1] CRAN (R 4.1.0)
 labeling       0.4.2   2020-10-20 [1] CRAN (R 4.1.0)
 lattice      * 0.20-45 2021-09-22 [1] CRAN (R 4.1.2)
 lifecycle      1.0.1   2021-09-24 [1] CRAN (R 4.1.0)
 lubridate      1.8.0   2021-10-07 [1] CRAN (R 4.1.0)
 magrittr       2.0.1   2020-11-17 [1] CRAN (R 4.1.0)
 MASS           7.3-55  2022-01-13 [1] CRAN (R 4.1.2)
 Matrix         1.4-0   2021-12-08 [1] CRAN (R 4.1.0)
 memoise        2.0.1   2021-11-26 [1] CRAN (R 4.1.0)
 mgcv           1.8-38  2021-10-06 [1] CRAN (R 4.1.2)
 modelr         0.1.8   2020-05-19 [1] CRAN (R 4.1.0)
 munsell        0.5.0   2018-06-12 [1] CRAN (R 4.1.0)
 nlme           3.1-155 2022-01-13 [1] CRAN (R 4.1.2)
 permute      * 0.9-5   2019-03-12 [1] CRAN (R 4.1.0)
 pheatmap     * 1.0.12  2019-01-04 [1] CRAN (R 4.1.0)
 pillar         1.6.4   2021-10-18 [1] CRAN (R 4.1.0)
 pkgbuild       1.3.1   2021-12-20 [1] CRAN (R 4.1.0)
 pkgconfig      2.0.3   2019-09-22 [1] CRAN (R 4.1.0)
 pkgload        1.2.4   2021-11-30 [1] CRAN (R 4.1.0)
 plyr           1.8.6   2020-03-03 [1] CRAN (R 4.1.0)
 png          * 0.1-7   2013-12-03 [1] CRAN (R 4.1.0)
 prettyunits    1.1.1   2020-01-24 [1] CRAN (R 4.1.0)
 processx       3.5.2   2021-04-30 [1] CRAN (R 4.1.0)
 ps             1.6.0   2021-02-28 [1] CRAN (R 4.1.0)
 purrr        * 0.3.4   2020-04-17 [1] CRAN (R 4.1.0)
 R6             2.5.1   2021-08-19 [1] CRAN (R 4.1.0)
 RColorBrewer * 1.1-2   2014-12-07 [1] CRAN (R 4.1.0)
 Rcpp           1.0.8   2022-01-13 [1] CRAN (R 4.1.2)
 readr        * 2.1.1   2021-11-30 [1] CRAN (R 4.1.0)
 readxl         1.3.1   2019-03-13 [1] CRAN (R 4.1.0)
 remotes        2.4.2   2021-11-30 [1] CRAN (R 4.1.0)
 reprex         2.0.1   2021-08-05 [1] CRAN (R 4.1.0)
 reshape      * 0.8.8   2018-10-23 [1] CRAN (R 4.1.0)
 rlang          0.4.12  2021-10-18 [1] CRAN (R 4.1.0)
 rmarkdown      2.11    2021-09-14 [1] CRAN (R 4.1.0)
 rprojroot      2.0.2   2020-11-15 [1] CRAN (R 4.1.0)
 rsconnect      0.8.25  2021-11-19 [1] CRAN (R 4.1.0)
 rstatix      * 0.7.0   2021-02-13 [1] CRAN (R 4.1.0)
 rstudioapi     0.13    2020-11-12 [1] CRAN (R 4.1.0)
 rvest          1.0.2   2021-10-16 [1] CRAN (R 4.1.0)
 scales         1.1.1   2020-05-11 [1] CRAN (R 4.1.0)
 sessioninfo    1.2.2   2021-12-06 [1] CRAN (R 4.1.0)
 statmod      * 1.4.36  2021-05-10 [1] CRAN (R 4.1.0)
 stringi        1.7.6   2021-11-29 [1] CRAN (R 4.1.0)
 stringr      * 1.4.0   2019-02-10 [1] CRAN (R 4.1.0)
 svglite        2.0.0   2021-02-20 [1] CRAN (R 4.1.0)
 systemfonts    1.0.3   2021-10-13 [1] CRAN (R 4.1.2)
 testthat       3.1.2   2022-01-20 [1] CRAN (R 4.1.2)
 tibble       * 3.1.6   2021-11-07 [1] CRAN (R 4.1.0)
 tidyr        * 1.1.4   2021-09-27 [1] CRAN (R 4.1.0)
 tidyselect     1.1.1   2021-04-30 [1] CRAN (R 4.1.0)
 tidyverse    * 1.3.1   2021-04-15 [1] CRAN (R 4.1.0)
 tinytex        0.36    2021-12-19 [1] CRAN (R 4.1.0)
 tzdb           0.2.0   2021-10-27 [1] CRAN (R 4.1.0)
 usethis      * 2.1.5   2021-12-09 [1] CRAN (R 4.1.0)
 utf8           1.2.2   2021-07-24 [1] CRAN (R 4.1.0)
 vctrs          0.3.8   2021-04-29 [1] CRAN (R 4.1.0)
 vegan        * 2.5-7   2020-11-28 [1] CRAN (R 4.1.0)
 viridisLite    0.4.0   2021-04-13 [1] CRAN (R 4.1.0)
 webshot        0.5.2   2019-11-22 [1] CRAN (R 4.1.0)
 withr          2.4.3   2021-11-30 [1] CRAN (R 4.1.0)
 xfun           0.29    2021-12-14 [1] CRAN (R 4.1.0)
 xml2           1.3.3   2021-11-30 [1] CRAN (R 4.1.0)
 yaml         * 2.2.1   2020-02-01 [1] CRAN (R 4.1.0)

 [1] /Library/Frameworks/R.framework/Versions/4.1/Resources/library

─────────────────────────────────────────────────────────────────────────────────────────────────────