Genomic and metagenomic analysis of microbes in a soil environment affected by the 2011 Great East Japan Earthquake tsunami

Sample collection

Soil samplings were conducted at Hiyoriyama (38°15’20”N, 141°0’42”E) and Amamiya (38°16’35”N, 140°52’16”E) in Sendai city, Miyagi, Japan, in July 2012 (Fig. 1). If needed, the owners of the lands gave permission to conduct the study on these sites. We confirm that the study did not involve endangered or protected species. The Hiyoriyama site is 0.5 km off the coastline and was affected by the tsunami, whereas the Amamiya site is 12 km off the coastline and was not affected; the two sites are 13 km apart. The surface soil was removed to a 5 cm depth before sampling. Intermingled plants were carefully removed using tweezers, and soil that passed through a 2.0-mm pore-sized sieve was collected. The collected soil samples were transported to the laboratory at 4 °C and immediately stored at -80 °C until the subsequent analysis.

Seawater sampling was conducted at St.5 (38°06’00”N, 142°15’00”E) and St.6 (38°22’59”N, 142°43’01”E) off the coast of Miyagi, Japan in the Pacific Ocean, in August 2012, during the KT-12-21 cruise of R/V Tansei-Maru (JURCAOS, JAMSTEC). The St.5 and St.6 stations are located 110 km and 150 km from Sendai city, respectively. Surface seawater was collected in a prewashed bucket and immediately spread onto agar plates on a research vessel.

Isolation and 16S rRNA sequencing

R2A medium (Wako Pure Chemical Industries) was used to cultivate microbial strains that grow under general, low-nutrient condition, and ZoBell marine medium (Becton Dickinson and Company) was used to cultivate microbial strains that have adapted to a seawater-affected condition. Soil samples were thawed at 4 °C overnight, suspended in R2A or ZoBell liquid medium, and plated to the corresponding agar medium at a density of 10 ⁻⁴ g soil per plate with five replicates. The plates were incubated at 20 °C for 7 days before colony counting and picking. To obtain strains that could grow in both media, single colonies on the R2A agar plates were transferred to ZoBell agar plates with sterilized sticks, incubated at 20 °C for 7 days, and isolated by spread-plating on ZoBell agar at 20 °C. The seawater samples were plated to R2A and ZoBell agar at a volume of 100 μL seawater per plate with three replicates. The plates were incubated at 20 °C for 7 days before colony counting and picking.

To sequence 16S rRNA genes, seven strains that were isolated from the Hiyoriyama site and could actively grow in both R2A and ZoBell media were randomly selected. After incubation in ZoBell liquid medium, DNA was extracted using Wizard Genomic DNA Purification Kit (Promega). The 16S rRNA genes were amplified using a standard polymerase chain reaction protocol with the primers 27F (5’-AGAGTTTGATCMTGGCTCAG-3’) and 1492R (5’-GGCTACCTTGTTACGACTT-3’) [17], and sequenced by the Sanger method.

Whole-genome sequencing and analysis

Four strains of the Arthrobacter genus that were isolated from the Hiyoriyama site and were cultivable in both R2A and ZoBell media were targeted for whole-genome sequencing. Genomic DNA was extracted by the phenol-chloroform method. Two strains (named Hiyo4 and Hiyo8) were sequenced using PacBio RS II (Pacific Biosciences) according to the manufacturer protocols. De novo genome assembly of the 62,608 (Hiyo4) and 65,240 (Hiyo8) raw reads obtained from the Sprai pipeline (http://zombie.cb.k.u-tokyo.ac.jp/sprai/) successfully generated one and three circular contigs, respectively, after manual curation. The other two strains (Hiyo1 and Hiyo6) were sequenced using GS FLX+ System (Roche) and Ion PGM (Thermo Fisher Scientific) according to the manufacturer protocols. De novo genome assembly of the 301,881 (Hiyo1) and 267,295 (Hiyo6) reads obtained from the Newbler assembler [18] generated 38 and 630 scaffolds, respectively.

Coding sequences (CDSs) were predicted by applying Prodigal [19] to the contig sequences. Functional annotation was performed by blastp searches [20] against the Swiss-Prot [21] and eggNOG v4.0 [22] databases with a cut-off e-value ≤1E-5. Transfer RNA (tRNA) and rRNA sequences were predicted using tRNAscan-SE [23] and RNAmmer [24], respectively, with default settings.

For comparative genome analysis, all 21 publicly available genome sequences (6 complete and 15 draft sequences) of the Arthrobacter genus were downloaded from GenBank [25] via EzGenome (http://www.ezbiocloud.net/ezgenome) in January 2015 (Additional file 1). The CDSs of the four isolated and 21 downloaded genomes were subjected to blastp searches against the eggNOG database [22] with cut-off e-value ≤1E-5 and identity ≥90 %.

For construction of a phylogenetic tree, the 16S rRNA sequences of 56 Arthrobacter type strains and Streptomyces coelicoflavus NBRC 15399^T were additionally downloaded from the RDP webserver [26]. Streptomyces coelicoflavus NBRC 15399^T was used as an outgroup [27]. The 16S rRNA sequences of the total 82 strains were subjected to multiple alignment using MUSCLE [28] with default settings. A maximum-likelihood (ML) tree was generated by MEGA 6 [29] with the K80 substitution model with a gamma distribution and invariant sites (K2+G+I), which was the AIC-selected model, and 1000 bootstrap replicates. An ML tree of the total 17 genome-available strains was constructed on the basis of the set of 400 conserved bacterial marker genes using PhyloPhlAn [30] and MEGA 6 [29] with the WAG substitution model that incorporates gamma distribution and the amino-acid frequencies of the dataset (WAG+G+F), which was the AIC-selected model, and 1000 bootstrap replicates.

Culture assays of iron dependency

To determine the difference in iron tolerance among strains in relation to the genetic analysis results, culture assays were conducted at different iron concentrations. In addition to the four isolated Arthrobacter strains, we cultivated four closely related and genome-sequenced species, A. aurescens Phillips 1953^T (JCM 1330^T), A. chlorophenolicus A6^T (JCM 12360^T), A. globiformis Conn 1928^T (JCM 1332^T), and A. phenanthrenivorans Sphe3^T (JCM 16027^T). All four species had intact siderophore-synthesis genes in their genomes. These strains were provided by the Japan Collection of Microorganisms, BioResource Center, RIKEN and National BioResource Project of Ministry of Education, Culture, Sports, Science and Technology, Japan.

Iron-controlled, modified MM9 medium was prepared as follows. A solution containing 0.3 g/L KH₂PO₄, 0.5 g/L NaCl, 1.0 g/L NH₄Cl, 6.0 g/L NaOH, and 30.24 g/L PIPES was adjusted to pH 7.0 with NaOH. After autoclaving, separately sterilized solutions of 10 mL of 20 wt% glucose, 1 mL of 1 M MgCl₂, and 0.1 mL of 1 M CaCl₂ were added to 1 L of the solution [31]. Then, the iron concentration was adjusted to 0, 0.1, 1, or 10 μM with a FeCl₃-containing solution that was prepared in the same manner.

Each strain was precultured until its optical density at 660 nm (OD₆₆₀) reached 0.1 in the iron-free modified MM9 liquid medium. Then, 100 μL of the suspension was inoculated to 50-mL tubes containing 10 mL of the iron-controlled, modified MM9 medium. Among the additional four strains, only A. phenanthrenivorans Sphe3 showed growth in the modified MM9 medium. The tubes were incubated at 30 °C on a linear shaker at 200 rpm for 3 days, and the OD₆₆₀ was measured periodically during the incubation period. The growth curve was fitted to the logistic model to calculate the maximum growth rate.

Soil chemical analysis

The soil samples were subjected to chemical analysis for pH, electrical conductivity, and concentrations of total organic carbon, total nitrogen, nitrate, nitrite, ammonium, effective phosphate, exchangeable ions (K⁺, Ca²⁺, Mg²⁺, Na⁺, and Mn²⁺), available iron (Fe), chloride ion (Cl^-), sulfate ion (SO\(_{4}^{2-}\)), eluted heavy metals (Cd, Cr (VI), total Hg, alkyl mercury, Pb, As, and B), and contained heavy metals (Cd, Cr (VI), Hg, Pb, As, B, Cu, Zn, and Ni). The analysis was conducted by Createrra Inc. (Tokyo, Japan).

Shotgun metagenome sequencing and analysis

Metagenomic DNA was extracted using PowerSoil DNA Isolation Kit (MoBio Laboratories). Shotgun metagenome sequencing was performed using the GS FLX+ System according to the supplier’s protocol. Duplicated reads were removed by CD-HIT-454 [32].

Taxonomic assignment was performed using Kraken [33] against complete prokaryotic genomes from RefSeq [34]. CDS prediction was performed using MetaProdigal [35]. CDSs less than 30 amino acids in length were excluded from further analysis. Functional annotations were based on blastp searches against the eggNOG [22] and Swiss-Prot [21] databases with a cut-off e-value ≤1E-5.

SortMeRNA [36] was applied to the shotgun metagenome data to extract 16S rRNA sequences. For each extracted 16S rRNA sequence, a blastn search was performed against MetaMetaDB [37] and the top hit sequences with e-value ≤1E-10 and identity ≥90 % were retrieved. Microbial habitability index (MHI) scores were calculated as described previously [37].

Data deposition

The whole-genome and plasmid sequence data of Hiyo1, Hiyo4, Hiyo6, and Hiyo8 were deposited in the DDBJ/ENA/GenBank database under the BioSample ID SAMD00024042, SAMD00024043, SAMD00024044, and SAMD00024045, respectively. The shotgun metagenome sequence data of Hiyoriyama and Amamiya were deposited in the DDBJ/ENA/GenBank database under BioSample ID SAMD00023516 and SAMD00023517, respectively. All data were registered under BioProject ID PRJDB3373.

Isolation of microbial strains

To investigate whether the microbial community at the Hiyoriyama (tsunami-affected) site contained more microbes that are adapted to seawater-affected conditions than that at the Amamiya (tsunami-unaffected) site, we conducted culture experiments using R2A (general low-nutrient) and ZoBell (seawater-based) media. At Hiyoriyama, the mean (± standard deviation) numbers of colony forming unit (CFU) per gram of soil were 7.0 ± 3.9 ×10⁵ and 3.0 ± 2.0 ×10⁵ on R2A and ZoBell, respectively. At Amamiya, these numbers of CFU were 21.8 ± 4.7 ×10⁵ and 3.6 ± 2.3 ×10⁵. The ZoBell/R2A CFU ratios were 0.43 and 0.17 at Hiyoriyama and Amamiya, respectively, indicating that the Hiyoriyama site would be comparatively enriched with microbes adapted to a seawater-affected condition at 10 months after the tsunami. For comparison, surface seawater samples collected at St. 5 and St. 6 in the offshore were spread onto both agar plates. The numbers of CFU per milliliter of seawater were 12.7 ± 13.3 ×10¹ and 81.7 ± 43.6 ×10¹ on R2A and ZoBell, respectively. As expected, the ZoBell/R2A CFU ratio (6.4) was significantly higher at the offshore sites than at Amamiya and Hiyoriyama (p-value <0.05, t-test).

To isolate microbial strains that are potentially adapted to both types of environments from Hiyoriyama, we aseptically transferred the microbial colonies grown on R2A to ZoBell agar plates. Seven isolated colonies were randomly picked up and their 16S rRNA genes were sequenced. Unexpectedly, all of the strains were found to belong to a single genus, Arthrobacter. The genus Arthrobacter is an aerobic, gram-positive member of the family Micrococcaceae, Actinobacteria [38, 39]. This genus is broadly found in soils, as well as in extreme environments, including the deep subsurface [40], arctic ice [41], radioactive sites [42], and heavy metal-contaminated sites [43]. Some Arthrobacter species were reported to tolerate drastic environmental stresses, e.g., desiccation [44], starvation [45], heavy metals [46, 47], and radioactivity [48]. Furthermore, at the time of analysis, 6 complete and 15 draft genome sequences were available for comparative genome analysis. Because of these characteristics, the isolated Arthrobacter strains were targeted as a possible platform for exploring genomic features that may be related to microbial adaptation to drastically changed environments.

Whole-genome sequencing of the isolated Arthrobacter strains

The whole-genome sequences of four Arthrobacter sp. strains were determined (Table 1). Assembly using the reads from PacBio RS II showed the complete genomes of two strains: Hiyo4 with one circular chromosome (3.8 Mbp), and Hiyo8 with one circular chromosome (4.7 Mbp) and two plasmids (0.3 Mbp and 15 kbp) (Fig. 2). Assembly using reads from GS FLX+ and Ion PGM System produced 38 and 630 scaffolds for two strains, Hiyo1 and Hiyo6, respectively.

	Hiyo1	Hiyo6	Hiyo4	Hiyo8
Sequencing platform	GS FLX+ & Ion PGM	GS FLX+ & Ion PGM	PacBio RS II	PacBio RS II
Scaffolds	38	630	1	3
Contigs	1,685	2,450	1	3
Total genome size (bp)	5,543,883 a	2,594,729 a	3,779,248	4,698,617 a
N50 (bp)	4,950	2,656	-	-
Coverage	20x	24x	79x	42x
GC content (%)	63.2	63.3	65.0	63.8
CDSs	5,292	3,767	5,120	7,041
rRNAs	2	3	12	15
tRNAs	51	33	50	53

^a Plasmid sequences were not excluded

The total genome sizes of the four strains ranged from 2.6 to 5.5 Mbp. Although the genome sizes of Hiyo4 and Hiyo8 were within the range of the previously reported genomes (Additional file 1), their CDS numbers were exceptionally large, possibly because of additional genes that facilitate adaptation to different environmental conditions. The functional categories of eggNOG were assigned to 61–66 % of the CDSs, and the GC content was 63–65 %, which is similar to that of the previously reported genomes (59–67 %) (Additional file 1).

Phylogenetic analysis and comparative genomics

To reveal the phylogenetic relationships among the four strains Hiyo1, Hiyo4, Hiyo6, and Hiyo8, we constructed a maximum-likelihood phylogenetic tree of the Arthrobacter genus based on 16S rRNA gene sequences (Additional file 3). The tree reliably placed the four isolated strains within this genus. There was only one nucleotide base gap between the 16S rRNA sequences of Hiyo1 and Hiyo8, suggesting their close relationship. Hiyo4 and Hiyo6 were classified into different clades in the tree.

Subsequently, we conducted comparative genome analysis with 21 publicly available Arthrobacter genomes. The relative abundance of the CDSs assigned to each eggNOG functional category in each genome (Additional file 2) shows small difference among these Arthrobacter strains, i.e., their genomes have similar functional composition overall. The most striking difference between the Arthrobacter genomes isolated from the tsunami-affected soil and those isolated from other environments was that desferrioxamine B biosynthesis genes were independently lost in each of the former genomes. Within 14 publicly available, high-quality Arthrobacter genomes, the desferrioxamine B biosynthesis gene cluster and surrounding synteny structures were found to be highly conserved (Fig. 3). On the other hand, the desferrioxamine B biosynthesis gene cluster was entirely absent in the completed Hiyo1 and Hiyo8 genomes: the desA (pyridoxal-dependent decarboxylase) and desB (L-lysine 6-monooxygenase) genes of the cluster had nonsense mutations in the Hiyo4 genome, and neither the cluster nor the surrounding synteny structure was found in the Hiyo6 genome.

Desferrioxamine B is a member of the siderophores family of molecules, which are low-molecular-weight, iron-chelating compounds secreted by many microbes and plants for the uptake of iron [14–16, 49, 50]. The ability to use siderophores confers an ecological advantage when iron is limited [51]. Many Arthrobacter strains have a desferrioxamine B biosynthesis gene cluster, which is composed of four genes, named desABCD, for siderophore production [52, 53]. It should be noted that no iron-rich media were used during the isolation procedures.

The independent losses of the siderophore-synthesis genes are not likely to have occurred by chance but likely because of natural selection. Thus, these Arthrobacter strains were assumed to have been under weak selection pressure for iron uptake and to be at a growth disadvantage under low iron concentrations. To evaluate the growth potentials of these strains under various iron concentrations, culture experiments in iron-controlled media were conducted (Fig. 4; Additional file 4). Two of the isolated strains (Hiyo1 and Hiyo8) required 10 μM Fe³⁺ iron for rapid growth, whereas a control strain (A. phenanthrenivorans Sphe3) that has a desferrioxamine B biosynthesis gene cluster required 1 μM Fe³⁺ iron (Fig. 4). Notably, under the 1 μM Fe³⁺ iron concentrations, the maximum growth rates of Hiyo1 and Hiyo8 (1.08 ± 0.14 ×10⁻² and 1.15 ± 0.32 ×10⁻², respectively) was significantly smaller than that of Sphe3 (2.34 ± 0.20 ×10⁻²) (p-value <0.05, t-test with Bonferroni correction). Hiyo4 and Hiyo6 showed very weak growth even with 10 μM Fe³⁺ iron, possibly because these two strains require additional nutrients for growth (Additional file 4 C and D).

Based on these results, we hypothesized that strains with de novo mutations in siderophore synthesis genes or those that originally lacked these genes would be selected under iron-enriched conditions. Notably, siderophore production by soil-living microbes has been reported to help various plants absorb iron (e.g., tomato, cucumber, barley, and corn) [54–56] and has been associated with N₂ fixation (pigeon pea) [57]; therefore, these observed genomic changes in the bacterial communities might also relate to plant growth.

Soil chemical analysis

Except for these chemicals, the characteristics of the two samples were similar overall, suggesting that the two soil samples share a similar geological origin. In particular, the absence of heavy metals such as Pb, Hg, and Cu in Hiyoriyama might indicate that the soil was not completely covered or replaced with marine sediments. In addition, the similarities in electrical conductivity and Na⁺ and Cl^- content between samples can be attributed to the effects of rain; in the case of the 2004 Indian Ocean tsunami, water-soluble salts derived from the tsunami were strongly reduced after a rainy season in a coastal area in Thailand [64].

Shotgun metagenome sequencing

Using Kraken [33], 114,838 (14.0 %) and 112,459 (11.7 %) shotgun reads from Hiyoriyama and Amamiya were taxonomically classified, respectively. Almost all reads were assigned to Bacteria (99.01 and 98.99 %), whereas few reads were assigned to Archaea (0.81 and 0.63 %) and Viruses (0.18 and 0.38 %). The microbial composition of abundant genera is shown in Additional file 5. The most abundant genus in both samples was Burkholderia (4.85 and 7.20 %), followed by Bradyrhizobium (4.66 and 6.30 %), Rhodopseudomonas (3.27 and 3.46 %), and Pseudomonas (2.98 and 3.13 %). The similar composition of the major taxonomic groups reflects the similar overall chemical characteristics between the two soil samples. In addition, we estimated the typical habitats of the contained microbes by querying the extracted 16S rRNA genes against MetaMetaDB [37], a database that links 16S rRNA gene sequences to environments based on comprehensive analysis of published metagenomic and amplicon-sequencing datasets. The estimated habitats quantified as MHI values [37] showed that the top habitat was soil in both communities; however, the marine habitat was estimated to be modestly more abundant in Hiyoriyama, whereas the soil habitat was more abundant in Amamiya, as expected (Fig. 5).

Figure 6 displays the genera whose relative abundance substantially differed between the two samples, including only those whose abundance in one sample was more than three times greater than that in the other. Notably, Arthrobacter was the only genus that was both abundant in and differed substantialy between the two samples. Considering the fact that Arthrobacter was the genus cultivated in both the R2A and ZoBell media, we propose that this genus likely shows a greater potential for adaptation to tsunami-affected soils. The other genera that were more abundant in Hiyoriyama included Erysipelothrix, where all reads were assigned to a single species, Erysipelothrix rhusiopathiae [65], which is known to cause erysipelas, a bacterial skin infection, in animals [66]. Although previous culture-based studies reported several pathogen species (Mycobacterium elephantis, Massilia timonae, Vibrio ichthyoenteri, V. natriegens, and V. fluvialis) in sludge derived from tsunami-affected soil in Tohoku [13, 67], these species were not detected in the present dataset. It may also be notable that the genera detected only in Hiyoriyama included typical marine-living groups such as Croceibacter [68], Marinitoga [69], and Pyrococcus [70–72], implying that the tsunami facilitated microbial immigration.

We annotated the CDSs and investigated the relative abundance of nitrogen cycle-related genes, because the taxonomic analysis identified genera known to metabolize inorganic nitrogens, such as Bradyrhizobium, Azospirillum, Frankia, Mesorhizobium, Rhizobium, and Sinorhizobium (Additional file 5), and the chemical analysis revealed differences in the amount of nitrogen compounds (Table 2). The abundance of functional genes showed that genes related to denitrification and nitrogen fixation were more abundant in Hiyoriyama and genes related to nitrite reduction were more abundant in Amamiya (Fig. 7). In addition to the oxidization of iron-sulfur compounds, this dominance of denitrification-related genes at Hiyoriyama may be another cause of the relatively small amount of nitrate observed in Hiyoriyama (Table 2), which might affect terrestrial vegetation indirectly.

We also investigated the abundance of siderophore-synthesis genes in the shotgun metagenome data, but only three and four reads of genes that are involved in this process (bactNOG07545, bactNOG14638, bactNOG30540) were detected in Hiyoriyama and Amamiya, respectively, which is not a sufficient sample size for statistical analysis. Differences in sulfur metabolism genes were not as large as those of nitrogen metabolism genes. A substantial difference was observed in the numbers of cation transporter genes, where 127 and 59 monovalent cation/H⁺ antiporter subunits, and 96 and 16 Na⁺/Ca²⁺ antiporter family proteins (bactNOG00892) were detected in Hiyoriyama and Amamiya, respectively. These genes may have facilitated salt tolerance in the tsunami-affected soil, because cation transporters are known to function in bacterial salt tolerance [73, 74].