Data Package Metadata   View Summary

Site selection

Wildfires have been suppressed in GSMNP since 1934 following a national policy of fire suppression within national parks and, within recent decades, concern about the increasing wildland-urban interface. During late November and early December 2016, the Chimney Tops 2 wildfire burned approximately 44.4 km² within the GSMNP. On November 28^th, the fire moved into the city of Gatlinburg, TN and burned a further 28.2 km². This wildfire was extremely heterogeneous, with unburned sites adjacent to severely burned sites, often within meters of each other. From within the burn matrix we randomly selected three sites from “high burn” areas, three “low/medium burn” sites, and three “no burn” sites both within natural areas of the GSMNP and the exurban area of Gatlinburg, TN. Burn levels were determined using a Burned Area Reflectance Classification (BARC) map generated by US Forest Service Remote Sensing Application Center from 2015 and 2016 EO ALI and Landsat satellite images. In brief, the map was generated by calculating a normalized burn ratio from two spectral bands, near infrared and shortwave infrared, pre- and post-fire, and represents the soil burn severity. Burned categories typically fall within 3 classes of fire severity, with Unburned/Very Low = <75; Medium = 76–187; and High = > 187. Satellite imagery data was collected from the Forest Service Remote Sensing Applications Center (RSAC) and the U.S. Geological Survey Center for Earth Resources Observation and Science (EROS). While BARC classes matched visual estimates of burn severity at our sites (e.g., amount of bark scorch, plant mortality, soil organic matter and litter layer accumulation; Hughes et al. 2020), satellite data can only resolve burn severity at coarser scales (30 m × 30 m) than our 1 m × 1 m sampling plots. Therefore, these estimates should only be interpreted as approximations of burn severity. More details on soil responses to fire are presented below.

Over-story plant community at all sites is primarily mixed deciduous hardwood with oak species dominant. Elevation ranges from 407 to 630m asl. Precipitation varies between mountains and valleys with an annual mean between 1397–2159 mm and mean daily temperatures range from 10 °C in January to 31 °C in July at low elevations and 1 °C to 18 °C at high elevations (Shanks 1954).

Sampling

In September 2018 we collected fine roots from one individual of each of the five dominant understory species co-occurring in established two 1 m × 1 m plots at each of 18 sites. All plant species were growing with active AM fungi within plant roots. Plant surveys occurred in these sites every four months following the Chimney Tops 2 fire; the sites showed minimal differences in composition or diversity among fire severities or location (exurban versus natural) in 2018 (Hubert et al. in prep.). We only sampled plant host species that were present in at least two burn severity categories and both exurban and natural settings to minimize the confounding effects of plant host and these other parameters (Supplemental Table 1). Species included Acer saccharum, Amphicarpaea bracteata, Kalmia latifolia, Lactuca canadensis, Lysimachia quadrifolia, Panicum milliaceum, Parthenocissus quinquefolia, Phytolacca americana, Smilax glauca, and Viola palmata. Species ranged from annual (Amphicarpaea, Lactuca, Panicum) to perennial (Acer, Kalmia, Lysimachia, Parthenocissus, Phytolacca, Smilax, Viola) life cycles and woody tree and shrub saplings (Acer, Kalmia), woody vines (Parthenocissus, Smilax), herbaceous forbs (Amphicarpaea, Lactuca, Lysimachia, Phytolacca, Viola), and grasses (Panicum).

In the field, fine root samples (< 2 mm in diameter) were immediately placed in sterile whirlpak bags and placed on ice until transport to the laboratory (within 24 hours). Once in the laboratory, roots were surface-sterilized with 70% EtOH and 10% NaClO followed by a sterile water rinse and stored frozen at -80°C.

In October 2018 we also collected bulk soils from all 18 sites. At each site we collected the organic and mineral horizons from ten 5 cm × 10 cm cores. Cores were homogenized into one aggregate sample for each horizon at each of the 18 sites. Following fire, only mineral horizon soils remained in high burn sites. From these homogenized samples we measured soil organic matter content via mass loss on ignition at 360°C for 3 hours (Davies 1974) and soil nutrients following extraction in 0.5 M K₂SO₄ (ammonium and nitrate) or 0.5M HCl (phosphate). All soil nutrients were quantified colorimetrically (see methods in D’Angelo et al. 2001, Doane and Horwath 2003).

