Soil exposed to silver nanoparticles reveals significant changes in community structure and altered microbial transcriptional profiles

Journal Pre-proof Soil exposed to silver nanoparticles reveals significant changes in community structure and altered microbial transcriptional profiles Matthew J. Meier, Annette E. Dodge, Ajith Dias Samarajeewa, Lee A. Beaudette PII:

S0269-7491(19)35353-9

DOI:

https://doi.org/10.1016/j.envpol.2019.113816

Reference:

ENPO 113816

To appear in:

Environmental Pollution

Received Date: 19 September 2019 Revised Date:

9 December 2019

Accepted Date: 13 December 2019

Please cite this article as: Meier, M.J., Dodge, A.E., Samarajeewa, A.D., Beaudette, L.A., Soil exposed to silver nanoparticles reveals significant changes in community structure and altered microbial transcriptional profiles, Environmental Pollution (2020), doi: https://doi.org/10.1016/ j.envpol.2019.113816. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.

Traditional soil microbial health tests Dose-response modeling

Soil samples Pristine + Silver nanoparticles (multiple doses)

(Previous study)

Genomicsbased tests microbial sequencing using 16S and metatranscriptomics (This study)

Genomics-based tests are more sensitive than traditional soil microbial health tests

1 2 3 4 5

Soil Exposed to Silver Nanoparticles Reveals Significant Changes in Community Structure

6

and Altered Microbial Transcriptional Profiles

7 8

Matthew J. Meier, Annette E. Dodge, Ajith Dias Samarajeewa, Lee A. Beaudette

9 10

Biological Assessment and Standardization Section, Environment and Climate Change Canada,

11

335 River Road, Ottawa, Ontario, Canada K1V 1C7.

12 13

Corresponding author:

14

M.J. Meier

15

Biological Assessment and Standardization Section,

16

Environment and Climate Change Canada, 335 River Road,

17

Ottawa, Ontario, Canada, K1V 1C7.

18

Tel. (613) 949-9283

19

Email: [email protected]

20

21

Abstract

22

Anthropogenic activities can disrupt soil ecosystems, normally resulting in reduced soil

23

microbial health. Regulatory agencies need to determine the effects of uncharacterized

24

substances on soil microbial health to establish the safety of these chemicals if they end up in the

25

environment. Previous work has focused on measuring traditional ecotoxicologial endpoints

26

within the categories of microbial biomass, activity, and community structure/diversity. Because

27

these tests can be labor intensive, lengthy to conduct, and cannot measure changes in individual

28

gene functions, we wanted to establish whether metatranscriptomics could be used as a more

29

sensitive endpoint and provide a perspective on community function that is more informative

30

than taxonomic identification of microbes alone. We spiked a freshly collected sandy loam soil

31

(Vulcan, Alberta, Canada) with 0, 60, 145, 347, 833, and 2000 mg kg-1 of silver nanoparticles

32

(AgNPs), a known antagonist of microorganisms due to its propensity for dissolution of toxic

33

silver ions. Assessments performed in our previous work using traditional tests demonstrated the

34

toxicity of AgNPs on soil microbial processes. We expanded this analysis with genomics-based

35

tests by measuring changes in community taxonomic structure and function using 16S rDNA

36

profiling and metatranscriptomics. In addition to identifying bacterial taxa affected by AgNPs,

37

we found that genes involved in heavy metal resistance (e.g., the CzcA efflux pump) and other

38

toxicity response pathways were highly upregulated in the presence of silver. Dose-response

39

analysis using BMDExpress2 software successfully modeled many physiologically relevant

40

genes responding to low concentrations of AgNPs. We found that the transcriptomic point of

41

departure (BMD50) was lower than the IC50s calculated using the traditional tests in our previous

42

work. These results suggest that dose-response modeling of metatranscriptomic gene expression

43

is a useful tool in soil microbial health assessment.

44

Keywords: soil; bacteria; genomics; silver nanoparticles; dose-response modeling

45 46

Summary: Genomics-based endpoints for the assessment of soil microbial health can be used to

47

perform quantitative dose-response modeling, and soil-based RNAseq adds functional insights.

48

49

Introduction

50

Microbial activity in soil is critical for cycling nutrients and degrading organic

51

compounds (Schloter et al. 2018). Ensuring that soils contain healthy populations of microbes is

52

important from the perspective of numerous disciplines, such as agriculture (Sergaki et al. 2018)

53

and climate change (Cavicchioli et al. 2019), since healthy soils are the foundation of healthy

54

ecosystems. Pollution resulting in poor soil quality can have broad economic impacts – not only

55

because of remediation costs, but also because of the loss of primary productivity and the effects

56

on global climate. Therefore, it is important to understand and quantify the effects of xenobiotic

57

pollutants on soil microbes so that soil microbial health can be maintained.

58

Biological disturbances of soil microbial communities resulting from xenobiotic

59

pollutants can be quantified using ecotoxicological endpoints based on soil microbial functional

60

and structural measurements (Kvas et al. 2017). These methods are used for evaluating 1) the

61

health of soil samples taken from the field, or 2) the effects of uncharacterized chemicals on soil

62

microbial health. Because many of the tests are based on standardized methods, they may be

63

used by regulatory agencies in ecotoxicological risk assessments. However, since the soil

64

microbiome is extremely diverse and can vary drastically with biogeography, the incorporation

65

of genomics-based test methods would expand our ability to assess microbial responses to

66

xenobiotics, and perhaps pinpoint highly specific bioindicators for different types of adverse

67

effects on soil microbial health (Schloter et al. 2018).

68

Metals are a xenobiotic pollutant of significant concern in soils (Wuana and Okieimen

69

2011). In particular, various metal-based nanomaterials have emerged as industrially useful

70

materials in recent years. Silver nanoparticles (AgNPs) have a wide variety of applications,

71

stemming from their antimicrobial effects, including water treatment, textile antimicrobials, food

72

preservation, and medical sterilization (León-Silva et al. 2016). Although AgNPs may have

73

useful applications because of their antimicrobial effects (e.g., preventing disease or spoilage),

74

their presence in consumer products has led to the unintended release of silver nanoparticles into

75

ecosystems (Stensberg et al. 2011). The most important route through which AgNPs enter the

76

environment is through wastewater treatment plants. Almost 90% of AgNPs in the wastewater

77

are retained in sewer sludge, which is often used as fertilizer for crops or sent to landfills

78

(McGillicuddy et al. 2017). This raises concerns about the potential toxicity of silver

79

nanoparticles to the soil microbial community, especially in agriculture soils. Although AgNPs

80

are toxic to bacteria, a variety of stress response pathways have evolved to eliminate silver ions

81

(and other metals) from their cells (Hobman and Crossman 2015). Some genetic elements have

82

been characterized that confer silver resistance (e.g., the sil genes on pMG101); furthermore,

83

efflux pumps for other metal ions (e.g., copper) are able to provide some level of silver

84

resistance (Hobman and Crossman 2015). Previous work in our laboratory has demonstrated that

85

silver ions (Ag+) dissolved from AgNPs are toxic to soil microbial communities (Samarajeewa et

86

al. 2017); therefore, we are interested in understanding how treatment with AgNPs can influence

87

microbial functional dynamics in soil.

88

Since next-generation sequencing (NGS) has become common, some groups have used

89

