Effects of biodegradable plastic film mulching on soil microbial communities in two agroecosystems

Plastic mulch films

Three commercially available biodegradable mulch films (BioAgri®, Naturecycle, and Organix A.G. Film™) and one experimental film comprised of a blend of polylactic acid (PLA) and polyhydroxyalkanoates (PHA) were tested alongside a polyethylene (PE) mulch (negative control), and cellulose paper mulch (WeedGuard Plus®, positive control). The paper mulch used in the experiment was a 100% biobased product and was chosen as the positive control because its major constituent is cellulose which is known to rapidly disintegrate in the field (Li et al., 2014b). Physicochemical properties of mulches are reported in Table 1.

Field trial description

Field experiments were set up in two locations: East Tennessee Research and Education Center (ETREC), University of Tennessee, Knoxville, TN and the Northwestern Washington Research & Extension Center (NWREC), Washington State University, Mount Vernon, WA. The soil at Knoxville is a sandy loam (59.9% sand, 23.5% silt, and 16.6% clay), classified as a fine kaolinitic thermic Typic Paleudults. The soil at Mount Vernon is a silt loam (14.2% sand, 69.8% silt, and 16% clay), classified as a fine-silty mixed nonacid mesic Typic Fluvaquents. Henceforth in the paper, Knoxville, TN will be referred to as TN and Mount Vernon, WA will be referred to as WA. The mulches were tested in the field over two years (2015 to 2016) under pie pumpkin (Cucurbita pepo) as a test crop, with full experimental details described in Sintim et al. (2019) and Ghimire et al. (2018b). Briefly, before mulch application began in TN and WA, the plots were under winter wheat (Triticum aestivum) cover crop in TN and clover (Trifolium spp.) at WA. Each field site was arranged as a randomized complete block design with four replications of seven main plot treatments (six mulch treatments described above and one no mulch control). Mulches were machine-laid on raised beds at the end of May to early June, and harvest was completed in September-October. PE mulch was removed soon after harvest, while BDMs were tilled in; all beds were tilled within two weeks of harvest. The two sites were planted with a winter wheat cover crop following harvest in the fall. Pie pumpkin (Cucurbita pepo) was used as the test crop because it met requirements of a large-scale, multi-location field experiment; namely, pie pumpkin is economically important, commonly grown in both our field experimental locations, and has sufficient season length to maximize treatment exposure. Additionally, pie pumpkins, like other cucurbits such as cucumbers, melons and squash, are commonly grown on plastic mulch films throughout the United States (Inglis, Miles & Wszelaki, 2015).

Soil water content and temperature were monitored as described in Sintim et al. (2019). Briefly, sensors (5TM, Decagon Devices Inc., Pullman, WA) installed in the center of each mulch treatment at 10-cm and 20-cm soil depths for one field block were connected to data loggers (EM50G, Decagon Devices Inc., Pullman, WA) that recorded the soil water and temperature data hourly. Soil water content and temperature data is reported in Sintim et al. (2019). Air temperature, precipitation, relative humidity, wind, and solar radiation were collected from a meteorological station located at the field site at TN (Decagon Devices Inc. Weather Station, Pullman, WA), and about 100 m away from the field site at WA (WSU AgWeatherNet Station, Mount Vernon, WA). Weather data for the two locations were continuously collected from 2015–2017.

Soil physical, chemical, and biological properties were assessed over the two-year study for this site, in order to assess changes in soil health. Detailed protocols for these measurements and raw data is provided in Sintim et al. (2019).

Soil sampling

Soil samples were collected from each of the 28 plots (seven treatments, replicated four times) at both locations in the Spring (May) and Fall (September) of 2015 and 2016. Spring samples were collected approximately 2 weeks prior to mulch application. Fall soil samples were collected while the mulches were still in the field, approximately 2 to 3 weeks before the mulch was tilled in (BDMs) or removed (PE). In our studies, BDMs did not full degrade over the winter, so the 2016 samples represent communities that have been exposed to (and are presumably still degrading) tilled-in mulch from the previous season; Spring 2016 soil had been exposed to tilled-in mulches from the previous season, and Fall 2016 soil had been exposed to both tilled-in mulches from the previous season and new mulch laid for the 2016 season. Soil was collected from the top 10 cm using a 2 cm diameter stainless steel auger. Thirty 10-cm soil cores were taken about 20 cm apart and composited for each of the plots. All sampling equipment was cleaned with 70% ethanol between plots to limit cross contamination. Roots and pebbles were removed by hand, and soils homogenized and stored in plastic bags for transport back to the lab. Soils were stored at −80 °C until DNA extraction and extracellular enzyme assays.