Microscopy

Plant roots from each sample were cleared in 10% KOH and stained with 0.01% acid fuchsin. Roots were scored for percent colonization of AM fungal hyphae, arbuscules, and vesicles with the gridline intercept method (McGonigle et al. 1990) on a Laxco 3000 microscope at 200x magnification.

DNA Extraction and amplification

Roots were ground to a fine powder using a sterile mortar and pestle with liquid nitrogen. DNA was then extracted from approximately 250 mg of roots using the Qiagen DNeasy Plant mini kit (Qiagen, Germantown, MD, USA). DNA was quantified with the Qubit high sensitivity kit (Qubit Fluorometer, Life Technologies, Carlsbad, CA, USA) and diluted to ~10 ng/µl in sterile water. Due to limited AM fungal DNA, we then performed a nested PCR reaction. The first reaction amplified a ~800 b region of AM fungi and plants in the 18S region using the NS1 – NS4 primers (White et al. 1990), the preferred marker gene for AM fungi (Lekberg et al. 2018). The nested reaction amplified a ~400 b region of 18S AM fungal DNA with barcoded Illumina TruSeq V3 indices (Illumina, San Diego, CA, USA) linked to the NS31 -AML2 primers (Morgan and Egerton-Warburton 2017). Each reaction contained: 21.5 µl of Platinum PCR Supermix (Invitrogen, Carlsbad, CA, USA), 1.25 µl of each primer (10 µM), 0.5 µl of BSA (20 mg/ml), and 2 µl (~20 ng) of DNA. The first reaction ran at 94 °C for 3 min, followed by 30 cycles of 94 °C for 30 sec, 40 °C for 1 min, and 72 °C for 1 min and the nested reaction at 94 °C for 5 min, followed by 40 cycles of 94 °C for 45 sec, 63.1 °C for 1 min, and 72 °C for 1.5 min. Triplicate nested reactions were combined, cleaned with Agencourt AMPure XP magnetic beads (Beckman Coulter, Brea, CA, USA), and quantitated fluorometrically (Qubit Fluorometer, Life Technologies, Carlsbad CA, USA). Samples were pooled into equal amounts and run on an Illumina MiSeq v3 sequencer in a 2 × 275 b run at the University of Tennessee Center for Environmental Biotechnology Core Facility.

Bioinformatics

AM fungal sequences were processed in the DADA2 pipeline in R (Callahan et al. 2016). First, primers were trimmed from all sequences and sequence error rates were calculated. Sequences were then merged into unique amplicon sequence variants (ASVs). Finally, chimeras were removed using a denovo chimera checker. Because the NS31-AML2 primers may amplify some non-AM fungal fungi, we then BLASTed representative sequence reads from each ASV against the MaarjAM database (Opik et al. 2010) and only retained reads that matched a known AM fungal virtual taxonomic unit by at least 97%. Sequences are deposited in the NCBI Sequence Read Archive (BioProjectID: PRJNA771625). All other data is available upon request from the authors.

Statistics

All statistics were performed in R v. 3.5.1 (R Core Team, 2019). Soil nutrients (n= 26) were analyzed using a general linear model with the fixed effects of soil horizon, burn severity, and urbanization status (i.e., exurban sites in Gatlinburg or natural sites from GSMNP), and their interactions, and the random effect of site using the lmer function in the lme4 package. When soil horizons differed in nutrient concentrations (NH₄ ⁺ and NO₃ ^-) we also performed an ANOVA on each soil horizon separately with the fixed effects of burn severity, and urbanization and their interaction using the aov function in base R.