soil microbiome sequencing to measure changes in response to environmental pollutants at the

90

community structure level (Ding et al. 2012; Fajardo et al. 2019) or at the gene function level

91

(Yuan et al. 2019; Sato et al. 2019; Yuan et al. 2017; Yergeau et al. 2012, 2014). Metagenome-

92

based methods are less biased than culture-based methods, but many studies to date focus only

93

on what taxa are present, because it is more expensive and technically challenging to measure

94

gene functions. However, understanding the distribution of gene functions is critical to our

95

understanding of dynamic changes in the microbial communities (Fondi et al. 2016). In marker

96

gene studies that attempt to impute gene functions, functional redundancy across taxa might

97

confound bacterial responses to pollutants, and provide an incomplete picture of community

98

function due to poor database annotation and difficulty in obtaining species- or strain-level

99

identifications. Furthermore, the quantification of gene functions using whole-metagenome

100

shotgun sequencing cannot measure what portions of the genome are actively transcribed in

101

response to environmental stimuli. To establish which genes are differentially regulated in

102

response to chemicals, the ideal method is to measure changes in gene expression by sequencing

103

RNA isolated directly from the environment. Uptake of this technology has been slow due to

104

difficulties in obtaining sufficient quantities of clean RNA and removing the rRNA fraction from

105

samples prior to building sequencing libraries. There are currently limited examples of

106

metatranscriptomics studies carried out in soils, though the list is continually growing (Sato et al.

107

2019; de Menezes et al. 2012; Carvalhais et al. 2012; Wegner and Liesack 2016; Yergeau et al.

108

2018; Jacquiod et al. 2018; Romero-Olivares et al. 2019; Barboza et al. 2018; Albright et al.

109

2018).

110

Toxicogenomics has emerged as a method to understand mechanisms of action driving

111

the disturbance of biological pathways, and often provides an accurate quantitative estimate of

112

the doses at which physiological endpoints are impacted. This type of data could be used to

113

develop risk assessments and regulations for chemicals of interest. Regulatory toxicologists use

114

benchmark dose (BMD) modeling to establish the dose at which a test chemical elicits a

115

response with a chosen magnitude (the “benchmark response”, e.g., 1 standard deviation relative

116

to control). In the context of toxicogenomics, this has been applied in a high-throughput manner

117

for gene expression data (Pagé-Larivière et al. 2019; Labib et al. 2017; Chepelev et al. 2016;

118

Webster et al. 2015) using BMDExpress2 (Phillips et al. 2018) where each gene is considered as

119

a separate endpoint to be modeled. Thus, genes with the lowest modeled BMDs are interpreted

120

as the most sensitive to the stimulus being measured. This information can be used to deduce

121

mode of action (Bourdon et al. 2013) or establish points of departure (Webster et al. 2015). This

122

approach has been widely used in mammalian toxicology studies, both in vivo and in vitro, as

123

well as ecotoxicology, but has not yet been adopted to measure the effects of environmental

124

pollutants on soil microorganisms.

125

The goal of this study was to determine if 16S community profiling and

126

metatranscriptomic sequencing methods could be used to quantitatively measure functional

127

changes within soil microbes in response to AgNPs. We also aimed to establish whether dose-

128

response modeling and pathway analysis could be applied to metatranscriptomic data to evaluate

129

soil microbial health in a manner similar to other fields of toxicology.

130

Methods

131

Soil treatment

132

Sandy-loam soil (0-15cm) was collected near the town of Vulcan, located in Alberta,

133

Canada in 2014, air dried briefly then sieved (2mm). This site was chosen because it is in the

134

boreal forest region of Canada and this soil is also used within our organization for invertebrate

135

soil toxicity testing (Velicogna et al. 2016; Jesmer et al. 2017). The test soils were amended with

136

increasing concentrations of pristine PVP coated AgNPs: 0, 60, 145, 347, 833, and 2000 mg Ag

137

kg-1 dry soil (nominal concentrations). Humic acid (1.8%) was used as a suspension agent for

138

AgNPs , Additional control amendments of PVP (equal to the amount of PVP present in 2000

139

mg Ag kg-1 samples) and humic acid (1.8%) were also included. Soil sample preparation and

140

treatment is described in detail in Samarajeewa et al. (2017). Briefly, soils amended in triplicates

141

were incubated for 63 days in the dark at 20°C. Samples collected from respective treatments at

142

the end of tests (i.e. 63 days after incubation) were stored at -20°C until DNA and RNA was

143

extracted.

144

Illumina MiSeq 16S Sequencing

145

DNA was extracted from 0.5 g of soil using the PowerSoil DNA Isolation Kit (MoBio

146

Laboratories Inc., CA, USA). The 16S V4 region was amplified using the universal forward

147

primer 515F-Y (5’- TCGTCGGCAGCGTCAGATGTGTATAAGAGA-

148

CAGGTGYCAGCMGCCGCGGTAA) and the reverse primer 926R (5’-

149

GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGCCGYCAATTYMTTTRAGTTT).

150

HotStarTaq Plus Master Mix (QIAGEN Inc., Germantown, MD, USA) was used with the

151

following theromocylcer program: 95 °C for 5 min; 40 cycles at 94 °C for 30 s, 50 °C for 30 s,

152

72 °C for 1 min; and a final elongation at 72 °C for 10 min. Amplicons were indexed using the

153

Nextera XT Index Primers (Illumina Inc., San Diego, CA, USA) and the KAPA HIFI HotStart

154

polymerase with the following thermocylcer program: 98 °C, 8 cycles of 30 sec at 98 °C, 30 sec

155

at 55 °C, 30 sec at 72 °C, and a final elongation step of 5 min at 72 °C. Indexed amplicons were

156

purified, pooled, and then quantified using a Qubit dsDNA BR Assay Kit (Life Technologies)

157

and combined at equimolar concentration before sequencing on an Illumina MiSeq (paired-end 2

158

X 250 bp). Sequencing was performed by the National Research Council Canada, Saskatoon.

159

Ion Torrent RNA Sequencing (Metatranscriptomics)

160

The concentrations 0, 60, 833, and 2000 mg Ag kg-1 were chosen for whole

161

metatranscriptome sequencing to ensure that effects at the lowest and highest doses were

162

captured in the study. Total RNA was extracted from 2 g of soil using the RNeasy Power Soil

163

Total RNA kit (Qiagen) and quantified using a Bioanalyser (Agilent® 2100, Agilent

164

Technologies, CA, USA) with the Agilent RNA 6000 Pico kit. Genomic DNA was removed by

165

treatment with DNase I at 37°C for one hour, and the reaction stopped by addition of EDTA and

166

heating to 65°C. RNA was purified with Agencourt AMPure XP PCR Purification magnetic

167

beads. The rRNA was removed from the samples using the RiboMinus Eukaryote System v2

168

(Life Technologies). A bacterial probe mix was used in place of the Eukaryote probe mix. The

169

ribodepleted RNA was concentrated using the Ion Torrent Magnetic Beads Cleanup Module.

170

Libraries were built using the Ion Total RNA-Seq Kit v2 for Whole Transcriptome Libraries

171

(Life Technologies) without the fragmentation step (the isolated RNA was already highly

172

fragmented, so further fragmentation was not required) and using unique barcodes for each

173

sample (IonXpress™ RNA Seq Barcode 01–16 Kit). Library quantification was done using a