Soil DNA extraction and quantification

Extraction of DNA from soil samples was completed using the MoBio™ PowerLyzer™ Power Soil DNA isolation kit (now branded under Qiagen™) with inhibitor removal technology, as per manufacturer’s instructions. 0.25 grams of soil were used for the extractions, and the DNA obtained after the final elution step was stored at −20 °C until further analyses.

Quantification of the DNA extracted from soil was completed using the Quant-It™ PicoGreen™ dsDNA Quantification Kit (ThermoFisher Scientific) per manufacturer’s instructions and quality of DNA was measured by 260/280 ratios in a NanoDrop™ Spectrophotometer (Table S1).

Quantitative PCR for bacterial and fungal abundances

As a proxy for bacterial and fungal abundances, 16S rRNA (bacteria) and ITS (fungi) gene copy abundances were quantified from soil DNA samples using Femto™ Bacterial DNA quantification kit (Zymo Research) and Femto™ Fungal DNA quantification kit (Zymo Research) following the manufacturer’s protocol. DNA extracts were diluted 1:10 prior to quantification and 1 µl of the diluted samples was used for each qPCR reaction. All samples were analyzed in triplicate. No-template negative controls were included in each run. Bacterial and fungal DNA standards were provided in the kit and the ng DNA standard per well was converted to copy numbers which were used for final calculations. qPCR reactions were performed in a CFX Connect Real-Time PCR Detection System (BioRad). qPCR efficiencies averaged around 85% and 90% for bacterial and fungal assays, respectively.

DNA amplification and sequencing

16S rRNA amplicon sequencing of DNA extracts was conducted by the Genomic Services Laboratory (GSL) at Hudson Alpha, Huntsville, AL, following their standard operating procedures. Extracted DNA samples were shipped frozen in 96 well plates. The V4 region of the 16S rRNA gene was amplified using primers 515F (GTGCCAAGCAGCCGCGGTAA) and 806R (GGACTACHVGGGTWTCTAAT) (Caporaso et al., 2012). The first PCR used V4 amplicon primers, Kapa HiFi master mix, and 20 cycles of PCR. All aliquots and dilutions of the samples were completed using the Biomek liquid handler. PCR products were purified and stored at −20 °C until further processing was completed. PCR indexing was completed for 16S (V4) amplicons using GSL3.7/PE1 primers, Kapa HiFi master mix, and 12 cycles of PCR. Products were purified using magnetic beads using the Biomek liquid handler. Final libraries were quantified using Pico Green. V4 amplicon size obtained was 425 bp for the soil samples. The amplified 16S rRNA genes were sequenced using 250 paired-end reads on an Illumina MiSeq platform. Sequence reads were deposited in the NCBI sequence read archive (Accession PRJNA564156).

Raw sequence data was processed using mothur v.1.39.5 following the MiSeq SOP (Schloss et al., 2009). Before aligning to the reference database (SILVA release 102), unique sequences were identified, and a count table generated. After alignment to SILVA database, sequences were filtered to remove overhangs at both ends, and sequences de-noised by pre-clustering sequences with up to two nucleotide differences. Chimeras were removed using the VSEARCH algorithm. Sequences were classified using the Bayesian classifier (Wang et al., 2007) against the Ribosomal Database Project training set (version 9) with a bootstrap value of >80% (Wang et al., 2007). Following this step, untargeted (i.e., non-bacterial) sequences classified as Eukaryota and Arachaeota were removed. Sequences were binned into phylotypes according to their taxonomic classification at the genus level. A consensus taxonomy for each OTU was generated by comparison to the RDP training set. The resulting OTU count table and taxonomy assignments were imported into R (v. 3.4.0) (R Core Team, 2018) for further downstream statistical analyses. Mothur code, R code and associated input files are available at: https://github.com/jdebruyn/BDM-Microbiology.

Extracellular enzyme assays

Fluorescence microplate enzyme assays were conducted using fluorescently labelled substrates to assess enzyme activities in soil (Bell et al., 2013). Seven enzymes were assayed using their respective fluorescent substrates and standards (Table S2).