AM fungal colonization rates (n = 100) were logit transformed to improve normality before analysis. We used general linear mixed effects models to test how AM fungal colonization rates (hyphae, arbuscules and vesicles) and alpha diversity (richness, Shannon’s Diversity, Simpson’s Diversity) varied with the fixed effects of plant host identity, burn severity, and urbanization status with the random effect of site using the lmer function in the lme4 package (Bates et al. 2015). We tested for all pairwise interactions but did not have enough statistical power to test for the three-way interaction of plant host identity × burn severity × urbanization status. When significant, we tested for post-hoc differences within fixed effects while correcting for multiple comparisons with the false discovery rate of α = 0.05 in the lsmeans package (Lenth 2016).

AM fungal composition (n = 118) was visualized in non-dimensional multivariate space using the Vegan package (Oksanen et al. 2009). Variation in AM fungal composition measured as beta-diversity (turnover among individual sites) was partitioned into fractions owing to plant host identity, burn severity, space, and urbanization using the VarPart function in Vegan (Oksanen et al. 2009). Since urbanization status did not explain any of the variation in AM fungal composition, it was dropped from the analysis. Overlap in AM fungal taxa among burn severity classes and in exurban versus natural samples were visualized using Venn diagrams with the VennDiagram package in R (Chen and Boutros 2011). Finally, heterogeneity in AM fungal composition was tested for each factor (i.e., burn severity, urbanization status, plant host identity) with the betadispr function in Vegan (Oksanen et al. 2009).

References

Bates D, M Machler, B Bolker, and S Walker. 2015. Fitting linear mixed-effects models using lme4. Journal of Statistical Software 67:1-48.

Callahan BJ, PJ McMurdie, MJ Rosen, AW Han, AJA Johnson, and SP Holmes. 2016. DADA2:High-resolution sample inference from Illumina amplicon data. Nature Methods 13:581-583.

Chen H and PC Boutros. 2011. VennDiagram: a package for the generation of highly-customizable Venn and Euler diagrams in R. BMC Bioinformatics 12:35.

D’Angelo E, J Crutchfield, and M Vandiviere. 2001. Rapid, sensitive, microscale determination of phosphate in water and soil. Journal of Environmental Quality 30:2206-2209.

Davies BE. 1974. Loss-on-ignition as an estimate of soil organic matter. Soil Science Society of America Journal 38:150-151.

Doane TA and WR Horwath. 2003. Spectrophotometric determination of nitrate with a single reagent. Analytical Letters 36:2713-2722.

Lekberg Y, M Vasar, LS Bullington, S-K Sepp, PM Antunes, R Bunn, BG Larkin, and M Opik. 2018. More bang for the buck? Can arbuscular mycorrhizal fungal communities be characterized adequately alongside other fungi using general fungal primers? New Phytologist 220:971-976.

Lenth RV. 2016. Least-square means: The R package lsmeans. Journal of Statistical Software 69:1-33.

McGonigle TP, MH Miller, DG Evans, GL Fairchild, and JA Swan. 1990. A new method which gives an objective measure of colonization of roots by vesicular arbuscular mycorrhizal fungi. New Phytologist 115:495-501.

Morgan BST and LM Egerton-Warburton. 2017. Barcoded NS31/AML2 primers for sequencing arbuscular mycorrhizal communities in environmental samples. Applied Plant Sciences 5:1700017.

Oksanen J, R Kindt, P Legendre, B O’Hara, GL Simpson, P Solymos et al. 2009. Vegan: community ecology package. R package version 1.15-4.

Opik M, A Vanatoa, E Vanatoa, M Moora, J Davison, JM Kalwij, U Reier, and M Zobel. 2010. The online database MaarjAM reveals global and ecosystemic distribution patterns in arbuscular mycorrhizal fungi (Glomeromycota). New Phytologist 188:223-241.

R Development Core Team. 2009. R: A language and environment for statistical computing. R foundation for statistical computing. Vienna, Austria. http://www.R-project.org/

Shanks RE. 1954. Climates of the Great Smoky Mountains. Ecology 354-361.

White TJ, TD Bruns, SH Lee, and JW Taylor. 1990. PCR protocols: a guide to methods and application. San Diego pg 315-320.