174

Bioanalyser using the Agilent® High Sensitivity 1000 DNA kit (Agilent Technologies, CA,

175

USA). Finally, next-generation sequencing was completed using the Ion Torrent™ platform. The

176

Ion Chef System and an Ion PI™ Hi-Q Chef Kit was used to complete the emulsion PCR,

177

enrichment of Ion Sphere Particles, and loading of the PI™ sequencing chip. DNA sequencing

178

was done on the Ion Proton (Life Technologies) using the Ion PI™ Hi-Q™ Sequencing 200 kit,

179

and base calling was done using the Torrent Suite software (version 5.8.0).

180

16S rDNA Data Analysis

181

Demultiplexed FASTQ files were processed using the DADA2 (Callahan et al. 2016)

182

pipeline (version 1.8.0) using settings recommended by the authors with the exception of setting

183

truncLen to 240 and 200 for forward and reverse reads, respectively. Exact sequence variants

184

(ESVs) were classified using the SILVA (version 132) database. The ESVs were then imported

185

into R (R Core Team 2018) for statistical tests and visualization using Phyloseq version 1.24.2

186

(McMurdie and Holmes 2013). Sequences without a taxa assigned at the Phylum level were

187

removed from analysis, as were those with a prevalence of less than 0.47 (0.01 X # samples).

188

Ordination analysis was carried out in Phyloseq using principal components analysis on square-

189

root transformed ESV counts (to normalize their dispersion) with the weighted UniFrac distance.

190

The adonis function from the R package VEGAN version 2.5-4 (Dixon 2003) was used to

191

conduct permutational multivariate analysis of variance on the square-root transformed data (999

192

permutations, using the model formula: data ~ silver concentration + amendment + replicate) and

193

determine coefficients for the top taxa separating groups. Homogeneity of group variances was

194

tested using the betadisper function from VEGAN on a Bray-Curtis distance matrix of square-

195

root transformed ESV counts. The Phyloseq object was imported into DESeq2 version 1.20.0

196

(Love et al. 2014) to test for differentially abundant ESVs. An adjusted p-value cutoff of < 0.01

197

was used for reporting results.

198

We used PICRUSt2 version 2.1.4 beta (Langille et al. 2013)

199

(https://github.com/picrust/picrust2/) to predict the functional genes present in each sample based

200

on the 16S sequences identified and their relative proportions. The full pipeline script was run

201

using default settings on the ESVs called by DADA2 for each sample. Next, we used the

202

STAMP (Parks and Beiko 2010) workflow to test for differentially abundant features between

203

AgNP-exposed and unexposed soils. This analysis was performed at the gene-level (KEGG

204

orthologs) and the pathway-level (MetaCyc). For both the gene- and pathway-level analysis,

205

enriched functions were screened using a two groups approach (i.e., to detect differentially

206

abundant features between control soil and silver-exposed soils) and a multiple groups approach

207

(i.e., to detect differentially abundant between any of the silver concentrations). For the two

208

group tests, Welch’s two-sided unequal variances t-test was performed on control samples vs.

209

silver-exposed samples. Confidence intervals (99%) for the differences between mean

210

proportions were constructed using Welch’s inverted method, and p-values were corrected for

211

multiple comparisons using Storey’s FDR. Features were filtered to exclude those with a q-value

212

of <0.01, as well as those with an effect size of <2 (fold-change) or an absolute difference of

213

<0.04. For the multiple groups analysis, we used the Kruskal-Wallis H-test and the Games-

214

Howell post-hoc test (with 99% confidence intervals). The Benjamini-Hochberg FDR was used

215

to correct for multiple comparisons, since the p-values were not normally distributed. Features

216

were filtered to exclude those with a q-value of <0.01, and those with an η2 effect size of <0.95

217

(for gene-level annotations) or <0.90 (for pathway-level annotations).

218

Metatranscriptomics Data Analysis

219

FASTQ files from multiple sequencing runs were combined by their IonXpress barcoding

220

index. Data was uploaded to the MG-RAST server (Metagenomics Rapid Annotation using

221

Subsystem Technology) (Meyer et al. 2008) to obtain functional assignments of

222

metatranscriptomic reads. The SEED subsystems annotations were downloaded in BIOM format

223

and used in downstream analysis. DESeq2 was used for differential gene expression analysis to

224

determine if the counts for each transcript were significantly different across the experimental

225

conditions. DESeq2 performed fitting to the negative binomial distribution and calculation of the

226

Wald statistic. Further comparisons (i.e., contrasts in the DESeq2 results function) were made

227

relative to controls (no AgNPs). An adjusted p-value of 0.01 was used to filter results.

228

BMDExpress2 Modeling of Metatranscriptomic Data

229

BMDExpress2 is a software package used by regulatory agencies to analyze gene

230

expression data generated by RNASeq or microarray in a dose-response analysis. Data for each

231

gene (or probe) is fit to one of several parameterized models; the best model is selected for each

232

gene, and a benchmark dose (BMD; i.e., a dose at which a pre-defined response occurs, such as 1

233

standard deviation) is computed. Then, genes responding at the lowest doses can be targeted for

234

downstream analysis or grouped into pathways to establish which biological systems may be

235

disrupted by the test chemical. Raw feature counts from the metatranscriptomic sequencing

236

annotated by MG-RAST was log2-transformed and formatted for BMDExpress2 using a custom

237

script (https://github.com/mattjmeier/Meier_et_al_Metatranscriptomics). Prior to BMD

238

modeling, a one-way ANOVA filter was applied to remove genes with a fold-change of ≤2.0 and

239

a p-value of 0.05. No multiple testing correction was applied as it has been suggested that this

240

can reduce sensitivity (Pagé-Larivière et al. 2019). ANOVA was chosen instead of the Williams

241

trend test because non-monotonic responses are observed in our data in which a large drop

242

occurs at higher doses (likely due to taxonomic changes in the microbial community). BMD

243

modeling was done using the linear, exponential 2, polynomial 2, Hill, and power models, with a

244

benchmark response of 50% (i.e., 2.6 standard deviations), and no assumption of constant

245

variance. The best model for each gene was selected using the Nested Chi Square test with a p-

246

value of 0.05. Hill modeling that produced a ‘k’ parameter of 1/3 the lowest dose were flagged

247

and removed from analysis, selecting the next best model instead. Genes were filtered to remove

248

those where: the goodness of fit p-value was > 0.1, the BMD was < the highest dose, and the

249

BMDU to BMDL ratio was > 40, as per the recommendations set out by the National Toxicology

250

Program Approach to Genomic Dose-Response Modeling (National Toxicology Program 2018).

251

We did not remove genes with a BMD < ten times lower than the lowest dose because we were

252

interested in characterizing the lowest-responding genes, regardless of how well they were

253

modeled.

254

For defined category analysis, the SEED subsystems hierarchy (levels 2 and 3) was used

255

to create a custom category file and probe file as described in the BMDExpress2 documentation;

256

this was done to enable pathway-level tests. Since bacterial pathways in the SEED subsystems

257

exist on a much larger scale than mammalian gene sets (i.e., there are over 1,600 SEED

258

subsystems comprising >185,000 FIGfams with >16 million protein encoding genes (Overbeek

259

et al. 2014)), we did not filter on percentage of the pathway populated. However, we did limit the

260

analysis to subsystems with three or more genes modeled, and at least 1 per gene pathway

261

passing all filters.

262

A transcriptomic point of departure (POD) was calculated from the filtered BMDs using

263

the following metrics: median BMDs of significant genes; median of the 20 lowest BMDs of

264

significant genes; 10th percentile BMD of significant genes; and the lowest median BMD for

265

SEED subsystems (levels 2 and 3) with ≥3 genes modeled. These methods for POD calculations

266

were selected based on the results of Pagé-Larivière et al. (2019), which is the most current

267

recommendation available for calculating transcriptomic PODs with ecotoxicology data.

268

Data Availability

269

All raw data for 16S sequencing and metatranscriptomic sequencing are available in

270

GenBank under BioProject PRJNA565927 and MG-RAST project mgp86217, respectively. R

271

code used to produce statistical inferences and figures are available at

272

https://github.com/mattjmeier/Meier_et_al_Metatranscriptomics.

273 274

275

Results

276

16S Community Profiling Reveals Taxonomic Shifts in Response to AgNPs

277

The high-throughput sequencing results are summarized in Supplementary Table S1. We

278

obtained a raw average per-sample read count of 13,678 from 16S sequencing, calling an average

279

6,162 ESVs per sample. The bacterial communities in samples exposed to AgNPs are highly

280

impacted at the taxonomic structure level (Fig. 1). The samples cluster strongly by dose in

281

principal coordinates analysis (Fig. 1A), with the first axis driven by whether or not AgNPs were

282

present, and the second axis driven mostly by the dose of AgNPs present in the soil. The genera

283

found to be most significantly impacted (adjusted p-value <0.01 from DESeq test) are shown in

284

Fig. 1B (full data in Supplementary Table S2). Primarily, many ESVs of Rhodanobacter sp. were

285

enriched at low to medium concentrations of AgNPs, while ESVs classified as Sphingomonas

286

and Bradyrhizobium were enriched at medium to high concentrations. The genus Terrimonas

287

appears to be the most sensitive to the presence of AgNPs. These results are supported by the top

288

20 coefficients from the Adonis test (Fig. 1C). We found that the concentration of AgNPs was

289

found to be a significant source of variation between the samples (R2= 0.18, p<0.001). The

290

permutation test for homogeneity indicated that differences in dispersion was not a strong driver

291

of these observations (p=0.46). Thus, the concentration of silver was found to be a strong

292

determinant of the microbial community composition.

293

PICRUSt2 Functional Extrapolation from 16S Sequencing

294

To determine whether the observed taxonomic shifts in microbial communities might

295

alter their functional capabilities, we used PICRUSt2 to infer genomic content combined with

296

STAMP to determine whether functions were differentially abundant. A total of 6,383 KEGG

297

orthologs and 389 MetaCyc pathways were tested. When considering two experimental groups,

298

i.e., controls and AgNP-treated samples, there were 31 KEGG orthologs and 29 MetaCyc

299

pathways found to be significantly altered (Supplementary Fig. S1A and B, respectively; top 20

300

most significant are shown, with full data available in Supplementary Tables S3 and S4). When

301

testing for differences between any experimental groups, 44 KEGG orthologs and 31 MetaCyc

302

pathways were identified (Supplementary Fig. S2A and B, respectively; top 20 most significant

303

are shown). Between the two group test and multi-group test, we found two KEGG orthologs in

304

common (ATP-dependent Clp protease ATP-binding subunit ClpC and proton-dependent

305

oligopeptide transporter, POT family), and 16 MetaCyc pathways in common. Several KEGG

306

orthologs identified in this analysis are associated with oxidative stress (e.g., glutathione S-

307

transferase) or multi-drug efflux. Based on these inferred functions, there are clear differences

308

between control soil samples and AgNP-exposed soil samples.

309

Metatranscriptomic Sequencing

310

Because the accuracy of interpolated functions from metagenomes relies entirely on how

311

complete and accurate the underlying genome annotations are, we pursued whole-

312

metatranscriptome shotgun sequencing to quantify the expressed genes directly from our AgNP-

313

exposed soils. The concentrations used for metatranscriptomic sequencing were 0, 60, 833 and

314

2000 mg kg-1 AgNPs. We obtained a raw average per-sample read count of 13,723,789

315

(Supplementary Table S1). MG-RAST annotation resulted in the identification of SEED

316

subsystems for 3,699 genes with a median read count of 33,464 per sample, with approximately

317

10% of the identified genes present as singletons. The majority of the reads in each sample were

318

filtered out at the rRNA removal stage of the annotation pipeline; this is an ongoing challenge

319

with metatranscriptomic datasets, because ribodepletion is not capable of removing all the

320

contaminating rRNA, and poly-A based methods are not an option in bacterial systems.

321

DESeq2 was able to identify a total of 45 non-redundant genes across all samples that

322

were differentially expressed with an adjusted p-value <0.01 (Fig. 2; Supplementary Table S5).

323

These genes were mainly up-regulated relative to controls and fall into the categories of

324

oxidative stress response (e.g., superoxide dismutase) and efflux pumps such as CusA (a

325

monovalent cation ion exporter that is capable of Ag+ efflux). There were also various genes

326

with metal-associated roles, including divalent cation efflux pumps (CzcA) and other genes

327

associated with cobalt-zinc-cadmium resistance. These results, though somewhat limited in the

328

number of genes identified, suggest that the metatranscriptomic sequencing was able to pinpoint

329

multiple genes with known roles in metal stress response.

330

Dose response modeling using BMDExpress2

331

From the initial 3,699 genes, there were 367 that passed the ANOVA filter showing some

332

differential expression across at least one sample. Of those, 39 genes (Supplementary Table S6,

333

bolded entries; Fig. 3A and B) passed the other filtering criteria described in the methods as set

334

out by the NTP (goodness of fit p-value <0.1, etc.). Many of the gene-level BMDs were more

335

than 10-fold below the lowest dose in the study and were accordingly removed them from

336

analysis. We likely observed this because the lowest dose elicited a strong transcriptional

337

response (and strong taxonomic changes). Individual genes that showed a dose-response were

338

related to microbial stress response to metals (CzcB) or more general microbial stress response

339

(e.g., YbbN thioredoxin/chaperone, CmeA efflux pump, etc).

340

We grouped the dose-response modeled genes into pathways using the defined category

341

analysis function of BMDExpress2 (Fig. 3), using less strict filtering criteria than those set out

342

by the NTP as described in methods. Filtering criteria was relaxed because of the relatively low

343

number of genes that were well-modeled and the comparatively large database size of microbial

344

genes. The pathways with a median BMD <6.0 are highlighted in red in Fig. 3C; these groups

345

show the strongest response related to the dose of AgNPs. As observed in the DESeq2 results,

346

biologically relevant groupings emerge, such as resistance to antibiotics and toxic compounds,

347

oxidative stress (in subsystem level 2, the broader hierarchical level tested), and cobalt-zinc-

348

cadmium resistance in subsystem level 3 (the more specific hierarchical level).

349

The calculated transcriptomic points of departure (POD) are shown in Table 1. For the

350

most part, very low PODs were calculated using this approach. The most conservative estimate

351

was obtained using the lowest median BMD from SEED subsystems level 3, which was the

352

subsystem of cobalt-zinc-cadmium resistance.

353

Discussion

354

In previous work, the study and assessment of soil microbial health has mainly been

355

conducted using methods that rely on simple proxies for the disruption of microbial activity and

356

diversity, such as altered substrate decomposition, enumeration of colony forming units on agar

357

plates, or changes in community structure (Samarajeewa et al. 2017; Kvas et al. 2017). In this

358

study, we compare the use of 16S community profiling to metatranscriptomic sequencing in the

359

application of soil microbial health toxicity testing of AgNPs. We found that, while community

360

profiling is a highly sensitive measure for alterations that occurred as a response to AgNPs (Fig.

361

1, Supplementary Figs. S1 and S2), metatranscriptomics provided a clear biological picture of

362

the functional alterations that arose from the molecular stress response to silver ions (Fig. 2).

363

Additionally, through dose-response analysis of the metatranscriptomic gene expression profiles,

364

we found that BMD modeling delivered accurate biological indicators of metal stress in pathway

365

analysis (using SEED subsystems). This study provides an example of genomic-based toxicity

366

testing for soil microbial health.

367

The results of 16S community profiling demonstrate that exposure of soil to AgNPs

368

resulted in a strong selective pressure on the bacterial community (Fig. 1). Principal components

369

analysis revealed several distinct clusters that group by the specific dose of AgNPs. This trend is

370

also seen in the metatranscriptomic data (Fig. 2), indicating that both taxonomic and functional

371

changes occurred. An advantage of using 16S community profiling to characterize microbial

372

diversity is that experiments can be performed on a large number of samples with high

373

sequencing depth at low cost. For toxicity testing, this means that more doses can be tested,

374

which also means that it is more likely to properly model effects that occur between low, mid,

375

and high concentrations. With metatranscriptomic sequencing, because of higher costs,

376

researchers are limited to fewer samples with lower depth. In our data, there is an obvious trade-

377

off between the technologies. With 16S community profiling, we identified specific taxa (e.g.,

378

Rhodanobacter, Bradyrhizobium) that were tolerant of the AgNPs. It is unknown why some

379

species in the study were highly selected by AgNPs, but they could carry high numbers of genes

380

involved in metal resistance (Randall et al. 2015; Hobman and Crossman 2015) and are subject

381

to a combination of undetermined ecological factors beyond the scope of this work. Furthermore,

382

recent improvements in denoising algorithms means that taxonomic identifications using 16S

383

sequencing are improving (Porter and Hajibabaei 2018), though still fundamentally dependent on

384

the completeness of databases (which, for soil microbes, lack many representatives). However,

385

while taxonomic identification of silver-tolerant microbes may be of academic interest, it has

386

limited applicability to toxicity testing, since different soils are likely to contain communities of

387

different taxonomic structure (Fondi et al. 2016). On the other hand, functional redundancy

388

across different environments could mean that the gene expression changes in response to

389

xenobiotics are well-conserved on a biogeographic scale.

390

Metatranscriptomic sequencing identified many genes that are biologically relevant to

391

metal ion stress; most differentially expressed features were up-regulated (Fig. 2). The most

392

obviously relevant gene that was induced by AgNP toxicity is CusA, part of a tetrapartite efflux

393

system on the periplasmic membrane that has a high specificity for copper and silver.

394

Transcription of the cus locus is highly dependent on silver and copper (Delmar et al. 2015).

395

Likewise, a copper-translocating P-type ATPase, another exporter of metals (Alexander et al.

396

2011), was upregulated. Interestingly, AgNP toxicity also induced divalent cation efflux pumps

397

such as the cobalt-zinc-cadmium resistance proteins CzcA and CzcB. CzcA works as a cation-

398

proton antiporter while CzcB works as a cation-binding subunit. However, these proteins are

399

mainly characterized as transporters of the divalent cations of cobalt, zinc, and cadmium (Bruins

400

et al. 2000; Nies and Silver 1995), suggesting a cross-sensitivity of regulatory elements or

401

uncharacterized specificity of these genes. Differentially expressed genes also fell into the

402

category of extracellular sequestration of metal ions, such as OpgC (induced at 2000 ppm). This

403

membrane protein succinylates periplasmic glucans, conferring a negative charge and thereby

404

potentially immobilizing toxic cations to prevent them from diffusing into the cell; it was also

405

induced in other studies of metals (Maynaud et al. 2013). We observed alterations in oxidative

406

stress response genes such as superoxide dismutase, coenzyme PQQ synthesis protein A (Misra

407

et al. 2004), and 4-hydroxyphenylpyruvate dioxygenase (which is involved in tocopherol

408

biosynthesis, a known plant antioxidant; (Das et al. 2015)). Proteins involved in the repair of

409

proteins damaged by oxidation were induced, including the peptide methionine sulfoxide

410

reductase proteins MsrA and MsrB. Methionine is particularly susceptible to oxidation, and the

411

Msr enzymes catalyze the reduction of damaged residues back to methionine (Vattanaviboon et

412

al. 2005). Additionally, several cell wall/membrane components including various lipoproteins

413

were induced, which are important for maintaining envelope architecture and stability. Serine

414

acetyltransferase was upregulated, which is involved in polysaccharide biosynthesis (biofilm

415

formation), which may hinder AgNPs from diffusing to the cell membrane. Finally, we observed

416

induction of genes involved in flagellar components and chemotaxis, likely representing a

417

general stress response to the AgNPs. Overall, these numerous transcriptional responses make

418

sense in the context of metal stress, supporting the idea that metatranscriptomics can accurately

419

evaluate functional changes in the soil microbiome.

420

16S data, while useful as a marker gene to detect shifts in communities, performed less

421

accurately than metatranscriptomic sequencing when attempts were made to infer functional

422

changes. Metagenome prediction using PICRUSt2 is a useful approach, and in this study, it did

423

yield several statistically significant associations (Fig. 1 and Supplementary Figs. S1 and S2), but

424

can be limited for cases where database annotations and reference genomes are lacking, as is the

425

case with soils. Nonetheless, we did observe several potentially relevant genes (e.g., multidrug

426

efflux pumps), and some overlap with the metatranscriptomic results (e.g., the clpB/C genes

427

were downregulated in both datasets). We conclude that 16S community profiling coupled with

428

PICRUSt2 or similar methods for metagenome prediction could be used in toxicity testing when

429

metatranscriptomic sequencing is not possible. The results obtained using this method may be

430

real associations, but they are more difficult to interpret. This is because they are based on

431

observed selective pressure of taxonomic groups and, thus, unrelated functions (i.e., what one

432

could call passenger genes) found in the genomes of organisms that are selected will also

433

produce statistically significant results, even though they might not be actively transcribed.

434

The results of dose-response modeling with BMDExpress2 indicate that both gene- and

435

pathway-level groupings using SEED subsystems annotations provide a reliable way to

436

summarize metatranscriptomic data. BMDExpress2 transcriptomic dose-response modeling has

437

gained significant attention in mammalian toxicology (Labib et al. 2017) and is beginning to see

438

use in wildlife ecotoxicology (Pagé-Larivière et al. 2019); however, we are unaware of any

439

studies that have used transcription-based dose-response modeling in the context of soil

440

microbial health. Our hypothesis was that it could be used in a completely analogous way as

441

done for single-species tests, such that the lowest responding genes in the metatranscriptome

442

could be used to derive transcriptomic points of departure. In our previous work (Samarajeewa et

443

al. 2017), we calculated IC50s (the concentration causing a 50% inhibition, relative to controls)

444

for AgNPs using standardized soil toxicity tests. The IC50s in that study ranged from 19.9

445

(dehydrogenase enzyme activity test) to 173.8 mg kg-1 AgNPs (organic matter decomposition

446

test). Using transcriptomic PODs (based on a 50% benchmark response to enable comparison to

447

our previous IC50s), we obtained a range of 7.5X10-17 (cobalt-zinc-cadmium in level 3 SEED

448

subsystem) to 106.2 mg kg-1 AgNPs (sugar alcohols in level 2 SEED subsystem). This is a broad

449

range: the calculated POD depends highly on the exact group of genes that were modeled. Many

450

of the genes were completely undetected in the controls, leading to very low BMD estimates for

451

those endpoints. A limitation of using transcriptomic PODs is that the dose range of the study

452

should be carefully selected to ensure that low-dose effects are well modeled. In our study, the

453

doses were selected based on the results from pilot studies performed using the previously

454

described standardized tests and, therefore, we found that the transcriptomic effects are not well

455

modeled at low doses. In future studies, it would be useful to capture the lower end of the dose

456

response curve to more accurately estimate transcriptomic PODs. Nonetheless, the use of

457

BMDExpress2 described in this work presents a novel perspective for determining PODs using

458

metatranscriptomics that could eventually be pursued as a new standard test. Generally, we

459

obtained lower estimates for a 50% effective concentration using metatranscriptomics than we

460

did using traditional tests, which means that dose-response modeling of metatranscriptomic

461

datasets represents a more sensitive (and thus more conservative) method for calculating PODs.

462

We expect that the cobalt-zinc-cadmium subsystem contains genes that are highly

463

relevant to the mode of action of AgNP toxicity and, therefore, it is not surprising that it had a

464

low median BMD. However, other groupings of genes that we suspect are also relevant showed

465

higher median BMDs and might be comprised of genes that are modeled better (e.g., resistance

466

to antibiotics and toxic compounds in SEED subsystem 2 has a BMD median of 157.4 based on

467

12 genes that passed all filters). Furthermore, it is unclear what relevance the sugar alcohols

468

pathway (which had the lowest median BMD within the level 2 SEED subsystems) might have

469

in AgNP toxicity. This suggests that some subjective manual curation of good functional

470

pathway-level biomarkers might be necessary to calculate metatranscriptomic PODs that reflect

471

a particular stress response in microbial communities.

472

Several improvements could be made to this study. First, normalizing metatranscriptomic

473

gene expression to whole metagenome shotgun DNA-sequenced assemblies might improve the

474

detection of genes. However, DESeq2 is specifically designed as a statistical test method that can

475

detect between-sample changes based on count data, without the need for RPK or similar metrics

476

(Love et al. 2014). Another issue is limitations in the correct identification and annotation of

477

genes. This has been a problem for a long time and will continue to be one for many years, since

478

it is a very slow process to discern the biochemical roles of genes from uncultivated microbes.

479

Furthermore, due to infrastructure limitations, the soils in this study were stored at -20°C. It

480

would be preferable for RNA integrity to store them at -80°C, but this would be challenging in

481

many facilities, where freezer space is often at a premium. Additionally, higher coverage in the

482

metatranscriptomic sequencing might uncover more relevant functions from lower-abundance

483

organisms. Further to this point, it remains a major challenge to obtain good quality RNA from

484

soils; however, in our experience, with the soil we used, the kit used in this study (RNeasy Power

485

Soil Total RNA kit) gave sufficient quantities of total RNA using 2 g of soil. Other soil types

486

with a varying organic content and biomass might give different results. When extracting RNA

487

from challenging samples such as soil, it is critical to keep the work area free of RNAse

488

contamination, so that long fragments of RNA persist in the sample, which will improve

489

ribodepletion (due to the nature of probes only binding to relatively short target sequences).

490

Another important consideration is that samples exposed to higher concentrations of toxicant will

491

usually end up with more cell death, and therefore more fragmented RNA. This means that

492

higher doses may not provide as high-quality samples (this is not necessarily a problem since

493

transcriptional alterations will often be visible at much lower concentrations). Finally, there are

494

no international standard test guidelines for the use of genomics in soil microbial health toxicity

495

testing; therefore, the study of soil microbial health using genomics-based technologies would be

496

drastically improved by providing a guidance document, as done by the NTP for genomics-based

497

dose-response analysis (National Toxicology Program 2018).

498

Although we have provided a cursory analysis of the known observed metal-responsive

499

genes, the data produced in this study may contain differentially expressed genes with unknown

500

functions in stress responses that other researchers could pursue in future studies. Furthermore,

501

the effect size in this study was large even at the lowest dose – it might be beneficial to test lower

502

doses of metals in future work (using pilot studies that perform range finding using genomics-

503

based endpoints instead of only traditional test methods). Finally, it is important to test other

504

chemicals and other soils in future studies, so that we can identify whether transcriptomic

505

responses are well-conserved across different soils, and whether chemicals of different classes

506

can produce a reliable signature of adverse soil microbial health.

507

In this study, in response to AgNP soil amendment, selection for organisms carrying

508

adaptive genes occurred, as did gene expression changes that constitute a response to metal

509

stress. We demonstrate that, for AgNP exposures, examination of the differentially expressed

510

genes clearly points to specific functional classifications involved in metal ion stress that are

511

disrupted. Based on our results, we suggest that metatranscriptomic sequencing could play two

512

roles in regulatory toxicology. First, it could generate mechanistic hypotheses that could later be

513

confirmed using chemical analysis for uncharacterized soil samples that may be adversely

514

impacted (given a reference site for comparison). Second, it could be used in the context of

515

toxicity testing to perform risk assessments for uncharacterized chemicals based on dose-

516

response analysis of gene expression and pathway analysis, and their calculated points of

517

departure. However, we caution that – despite ever-falling sequencing costs – the technology to

518

perform metatranscriptomics is not yet easily accessible to most labs, due to limitations in

519

analytical (i.e., computational) requirements, limited budgets, and a lack of standard test

520

methods. We anticipate that these methods will continue to improve as technological advances

521

continue and sequencing costs make these studies more affordable, enabling the characterization

522

of a wide range of soils and chemicals.

523

Acknowledgements

524

This work was supported by the Government of Canada’s Interdepartmetnal - Genomics

525

Research and Development Initiative (GRDI) EcoBiomics project. We are thankful for technical

526

contributions from Dr. Armand Seguin and Dr. Christine Martineau, as well as comments on the

527

manuscript from Dr. Rick Scroggins.

528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571

References Albright MB, Johansen R, Lopez D, Steven B, Kuske CR, Dunbar J, others. 2018. Short-term transcriptional response of microbial communities to nitrogen fertilization in a pine forest soil. Appl Environ Microbiol 84: e00598–18. Alexander SPH, Mathie A, Peters JA. 2011. Guide to Receptors and Channels (GRAC), 5th edition. Br J Pharmacol 164 Suppl 1: S1–324. Barboza ADM, Pylro VS, Jacques RJS, Gubiani PI, de Quadros FLF, da Trindade JK, Triplett EW, Roesch L. 2018. Seasonal dynamics alter taxonomical and functional microbial profiles in Pampa biome soils under natural grasslands. PeerJ 6: e4991. Bourdon JA, Williams A, Kuo B, Moffat I, White PA, Halappanavar S, Vogel U, Wallin H, Yauk CL. 2013. Gene expression profiling to identify potentially relevant disease outcomes and support human health risk assessment for carbon black nanoparticle exposure. Toxicology 303: 83–93. Bruins MR, Kapil S, Oehme FW. 2000. Microbial resistance to metals in the environment. Ecotoxicol Environ Saf 45: 198–207. Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. 2016. DADA2: High-resolution sample inference from Illumina amplicon data. Nat Methods 13: 581–3. Carvalhais LC, Dennis PG, Tyson GW, Schenk PM. 2012. Application of metatranscriptomics to soil environments. J Microbiol Methods 91: 246–51. Cavicchioli R, Ripple WJ, Timmis KN, Azam F, Bakken LR, Baylis M, Behrenfeld MJ, Boetius A, Boyd PW, Classen AT, et al. 2019. Scientists’ warning to humanity: microorganisms and climate change. Nat Rev Microbiol. Chepelev NL, Long AS, Williams A, Kuo B, Gagné R, Kennedy DA, Phillips DH, Arlt VM, White PA, Yauk CL. 2016. Transcriptional Profiling of Dibenzo[def,p]chrysene-induced Spleen Atrophy Provides Mechanistic Insights into its Immunotoxicity in MutaMouse. Toxicol Sci 149: 251–68. Das P, Nutan KK, Singla-Pareek SL, Pareek A. 2015. Oxidative environment and redox homeostasis in plants: dissecting out significant contribution of major cellular organelles. Frontiers in Environmental Science 2: 1–11. Delmar JA, Su C-C, Yu EW. 2015. Heavy metal transport by the CusCFBA efflux system. Protein Sci 24: 1720–36. Ding G-C, Heuer H, Smalla K. 2012. Dynamics of bacterial communities in two unpolluted soils after spiking with phenanthrene: soil type specific and common responders. Front Microbiol

572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617

3: 1–16. Dixon P. 2003. VEGAN, a package of R functions for community ecology. Journal of Vegetation Science 14: 927–930. Fajardo C, Costa G, Nande M, Bot’\ias P, Garc’\ia-Cantalejo J, Mart’\in M. 2019. Pb, Cd, and Zn soil contamination: Monitoring functional and structural impacts on the microbiome. Applied soil ecology 135: 56–64. Fondi M, Karkman A, Tamminen MV, Bosi E, Virta M, Fani R, Alm E, McInerney JO. 2016. “Every Gene Is Everywhere but the Environment Selects”: Global Geolocalization of Gene Sharing in Environmental Samples through Network Analysis. Genome Biol Evol 8: 1388– 400. Hobman JL, Crossman LC. 2015. Bacterial antimicrobial metal ion resistance. Journal of medical microbiology 64: 471–497. Jacquiod S, Nunes I, Brejnrod A, Hansen MA, Holm PE, Johansen A, Brandt KK, Priemé A, Sørensen SJ. 2018. Long-term soil metal exposure impaired temporal variation in microbial metatranscriptomes and enriched active phages. Microbiome 6: 1–14. Jesmer AH, Velicogna JR, Schwertfeger DM, Scroggins RP, Princz JI. 2017. The toxicity of silver to soil organisms exposed to silver nanoparticles and silver nitrate in biosolids-amended field soil. Environ Toxicol Chem 36: 2756–2765. Kvas S, Rahn J, Engel K, Neufeld JD, Villeneuve PJ, Trevors JT, Lee H, Scroggins RP, Beaudette LA. 2017. Development of a microbial test suite and data integration method for assessing microbial health of contaminated soil. J Microbiol Methods 143: 66–77. Labib S, Williams A, Kuo B, Yauk CL, White PA, Halappanavar S. 2017. A framework for the use of single-chemical transcriptomics data in predicting the hazards associated with complex mixtures of polycyclic aromatic hydrocarbons. Arch Toxicol 91: 2599–2616. Langille MGI, Zaneveld J, Caporaso JG, McDonald D, Knights D, Reyes JA, Clemente JC, Burkepile DE, Vega Thurber RL, Knight R, et al. 2013. Predictive functional profiling of microbial communities using 16S rRNA marker gene sequences. Nat Biotechnol 31: 814–21. León-Silva S, Fernández-Luqueño F, López-Valdez F. 2016. Silver nanoparticles (AgNP) in the environment: a review of potential risks on human and environmental health. Water, Air, \& Soil Pollution 227: 1–20. Love MI, Huber W, Anders S. 2014. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15: 550. Maynaud G, Brunel B, Mornico D, Durot M, Severac D, Dubois E, Navarro E, Cleyet-Marel JC, Le Quéré A. 2013. Genome-wide transcriptional responses of two metal-tolerant symbiotic

618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662

Mesorhizobium isolates to zinc and cadmium exposure. BMC Genomics 14: 292. McGillicuddy E, Murray I, Kavanagh S, Morrison L, Fogarty A, Cormican M, Dockery P, Prendergast M, Rowan N, Morris D. 2017. Silver nanoparticles in the environment: Sources, detection and ecotoxicology. Sci Total Environ 575: 231–246. McMurdie PJ, Holmes S. 2013. phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PLoS ONE 8: e61217. De Menezes A, Clipson N, Doyle E. 2012. Comparative metatranscriptomics reveals widespread community responses during phenanthrene degradation in soil. Environ Microbiol 14: 2577– 88. Meyer F, Paarmann D, D’Souza M, Olson R, Glass EM, Kubal M, Paczian T, Rodriguez A, Stevens R, Wilke A, et al. 2008. The metagenomics RAST server - a public resource for the automatic phylogenetic and functional analysis of metagenomes. BMC Bioinformatics 9: 386. Misra HS, Khairnar NP, Barik A, Indira Priyadarsini K, Mohan H, Apte SK. 2004. Pyrroloquinoline-quinone: a reactive oxygen species scavenger in bacteria. FEBS Lett 578: 26–30. National Toxicology Program. 2018. NTP research report on national toxicology program approach to genomic dose-response modeling. Nies DH, Silver S. 1995. Ion efflux systems involved in bacterial metal resistances. J Ind Microbiol 14: 186–99. Overbeek R, Olson R, Pusch GD, Olsen GJ, Davis JJ, Disz T, Edwards RA, Gerdes S, Parrello B, Shukla M, et al. 2014. The SEED and the Rapid Annotation of microbial genomes using Subsystems Technology (RAST). Nucleic Acids Res 42: D206–14. Pagé-Larivière F, Crump D, O’Brien JM. 2019. Transcriptomic points-of-departure from shortterm exposure studies are protective of chronic effects for fish exposed to estrogenic chemicals. Toxicol Appl Pharmacol 378: 114634. Parks DH, Beiko RG. 2010. Identifying biologically relevant differences between metagenomic communities. Bioinformatics 26: 715–21. Phillips JR, Svoboda DL, Tandon A, Patel S, Sedykh A, Mav D, Kuo B, Yauk CL, Yang L, Thomas RS, et al. 2018. BMDExpress 2: enhanced transcriptomic dose-response analysis workflow. Bioinformatics. Porter TM, Hajibabaei M. 2018. Scaling up: A guide to high-throughput genomic approaches for biodiversity analysis. Mol Ecol 27: 313–338.

663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708

R Core Team. 2018. R: A Language and Environment for Statistical Austria http://www.R-project.org.

Computing. Vienna,

Randall CP, Gupta A, Jackson N, Busse D, O’Neill AJ. 2015. Silver resistance in Gram-negative bacteria: a dissection of endogenous and exogenous mechanisms. Journal of Antimicrobial Chemotherapy 70: 1037–1046. Romero-Olivares AL, Meléndrez-Carballo G, Lago-Lestón A, Treseder KK. 2019. Soil metatranscriptomes under long-term experimental warming and drying: fungi allocate resources to cell metabolic maintenance rather than decay. Frontiers in Microbiology 10: 1914. Samarajeewa AD, Velicogna JR, Princz JI, Subasinghe RM, Scroggins RP, Beaudette LA. 2017. Effect of silver nano-particles on soil microbial growth, activity and community diversity in a sandy loam soil. Environ Pollut 220: 504–513. Sato Y, Hori T, Koike H, Navarro RR, Ogata A, Habe H. 2019. Transcriptome analysis of activated sludge microbiomes reveals an unexpected role of minority nitrifiers in carbon metabolism. Commun Biol 2: 179. Schloter M, Nannipieri P, Sørensen SJ, van Elsas JD. 2018. Microbial indicators for soil quality. Biology and Fertility of Soils 54: 1–10. Sergaki C, Lagunas B, Lidbury I, Gifford ML, Schäfer P. 2018. Challenges and Approaches in Microbiome Research: From Fundamental to Applied. Front Plant Sci 9: 1205. Stensberg MC, Wei Q, McLamore ES, Porterfield DM, Wei A, Sepúlveda MS. 2011. Toxicological studies on silver nanoparticles: challenges and opportunities in assessment, monitoring and imaging. Nanomedicine (Lond) 6: 879–98. Vattanaviboon P, Seeanukun C, Whangsuk W, Utamapongchai S, Mongkolsuk S. 2005. Important role for methionine sulfoxide reductase in the oxidative stress response of Xanthomonas campestris pv. phaseoli. J Bacteriol 187: 5831–6. Velicogna JR, Ritchie EE, Scroggins RP, Princz JI. 2016. A comparison of the effects of silver nanoparticles and silver nitrate on a suite of soil dwelling organisms in two field soils. Nanotoxicology 10: 1144–51. Webster AF, Chepelev N, Gagné R, Kuo B, Recio L, Williams A, Yauk CL. 2015. Impact of Genomics Platform and Statistical Filtering on Transcriptional Benchmark Doses (BMD) and Multiple Approaches for Selection of Chemical Point of Departure (PoD). PLoS ONE 10: e0136764. Wegner C-E, Liesack W. 2016. Microbial community dynamics during the early stages of plant polymer breakdown in paddy soil. Environ Microbiol 18: 2825–42.

709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731

Wuana RA, Okieimen FE. 2011. Heavy metals in contaminated soils: a review of sources, chemistry, risks and best available strategies for remediation. Isrn Ecology 2011. Yergeau E, Sanschagrin S, Beaumier D, Greer CW. 2012. Metagenomic analysis of the bioremediation of diesel-contaminated Canadian high arctic soils. PLoS ONE 7: e30058. Yergeau E, Sanschagrin S, Maynard C, St-Arnaud M, Greer CW. 2014. Microbial expression profiles in the rhizosphere of willows depend on soil contamination. ISME J 8: 344–58. Yergeau E, Tremblay J, Joly S, Labrecque M, Maynard C, Pitre FE, St-Arnaud M, Greer CW. 2018. Soil contamination alters the willow root and rhizosphere metatranscriptome and the root-rhizosphere interactome. ISME J 12: 869–884. Yuan K, Chen B, Qing Q, Zou S, Wang X, Luan T. 2017. Polycyclic aromatic hydrocarbons (PAHs) enrich their degrading genera and genes in human-impacted aquatic environments. Environ Pollut 230: 936–944. Yuan L, Li Z-H, Zhang M-Q, Shao W, Fan Y-Y, Sheng G-P. 2019. Mercury/silver resistance genes and their association with antibiotic resistance genes and microbial community in a municipal wastewater treatment plant. Science of The Total Environment 657: 1014–1022.

732

Table 1.

733

Transcriptomic points of departure for AgNP toxicity using BMDs calculated using

734

BMDExpress2 summarized by different approaches. Median BMD (all significant genes)

2.85E-04

735 736

Median 20 lowest BMDs

3.75E-07

10th percentile BMD

Lowest median BMD in SEED Level 3 (4 genes in cobalt-zinc-cadmium resistance) 8.25E-10 7.50E-17

Lowest median BMD for SEED Level 2 (7 genes in sugar alcohols) 106.19

737

Figure 1.

738

Introduction of silver nanoparticles (AgNPs) induced significant changes in the

739

taxonomic structure of the soil (measured by 16S rDNA profiling). This is consistent with

740

previous studies demonstrating that AgNPs adversely affect soil health.

741 742 743 744 745 746 747

A) PCoA plot (weighted UniFrac distance) demonstrating a clear separation (Axis 1) of silver-exposed vs. pristine soil conditions. B) Heatmap showing taxa with statistically significant differences in relative abundance (adjusted P-value <0.01). C) Coefficients from multivariate permutational analysis showing the top taxa contributing to variability between each experimental group.