Soil slurries were prepared in a sodium acetate trihydrate buffer whose pH was matched closely with the soil pH. 800 µl of soil slurry was pipetted into deep well 96 well plates. Separate plates were prepared for 4-methylumbelliferone (MUB) and 7-amino-4-methylcoumarin (MUC) standard curves for each sample. 200 µl of appropriate standards and substrates were added to the soil slurries. The plates were sealed and inverted to mix the contents. Incubation was done for 3 h at room temperature, after which the substrate and standard plates were centrifuged at 1,500 rpm (∼327 × g) for 3 min. The supernatants were pipetted into black 96 well plates and fluorescence measured at 365 nm excitation wavelength and 450 nm emission wavelength in a BioTek® Synergy plate reader.

Statistical analyses

Beta diversity was computed using Bray-Curtis distances of microbial community composition using the vegan package (v 2.4-3) in R version 3.4.0 (R Core Team, 2018) based on OTU tables, and were then visualized using non-metric multidimensional scaling (NMDS) using phyloseq package v.1.21.0 in R (McMurdie & Holmes, 2013). To determine whether significant differences existed in bacterial community composition between bacterial communities across different locations, seasons, and mulch treatments, a permutational multivariate analysis of variance (PERMANOVA) was performed using the ADONIS function implemented in R, based on the Bray–Curtis dissimilarity matrix. All libraries were subsampled to even depth (minimum sample read count, i.e., smallest library size, of 34,266) before analysis was performed. Similarity percentage analyses (SIMPER) was completed in R to reveal the most influential OTUs driving differences between soil bacterial communities in different locations, and across different seasons. Differences in relative abundances of taxa between locations and seasons were determined using Kruskal-Wallis rank sum non-parametric test. A post-hoc test was completed using pairwise Wilcoxon rank sum test if significant differences were reported using Kruskal-Wallis test. P values were adjusted using the method of Benjamini & Hochberg (1995) to control the false discovery rates (p < 0.05). Canonical analysis of principal coordinates (CAP) was done to relate environmental variables reported in Sintim et al. (2019) to changes in bacterial community composition. Our a priori hypothesis for conducting this statistical test was that we would see changes in microbial communities across the two locations driven by specific environmental variables. The ordination axes were constrained to linear combinations of environmental variables, then the environmental scores were plotted onto the ordination. A PERMANOVA was performed on the CAP axes. These analyses were completed in R following the online tutorial by Berry (2016).

Libraries were subsampled with replacement to equal size prior to computing alpha diversity metrics. The estimate_richness function was used in R phyloseq package to calculate observed richness and inverse Simpson indices (for diversity). A mixed model analysis of variance was completed using the generalized linear mixed model (GLIMMIX) procedure in SAS V. 9.3 to assess changes in richness and inverse Simpson over time. The fixed effects were location (TN and WA), mulch treatments (seven treatments) and time/season of soil sampling (four time points), while the random effect was block (three replicates). Repeated measures were incorporated in the model as sampling was done over time. The model was a completely randomized design (CRD) split-split-plot with repeated measures in the sub-sub plot. Normality of data was checked using Shapiro–Wilk test (W > 0.9) and equal variance using Levene’s test (α = 0.05). All data were normal and hence no transformations were performed.

To visualize differences in the functional profile of the communities; i.e., all seven enzyme rates, NMDS ordination of Bray-Curtis similarities was done in Primer 7 v. 7.0.13 (PRIMER-E). A mixed model analysis of variance with repeated measures was completed using the generalized linear mixed model (GLIMMIX) procedure in SAS V. 9.3 to assess changes in enzyme activities over time. Fixed and random effects were same as specified above. Boxplots for equal variance and outliers, reported in SAS, were used to remove outliers in the dataset. Normality was checked using Shapiro–Wilk test (W > 0.9) and probability plots for residuals, and equal variance using Levene’s test (α = 0.05). Data were log transformed as necessary when these conditions were not met. Type III tests of fixed effects and interaction effects are reported.

To assess for potential enrichment of bacteria and fungi, a paired t-test was conducted using initial and final 16S rRNA and ITS gene copy abundances (determined by qPCR) from Spring 2015 and Fall 2016 to see if there was a significant change over the course of the experiment. Initial 16S and ITS gene copy abundances from Spring 2015 were also subtracted from final abundances in Fall 2016 to get change in abundance over time. To determine if the enrichment or depletion of bacterial and fungal abundances was significantly different between treatments, a mixed model analysis of variance in SAS v. 9.3 using the GLIMMIX procedure was conducted on the differences. Significance level of all analyses were assessed at α = 0.05. All data were checked for normality using Shapiro–Wilk test (W > 0.9).

Environmental and soil physicochemical data

Environmental data collected during the experiment is reported in Sintim et al. (2019) and in Table S3. The mean daily air temperature in 2015 to 2016 was about 4 °C higher in Knoxville, TN than in Mount Vernon, WA (Table S3). The total annual precipitation during the experimental years was higher in Knoxville, TN than in Mount Vernon, WA.

In all plots with plastic mulching, fragments of the mulches (i.e., remnants from the previous season’s mulches) were visible in the soils throughout the experiment. Soil temperature, moisture and physicochemical properties were measured and reported previously by Sintim et al. (2019). In summary, significantly increased soil temperature was observed in the early growing seasons in the plastic mulch plots compared to the cellulose and no-mulch plots. On average, the monthly soil temperature was greater in TN than in WA. Overall, mulched plots had higher water content than the no mulch plots, with PE mulch having the highest soil water content for the greatest time. The soil health analysis revealed some effects of mulching on certain properties (namely aggregate stability, infiltration, soil pH, electrical conductivity, nitrate, and exchangeable potassium), but these were not consistent among BDMs, nor across sampling times and locations.

Soil bacterial community diversity and structure

For the 16S rRNA gene sequences, the percentage of bases with a Phred quality (Q) score ≥30 was 78% (24,090,356 total reads with 94% reads identified), and 90% (21,712,542 total reads with 93% reads identified) for the two flow cells used. The NMDS ordination revealed a clear difference in community structure between TN and WA when combining data from all four sampling seasons (Spring 2015 to Fall 2016) (Fig. 1A). Permutational ANOVA (PERMANOVA) tests confirmed significant differences between TN and WA soil microbial communities (Table 2, Table S4). The mean relative abundances of the most abundant classes of bacteria are shown in Fig. 1B. Similarity percentage tests (SIMPER) revealed the most influential OTUs contributing to the variation seen between location (Fig. 1B). The most influential OTUs belonged to several classes of microbes such as Acidobacteria_Gp7, Acidobacteria_Gp16, Acidobacteria_Gp4, Planctomycetacia and Spartobacteria. CAP analysis revealed that the differences in soil communities between TN and WA were most related soil moisture and organic matter content (Fig. S1).

In addition to locational differences, bacterial communities also differed significantly between the different seasons (Table 2, Table S4). For both locations, Spring communities were more similar to each other than Fall communities (Figs. 2A and 2B). SIMPER tests revealed that several genera of Acidobacteria, Planctomycetaceae, Spartobacteria and Actinobacteria (such as Streptomyces) were cumulatively responsible for 60% of the seasonal variance in bacterial communities (Figs. S2 and S3). Interestingly, relative abundance of Streptomyces increased over time from Spring 2015 to Fall 2016 in both TN and WA (Fig. S2).

Unlike location and season, the mulch treatments did not have a significant effect on bacterial community structure. Because of the locational and seasonal differences, we additionally analyzed each time-location set separately, and did not detect any significant effects of treatment on community structure except for Fall 2015 in WA (Fig. S4, Table 2, Table S4). We further used a pairwise.adonis function in R (Salazar, 2019) to determine pairwise differences for Fall 2015 in WA, but no significant patterns emerged.

Alpha diversity of the soil bacterial communities was estimated using observed species richness and inverse Simpson index of diversity (Table S5). The observed species richness estimator measures count of unique OTUs in each sample. There were significant differences between TN and WA (p < 0.05) in richness estimates (Table 3, Fig. 3A). TN had greater richness than WA throughout the experiment, ranging from 260 to 300 unique OTUs. WA richness estimates ranged from 250 to 280 OTUs over the two years. The locational differences in richness were due to a lower richness in Fall 2015, Spring 2016 and Fall 2016 in WA (Fig. 3A). The inverse Simpson diversity index ranges were similar between TN and WA, ranging from 7 to 11 (Fig. 3B).

For both TN and WA, there was a significant difference between seasons in both richness and inverse Simpson index (Table 3). In TN in 2016, PE had the lowest richness and BioAgri had the highest. However, treatment differences in richness estimates were not significant (Table 3) when analyzing data using a mixed model. Inverse Simpson diversity indices were also not significantly different between treatments (Table 3). Looking at the final time point in TN, diversity estimates were highest for Weedguard, and lowest for PE; in WA, the estimates were highest for Weedguard, followed by PE with BDMs having lower diversity than PE or Weedguard; however, these differences were not significant (Fig. 3B).

Microbial community abundances

As a proxy for bacterial and fungal abundances, bacterial (16S) and fungal (ITS) rRNA gene copies were quantified using qPCR assays for soil samples from all seasons. In general, bacterial (16S) gene copy numbers averaged to 2.5 × 10⁹ gene copies per gram dry weight soil in Spring 2015 and 4.2 × 10⁹ gene copies g⁻¹ in Fall 2016 in TN; 1.8 × 10⁹ gene copies g⁻¹ in Spring 2015 and 6.9 × 10⁹ gene copies g⁻¹ in Fall 2016 in WA. For fungal (ITS) abundances, values ranged from 2.9 × 10⁸ gene copies g⁻¹ in Spring 2015 to 3.8 × 10⁸ gene copies g⁻¹ in Fall 2016 in TN, and 4.5 × 10⁸ gene copies g⁻¹ in Spring 2015 to 8.4 × 10⁸ gene copies g⁻¹ in Fall 2016 in WA. In order to assess if gene abundances had significantly changed over the course of the experiment (Spring 2015 to Fall 2016) for each mulch treatment, a paired t-test was used to identify changes that are significantly different from zero (Table S6). There was a significant increase in bacterial gene copies under BDM and Weedguard treatments in WA, but no significant change for no mulch and PE treatments (Table S6). There was also a significant enrichment in fungal gene copies for two of the BDMs (PLA/PHA and Naturecycle) in WA. In TN, significant enrichment in bacterial gene copies was seen under Organix, PLA/PHA and PE treatments (Table S6) but no enrichment was seen in fungal gene copies. In order to determine if these changes were significantly different between treatments, the differences between the final (Fall 2016) and the initial (Spring 2015) abundances were analyzed using a mixed model analysis of variance and Tukey post hoc tests. In both locations, mulch treatments did not have a significant effect on the changes in either 16S or ITS gene copies over the course of the experiment (Figs. 4A and 4B).

Microbial community functions

To assess potential functional responses of the soil microbial communities to the mulching treatments, extracellular enzyme potential rate assays were conducted for common carbon, nitrogen, and phosphorus cycling enzymes in soil (Table S2). The data were combined over the two years to visualize Bray Curtis similarities of the enzyme rate profiles (Fig. 5). Locational differences in the enzyme profile were significant (p < 0.05), as were seasonal differences in both TN (p < 0.05) and WA (p < 0.05) evaluated using PERMANOVA (Fig. 5). However, mulch treatment did not have a significant effect on the enzyme profile for any of the seasons at either location (p < 0.05). NMDS ordination for the final sampling time point Spring 2017 is shown in Fig. S5, showing no clear treatment differences. In general, the enzyme activity rates oscillated between higher activities in the Spring and lower activities in the Fall. When analyzed separately for each enzyme, the data over the two years revealed a significant effect of sampling time in TN for all seven enzymes assayed. In WA, enzyme activities of β-xylosidase, β-glucosidase, α-glucosidase, N-acetyl β-glucosaminidase and phosphatase were significantly different between sampling times (Fig. 6). In WA, cellobiosidase and leucine amino peptidase activities remained unchanged across the seasons (10–22 nmol activity g⁻¹ dry soil h⁻¹ for cellobiosidase and 200–375 nmol activity g⁻¹ dry soil h⁻¹ for leucine amino peptidase) (Fig. 6).

When averaged across seasons, mulch treatment differences were not significant for any soil enzymes in WA (Table 4). However, in TN, an effect of mulch treatment was observed for N-acetylβ-glucosaminidase activities (Table 4). N-acetyl β-glucosaminidase activity was reduced under BDMs and PE compared to no mulch plots. Interaction effects of mulch treatment and time of sampling were not detectable for any of the enzymes assayed in TN or WA (Table 4).