Unearthing the Boreal Soil Resistome Associated with Permafrost Thaw

Understanding the distribution and mobility of antibiotic resistance genes (ARGs) in soil bacteria from diverse ecological niches is critical in assessing their impacts on the global spread of antibiotic resistance. In permafrost associated soils, climate and human driven forces augment near-surface thaw altering the overlying active layer. Physiochemical changes shift bacterial community composition and metabolic functioning, however, it is unknown if permafrost thaw will affect ARGs comprising the boreal soil resistome. To assess how thaw shifts the resistome, we performed susceptibility testing and whole genome sequencing on soil isolates from a disturbance-induced thaw gradient in Interior Alaska. We found resistance was widespread in the Alaskan isolates, with 87% of the 90 isolates resistant to at least one of the five antibiotics. We also observed positive trends in both the proportion of resistant isolates and the abundance of ARGs with permafrost thaw. However, the number of ARGs per genome and types of genes present were shown to cluster more strongly by bacterial taxa rather than thaw emphasizing the evolutionary origins of resistance and the role vertical gene transfer has in shaping the predominantly chromosomally encoded ARGs. The observed higher proportion of plasmid-borne and distinct ARGs in our isolates compared to RefSoil+ suggests local conditions affect the composition of the resistome along with selection for ARG mobility. Overall taxonomy and geography shape the resistome, suggesting that as microbial communities shift in response to permafrost thaw so will the ARGs in the boreal active layer. IMPORTANCE As antibiotic resistance continues to emerge and rapidly spread in clinical settings, it is imperative to generate studies that build insight into the ecology of environmental resistance genes that pose a threat to human health. This study provides insight into the occurrence of diverse ARGs found in Alaskan soil bacteria which is suggestive of the potential to compromise health. The observed differences in ARG abundance with increasing permafrost thaw suggest the role of soil disturbance in driving the distribution of resistant determinants and the predominant taxa that shape the resistome. Moreover, the high-quality whole genome assemblies generated in this study are an extensive resource for microbial researchers interested in permafrost thaw and will provide a steppingstone for future research into ARG mobility and transmission risks.

unknown if permafrost thaw will affect ARGs comprising the boreal soil resistome. To 23 assess how thaw shifts the resistome, we performed susceptibility testing and whole 24 genome sequencing on soil isolates from a disturbance-induced thaw gradient in Interior 25 Alaska. We found resistance was widespread in the Alaskan isolates, with 87% of the 90 26 isolates resistant to at least one of the five antibiotics. We also observed positive trends 27 in both the proportion of resistant isolates and the abundance of ARGs with permafrost 28 thaw. However, the number of ARGs per genome and types of genes present were 29 shown to cluster more strongly by bacterial taxa rather than thaw emphasizing the 30 evolutionary origins of resistance and the role vertical gene transfer has in shaping the 31 predominantly chromosomally encoded ARGs. The observed higher proportion of 32 plasmid-borne and distinct ARGs in our isolates compared to RefSoil+ suggests local 33 conditions affect the composition of the resistome along with selection for ARG mobility. 34 Overall taxonomy and geography shape the resistome, suggesting that as microbial 35 communities shift in response to permafrost thaw so will the ARGs in the boreal active 36 layer. 37 38 Permafrost thaw associated Boreal Soil Resistome 3 IMPORTANCE. 39 As antibiotic resistance continues to emerge and rapidly spread in clinical settings, it is 40 imperative to generate studies that build insight into the ecology of environmental 41 resistance genes that pose a threat to human health. This study provides insight into the 42 occurrence of diverse ARGs found in Alaskan soil bacteria which is suggestive of the 43 potential to compromise health. The observed differences in ARG abundance with 44 increasing permafrost thaw suggest the role of soil disturbance in driving the distribution 45 of resistant determinants and the predominant taxa that shape the resistome. Moreover, 46 the high-quality whole genome assemblies generated in this study are an extensive 47 resource for microbial researchers interested in permafrost thaw and will provide a 48 steppingstone for future research into ARG mobility and transmission risks.

Bacterial Culturing. 182
In the September of 2017 we collected two 10 cm wide by 20 cm deep soil cores 183 from each FPES site with a sterilized corer. Prior to coring the top layer of moss and 184 vegetation were removed. Soil cores were then extracted and immediately stored in a 185 cooler throughout sample collection. Soil cores were moved to a +4°C fridge until 186 processing 24 h later. In order to process the cores in a way that prevented 187 contamination from any exogenous cells on the exterior of the soil core, outer portion of 188 each core was removed using a sterile scalpel. The interiors of each core were then 189 sub-sample using sterile forceps along a depth gradient at intervals of about 2.5 cm for 190 20 cm to generate a total of 1 g of soil. The 1 g of soil was used to inoculate 100ml 191 tryptic soy broth (TSB) to produce an enrichment culture. After 48 h at 22°C, we plated 192 serial dilutions (1:10, 1:100, 1:1000) of the enrichment culture three times for each 193 sample and incubated the plates at three temperatures (+4˚C, +12˚C, and +20 ˚C), for a 194 total of twelve plates per sample, until distinct colony formation was observed. We used a Nextera XT library (Illumina, San Diego, CA, USA) prepared by the 218 Genomics Core Lab at the University of Alaska Fairbanks to sequence on an Illumina 219 MiSeq platform using version 3 reagents. We trimmed adapters with TrimGalore version 220 0.5.0 (https://github.com/FelixKrueger/TrimGalore). For long reads, a combination of 221 SQK-RBK004 and VSK-VSK002 library preparation (Oxford Nanopore Technology 222 [ONT]) (Table S2)  unlike any of the other databases was able to accurately identify the maximum number 250 of resistance genes from the whole genome sequences of 3 strains of methicillin 251 resistant Staphylococcus aureus (Xavier et al. 2006). In this study we annotated each 252 whole genome assembly with CARD version 3.0.9 using command line tool Resistance 253 Gene Identifier (RGI) version 5.1.0 Main specifying input type contig (-t contig) with 254 default parameters for BLAST alignment (-a BLAST) and strict and perfect hits only by 255 not choosing the include loose option. We chose to use the strict algorithm rather than 256 the perfect or loose because it can detect previously unknown homologs using detection 257 models with curated similarity cut-offs while still ensuring the detected variant is a 258 functional resistance gene rather than spurious partial hits (Jia et al. 2006). Results from 259 RGI were further quality controlled by removing any ARO hits defined as mutations or 260 ARO hits with less than 50% coverage of the reference sequence unless cutoff on the 261 edge of a contig in order to ensure hits were ARG homologs rather than spurious hits. 262 To determine if a hit was located on a chromosome or plasmid, contigs containing hits 263 were run through BLASTn (Zhang et al. 2000) (Figure 2). 304

Antibiotic Resistance Genes. 311
Across all FPES genomes RGI identified 379 significant hits comprising 27 312 CARD-based AROs (Table S1). Of these hits 30 had 100% sequence identity to CARD 313 AROs and another 32 were highly similar (sequence identity >90%). Four genes hits, 314 two encoding aminoglycoside inactivating enzymes AAC(6')-32 and AAC(6')-Ir and two 315 encoding bcrC, an undecaprenyl pyrophosphate related protein, had full length coverage 316 of the reference sequence along with 100% sequence identity. These findings show that 317 some ARGs present in these soil isolates are highly homologous to those present in the 318 CARD database rather than just ancient divergent homologs. However, overall mean 319 percent identity across all hits identified was 68.4 % ± 19 suggesting that although some 320 hits were highly similar to known resistance genes, there were also many hits identified 321 that diverge from ARGs present in the CARD database suggesting novel homologs of 322 known resistance genes are present in these isolates. Bacillus isolates across all thaw sites. We also observed beta-lactam inactivating genes 346 from four distinct beta-lactamase families (Table S1). Bc beta-lactamase had the highest 347 abundance of gene copies (n=20) all encoding BcII, a zinc metallo-beta-lactamase that 348 hydrolyzes a large number of penicillins and cephalosporins. BcII gene copies in our 349 samples were confined to the genus Bacillus and found to be both highly similar (% 350 identity = 90.755 ± 0.567) to BcII homologs in the CARD database and had high gene

Influence of Permafrost Thaw Distrubance 367
We found an increasing abundance of ARGs associated with increasing thaw 368 distrubance (UD=101, SD=133, MD=145). We found that the mean number of CARD 369 hits per isolate was highly significant by genus (Kruskal-Wallis p=1.9e-12; Figure 3B) 370 and significant between the undisturbed and most-disturbed sites (Wilcoxon p = 0.045). 371 The significant difference in ARG copies per isolate by genus suggests that the 372 taxonomy plays a significant role in shaping abundance of resistance genes and thus 373 with a shift in a community composition can affect the abundance and types of 374 resistance genes. When examining the types of genes present within isolates, there is a 375 distinct set of core ARGs found within each taxon rather than between taxa across sites 376 such as FosB, BcII in Bacillus isolates and adeF, armA, and soxR in Pseudomonas 377 ( Figure 4). The set of intrinsic ARGs within each taxon caused groups to cluster more strongly by phylogeny rather than FPES site (Figure 4). Overall Bacillus isolates' core 379 resistance genes were comprised of primarily antibiotic inactivation genes, 380 Pseudomonas was antibiotic efflux and target alteration, and Erwinia contained antibiotic 381 efflux, target alteration, and inactivation. Although there were many genes found across 382 all sites within some taxa, some genes were observed exclusively in one taxon and site 383 such as MCR-4.1 in Bacillus from the semi-disturbed site. Although rare, these less 384 abundant genes, and thus the ARGs with less of an influence on clustering, are more 385 likely resistance determinants that were acquired either through conjugation, 386 transformation, or transduction from members within the soil community and potentially 387 more of an interest for assessing clinical risk. 388

RefSoil+ Comparison 389
Across the genera present in our FPES isolates, RefSoil+ and FPES had 15 390 similar ARGs, 41 unique to RefSoil+ and 12 unique to FPES ( Figure 5A). The similar 391 ARGs were primarily genes determined to be the more abundant intrinsic ARGs in FPES 392 isolates whereas the distinct genes were often rare variants (i.e. only one copy across 393 isolate) within each database. When comparing by genus and database we found there 394 were significant differences in number of ARGs per isolate between Pseudomonas and 395 Bacillus isolates, which was higher in FPES for Pseudomonas and higher in RefSoil+ for 396 Bacillus ( Figure 5B). Across the 90 FPES isolates, plasmids contained 4 ARO types with 397 6 gene copies whereas there was no significant plasmid hits across all 127 RefSoil+ 398 isolates examined. The plasmid-borne ARGs from FPES included a BES-1 beta 399 lactamase, two antibiotic efflux pumps (TriC and KpnF), and an undecaprenyl 400 pyrophosphate related protein (bcrC) ( Table 2)