748 749

750

Figure 2.

751

Functional alterations at the transcriptome level are evident when AgNPs are added at

752

different concentrations. By annotating the gene functions from metatranscriptomic sequencing

753

of unamended soil and low, medium, and high dose groups of AgNPs, we detected statistically

754

significant differences between the treatments that correspond to biological adverse effects.

755

Numerous genes involved in metal efflux, oxidative stress, and membrane integrity were

756

differentially expressed in the presence of AgNPs.

757

Figure 3.

758

BMDExpress2 was used to model dose-response curves for each gene in the

759

metatranscriptomics data. A) An approximate 1:1 relationship was observed between BMD and

760

BMDL, indicating that genes passing filters were well-modeled. B) Expression profiles for the

761

39 genes that were well-modeled (excluding those with a BMD < 10 times the lowest dose).

762

Genes involved in metal stress response (e.g., cobalt-zinc-cadmium efflux, other efflux pumps)

763

were present. C) The accumulation plot at the pathway level (SEED subsystems 2 and 3)

764

indicates that biologically relevant groupings of genes (3 or more genes per pathway) were

765

present. Pathways highlighted in red have an average BMD <6.0 (ten times the lowest dose in the

766

study), which indicates that some genes in these pathways may be poorly modeled; however, in

767

this study they were still found to contain relevant groupings (e.g., oxidative stress, cobalt-zinc-

768

cadmium resistance, resistance to toxic compounds).

769

Highlights

• • •

We used RNAseq to measure the effect of silver nanoparticles on soil bacteria Gene expression in soil bacteria can be used to perform dose-response modeling Transcriptional profiles from silver exposure reveal mechanisms of stress response

Author Contributions MJM: Conceptualization, Investigation, Formal analysis, Data curation, Visualization, Supervision, Writing AED: Investigation, Writing – Original Draft ADS: Investigation, Writing – Review and Editing LAB: Conceptualization, Supervision, Funding acquisition, Project administration, Writing – Review and Editing

Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: