A Nitrate-Transforming Bacterial Community Dominates in the Miscanthus Rhizosphere on Nitrogen-Deficient Volcanic Deposits of Miyake-jima

The perennial gramineous grass Miscanthus condensatus functions as a major pioneer plant in colonizing acidic volcanic deposits on Miyake-jima, Japan, despite a lack of nitrogen nutrients. The nitrogen cycle in the rhizosphere is important for the vigorous growth of M. condensatus in this unfavorable environment. In the present study, we identified the nitrogen-cycling bacterial community in the M. condensatus rhizosphere on these volcanic deposits using a combination of metagenomics and culture-based analyses. Our results showed a large number of functional genes related to denitrification and dissimilatory nitrate reduction to ammonium (DNRA) in the rhizosphere, indicating that nitrate-transforming bacteria dominated the rhizosphere biome. Furthermore, nitrite reductase genes (i.e., nirK and nirS) related to the denitrification and those genes related to DNRA (i.e., nirB and nrfA) were mainly annotated to the classes Alpha-proteobacteria, Beta-proteobacteria, and Gamma-proteobacteria. A total of 304 nitrate-succinate-stimulated isolates were obtained from the M. condensatus rhizosphere and were classified into 34 operational taxonomic units according to amplified 16S rRNA gene restriction fragment pattern analysis. Additionally, two strains belonging to the genus Cupriavidus in the class Beta-proteobacteria showed a high in vitro denitrifying activity; however, metagenomic results indicated that the DNRA-related rhizobacteria appeared to take a major role in the nitrogen cycle of the M. condensatus rhizosphere in recent Miyake-jima volcanic deposits. This study elucidates the association between the Miscanthus rhizosphere and the nitrate-reducing bacterial community on newly placed volcanic deposits, which furthers our understanding of the transformation of nitrogen nutrition involved in the early development of vegetation.


Introduction
The soil-plant ecosystem is the main reservoir of elements such as carbon and nitrogen, and the main players in the element/material cycle are plants and microorganisms. Although carbon and nitrogen are not contained in the soil parent materials, these elements are introduced by physical deposition from the atmosphere and biological fixation. Primary succession on glacier forelands [1] and volcanoes [2] presents an ideal opportunity to understand the formation of terrestrial ecosystems through biological colonization.
Miyake Island (Miyake-jima) in Japan is an active volcanic island located on the western rim of the Pacific Ocean, with a humid subtropical climate. Since 1084, there have been fifteen eruptions on Miyake-jima, with the most recent eruptions occurring in 1874,1940,1962,1983, and July-August 2000 [3]. The effects on the ecosystem of the last eruption involved the formation of a new caldera and the deposition of ash to form acidic volcanic deposits (pH 3-4) that lack carbon and nitrogen but have high amounts of either exchangeable cation calcium or aluminum. The area around the crater was exposed to poisonous gases such as sulfur dioxide and hydrogen sulfide, which caused extensive damage to the vegetation in the area. After the eruption in 2000, almost all vegetation that covered around the crater was lost. Sato et al., (2009) clarified that pioneer microorganisms of volcanic ash deposits (pH 3) acidified with volcanic gas are chemotrophs near the mountaintop of Oyama, Miyakejima [4]. Fujimura et al., (2012 and2016) revealed the role of early microbial communities in supporting pioneering plant colonization in unvegetated volcanic deposits, suggesting that the bacterial community structure transitioned from chemotrophs (autotrophs) to heterotrophs [2,5]. In addition, Guo et al. (2014) and Lathifah et al. (2019) reported that successional change in the early soil microbial community occurred with primary vegetation that developed on the newly placed volcanic deposits, in which the plant benefit groups are mainly affiliated with the orders Burkholderiales and Rhizobiales, which are significantly increased with the vegetation succession [6,7]. Guo et al. (2021) further revealed that different plant benefit groups dominated the rhizosphere in different developmental phases on the recent Miyake-jima volcanic deposits, i.e., Paraburkholderia and Trinickia represented the major plant-beneficial groups in the early colonization phase, whereas Rhizobiales and Arthrobacter increased in the later colonization phase [8]. Miscanthus grass is considered tolerant to high Al 3+ concentrations and acidic conditions, as well as nitrogen deficiency [9], whereas associated arbuscular mycorrhizal fungi have been reported to support Miscanthus colonization on acidic soils [10,11]. In addition to the dominance of nitrogen-fixing bacteria, a considerable proportion of denitrification-related genes (i.e., nirK, nirS, norB, and nosZ) were detected as increasing in the 3.5-9.5-year-old Miyake-jima volcanic deposits at an unvegetated study site [2]. However, those genes associated with a competitive process in soils, i.e., dissimilatory nitrate reduction to ammonium (DNRA), have not yet been focused.
Denitrification and DNRA are thought of as the major microbial processes of nitrate transformation in terrestrial ecosystems [12]. The former induces gaseous losses of nitrogen from soils via the reduction of nitrate (NO 3 − ) to nitrous oxide (N 2 O) and/or dinitrogen (N 2 ) [13][14][15], whereas the latter reduces NO 3 − to ammonium (NH 4 + ), protecting NO 3 − from leaching and gaseous losses [16] and enriches soils with readily available NH 4 + for plants, such as N-fertilizer [12]. Therefore, the investigation of denitrification-and DNRAassociated genes and communities in the rhizosphere is expected to elucidate the supply of N-nutrition to pioneer plants and gain a comprehensive understanding of vegetation recovery on barren land.
Recent advances in high-throughput sequencing technologies revealed that nitrogencycling microbes are unexpectedly diverse in their functions and phylogenies [17]. Vegetation recovery after the volcanic eruption on Miyake-jima presents an ideal opportunity to increase understanding of the formation of the terrestrial ecosystem. However, there is still lack of information on the nitrogen-cycling microbial community associated with the pioneer plant Miscanthus, particularly, those associated with the processes of denitrification and DNRA. The objective of this study was to characterize and identify the nitrogen-cycling bacterial community associated with the Miscanthus rhizosphere on Miyake-jima volcanic deposits based on shotgun metagenomic and culture-based analyses.

Site Description and Sample Collection
Miyake-jima, designated as one of national parks of Japan (https://www.kankyo. metro.tokyo.lg.jp/naturepark/english/know/park/introduction/kokuritsu/fujihakone/ miyakejima/index.html, accessed on 24 November 2022), provides a natural ecological model to study the ecosystem restoration after the volcanic eruption. We were authorized to collect the Miyake-jima soil, volcanic deposit, and plant samples as the minimum required by this investigation; therefore, we collected one replicate of each sample for the sampling events in May 2016, September 2017, and September 2018. Site IG7 was established on the windward (north) side of Mt. Oyama (34 • 05.37 N and 139 • 30.84 E) of Miyake-jima at an altitude of 547 m a.s.l. (Figure 1). The climate type in this region is humid subtropical, with a mean annual rainfall of 3024.7 mm and a mean air temperature of 18.0 • C from 1991 to 2020, according to Japan Meteorological Agency (https://www.data.jma.go.jp/obd/stats/etrn/ view/nml_sfc_ym.php?prec_no=44&block_no=47677&year=&month=&day=&view=, accessed on 24 November 2022). Vegetation has gradually recovered halfway up the mountaintop since 2001 with the colonization of Miscanthus condensatus, belonging to the family Poaceae as one of the pioneer grasses ( Figure 1A-C) [18][19][20][21]. At site IG7, the surface coverage of M. condensatus increased as follows: <20% in 2009, <45% in 2013, and 90% in 2018. Complete Miscanthus grass samples (approximately 70 cm in height and 12 cm width at the root rhizosphere) were collected in September 2017 and September 2018. In March 2016, only volcanic deposit layer (referred to as the C layer), in 2017, a new organic soil layer (referred to as the A layer) on surface and C layer (>20 cm in depth), and in 2018, the A layer (0-10 cm in depth) and the upper (10-20 cm in depth) and lower (20-30 cm in depth) volcanic deposit layers (referred to as the C1 and C2 layer, respectively) were observed and randomly sampled from several plots and mixed in sterile plastic bags. All samples were kept at 4 • C and transported to laboratory within 2 days.

Site Description and Sample Collection
Miyake-jima, designated as one of national parks of Japan (https://www.kankyo.metro.tokyo.lg.jp/naturepark/english/know/park/introduction/kokuritsu/fujihakone/miyakejima/index.html, accessed on 24 November 2022), provides a natural ecological model to study the ecosystem restoration after the volcanic eruption. We were authorized to collect the Miyake-jima soil, volcanic deposit, and plant samples as the minimum required by this investigation; therefore, we collected one replicate of each sample for the sampling events in May 2016, September 2017, and September 2018. Site IG7 was established on the windward (north) side of Mt. Oyama (34° 05. 37' N and 139° 30. 84' E) of Miyake-jima at an altitude of 547 m a.s.l. (Figure 1). The climate type in this region is humid subtropical, with a mean annual rainfall of 3024.7 mm and a mean air temperature of 18.0 °C from 1991 to 2020, according to Japan Meteorological Agency (https://www.data.jma.go.jp/obd/stats/etrn/view/nml_sfc_ym.php?prec_no=44&block_n o=47677&year=&month=&day=&view=, accessed on 24 November 2022). Vegetation has gradually recovered halfway up the mountaintop since 2001 with the colonization of Miscanthus condensatus, belonging to the family Poaceae as one of the pioneer grasses ( Figure  1A-C) [18][19][20][21]. At site IG7, the surface coverage of M. condensatus increased as follows: <20% in 2009, <45% in 2013, and 90% in 2018. Complete Miscanthus grass samples (approximately 70 cm in height and 12 cm width at the root rhizosphere) were collected in September 2017 and September 2018. In March 2016, only volcanic deposit layer (referred to as the C layer), in 2017, a new organic soil layer (referred to as the A layer) on surface and C layer (>20 cm in depth), and in 2018, the A layer (0-10 cm in depth) and the upper (10-20 cm in depth) and lower (20-30 cm in depth) volcanic deposit layers (referred to as the C1 and C2 layer, respectively) were observed and randomly sampled from several plots and mixed in sterile plastic bags. All samples were kept at 4 °C and transported to laboratory within 2 days.

Soil Chemical Properties Analysis
A portion of C, C1, C2, and A samples (4-8 g in weight) was dried at 105 • C for over 12 h. Total carbon and total nitrogen in the dried samples were determined using a JM3000 CN counter with a JMA 3000 Auto Samper (J-Science Lab., Kyoto, Japan). The volumetric water content was calculated as reported previously [6]. A slurry consisting of a 1:2.5 mass ratio of wet sample and distilled water was used to determine the sample pH value. Fresh C1, C2, and A samples (5-15 g in weight) were used to determine nitrate (NO 3 − ) and ammonium (NH 4 + ) contents using an SPCA-6210 (Shimadzu, Kyoto, Japan). The measurement was initiated by 1 N KCl extract and mixed solution sample with either distilled water for NO 3 − or a combination of solution A (phenol, 20% NaOH/sodium hydroxide, acetone) and solution B (NaClO/sodium hypochlorite) for NH 4 + . All analyses were performed in triplicates except for water content measurements.

Sample Pretreatment of Miscanthus
The root-rhizosphere of Miscanthus grass at the IG7 site, which was sampled in 2017 and 2018, was subjected to pretreatment for metagenomic and culture-dependent analyses. To obtain the rhizosphere soil, we used soil-associated roots (tiny soil) with consideration for the root growth features of Miscanthus, which forms a net-like structure and binds loose soil. Due to the one plant colony that was obtained for each sampling event, the rhizosphere soils were recovered from the plant roots as fully as possible and were pooled together to reduce bias. The Miscanthus roots were cut and separated using sterilized scissors in sterile conditions. Approximately 3 g of root sample was transferred to 30 mL of sterile distilled water in a 50 mL Falcon tube [8]. In order to leach the rhizosphere soil from roots, the samples were sonicated using an EYELA ultrasonic cleaner machine (Tokyo Rikakikai, Kanagawa, Japan) for about 5 min. This operation was repeated several times as necessary. The suspension of rhizosphere soil was centrifuged at 10,000× g for 5 min to form pellets. After removing the supernatant, the pellets were treated as rhizosphere samples. Several preparations of rhizosphere soil pellets were pooled together and stored at −20 • C until DNA extraction for shotgun metagenomic analysis. In addition, the main and lateral roots were separated as shown in Figure S1, and then the adhered soil was recovered in sterile distilled water in accordance with the above-mentioned method. The suspensions of rhizosphere soil recovered from the main and lateral roots were stored at 4 • C and used for the bacterial isolation within 1 day.

Genome Sequencing, Assembly, and Annotation of Metagenomic Contigs
The extraction of total DNA from 0.75 g rhizosphere sample of Miscanthus was performed using a DNeasy PowerSoil kit (QIAGEN, Hilden, Germany) in accordance with the manufacturer's instructions. The obtained DNA (16.2 µg in total) was used to build a paired-end library with insert size of~550 bp sequenced on an Illumina HiSeq2500 platform at GeneBay, Yokohama, Japan (http://genebay.co.jp, accessed on 24 November 2022). The obtained reads were analyzed using a metagenomic data analyzing pipeline metaWRAP version 1.0 (https://github.com/bxlab/metaWRAP/releases/tag/v1.0, accessed on 24 November 2022) [22]. Briefly, the raw reads were trimmed using the metaWRAP Read_qc module with default parameters and were assembled using the metaWRAP Assembly module with the megahit option and default settings [23]. All the outputted contigs were used to identify their taxonomy using the metaWRAP Kraken module and to calculate their abundance in each sample using the metaWRAP Blobology module with default parameters. The ribosomal RNA genes and protein-coding sequences (CDSs) were predicted using Barrnap version 0.9 (https://github.com/tseemann/barrnap, accessed on 24 November 2022) and MetaGeneMark (http://exon.gatech.edu/meta_gmhmmp.cgi, accessed on 24 November 2022) [24], respectively. The gene function and taxonomic identification of the predicted CDSs was assigned using GhostKOALA version 2.2 with the KEGG database [25], and the COG categories were assigned using RPS-BLAST against the COG database [26]. To calculate the relative abundances of the CDSs, the reads outputted from the metaWRAP Read_qc module were mapped to contigs using the BWA version 0.7.17 options of "index" and "mem" with the default parameters, and then the resulting SAM file was converted to BAM format and sorted using the Samtools version 1.9 options of "view" and "sort" [27].
Finally, the relative abundances of the CDSs were calculated as gene copies per million mapped reads (CPM) using the tool "jgi_summarize_bam_contig_depths" of MetaBAT version 2.12.1 [28].

Isolation of the Miscanthus Root-Associated Bacteria
The soil rhizosphere suspension (1.0 mL) was transferred to 9.0 mL of sterile distilled water to make 10 −2 , 10 −3 , and 10 −4 dilutions [6]. Then, 100 µL of each diluted soil suspension was spread onto 1.5% agar plates of 100-fold diluted nutrient broth medium supplemented with 3.0 mM sodium nitrate (NaNO 3 ) and 4.4 mM sodium succinate (DNB-NS medium and pH 6.8) in triplicates. The plates were grown both aerobically and anaerobically at 30 • C for 2 weeks [29,30]. Additionally, to retard bacterial growth, DNB-NS medium supplemented with antibiotic nalidixic acid was also used in this study and grown anaerobically at 30 • C for 2 weeks [31]. Single colonies appearing on the 10 −3 and 10 −4 plates were selected according to unique morphotypes [32], and were subcultured repeatedly on the same media with the same incubation conditions in order to obtain pure cultures of isolates.

Terminal-Restriction Fragment Length Polymorphism (T-RFLP) Profiling
A portion of the 10 −2 plates containing grown isolates as described above was used in this analysis. Total genomic DNA was extracted from the bacterial cells grown on the plates using ISOIL for beads beating kit (Nippon Gene, Tokyo, Japan) in accordance with the manufacturer's procedures. Polymerase chain reaction (PCR) amplification of the bacterial 16S rRNA gene was performed as follows: the total volume of the PCR mixture was 30 µL, which contained 1 µL of template DNA (100 ng/µL), 0.5 µL of 10 µM forward primer (Q-10F, 5 -CAGTTTGATCCTGGCTCAG-3 ; J-Bio21, Kisarazu, Japan), 0.5 µL of 10 µM reverse primer (926r, 5 -CCGTCAATTCCTTTRAGTTT-3 ) [4], 3 µL of 10 × Ex Taq buffer (20 mM Mg 2+ plus), 3 µL of dNTP mixture (2.5 mM each), and TaKaRa Ex Taq polymerase (TaKaRa Bio, Shiga, Japan). Q-10F primer labelled with a quenching fluorescence was purchased from J-Bio21. PCR conditions were initially 95 • C (2 min), then 25 cycles of 95 • C (30 s), 54 • C (45 s), and 72 • C (90 s). PCR reactions were carried out in a TaKaRa PCR Thermal Cycler Dice ® Touch (TaKaRa Bio). The PCR products were purified using a QIAquick PCR purification kit (QIAGEN) in accordance with the manufacturer's protocol. AluI, HaeIII, HhaI, and MspI (TaKaRa Bio) were used for restriction enzyme digestion, and the fragments were purified using a QIAquick PCR purification kit (QIAGEN) [33]. Two microliters of terminal-restriction fragments (T-RFs) were mixed with 12 µL of Hi-Di™ Formamide Genetic-analysis grade (Applied Biosystems, Foster City, CA, USA) and 0.1 µL of GeneScan™ 500 Liz ® Size Standard. T-RFs were denatured at 96 • C for 2 min and quenched thereafter. The length of T-RFs was determined using an ABI PRISM 3130xl Genetic Analyzer (Applied Biosystem) in GeneMapper mode [33]. Analysis for the determination of T-RF lengths was performed by using GeneMapper ® Software version 3.7 (Applied Biosystems, Foster City, CA, USA). T-RF data were developed using T-REX software [34]. T-RFs whose sizes were under 50 bases and <0.5% of the total area were omitted from this analysis. T-RFs with a different size that was less than 0.5 bp, digested with the same enzyme, and obtained from the same sample, were grouped as one T-RF. Based on the T-RFs profile, a Bray-Curtis-based distance matrix of the cultivated bacterial communities was calculated using the vegan package of R software version 3.5.3 and visualized using a wild bootstrapping dendrogram [35].

Amplified Ribosomal DNA Restriction Analysis (ARDRA)
Direct PCR of intact bacteria (colony PCR) method was used to amplify the 16S rRNA gene of each obtained strain (i.e., a pure culture) using the 10F and 1541R primer pair (10 µM) [4]. The reaction mixture (10 µL) contained 5 µL of GoTaq Green Master Mix (Promega, Madison, WI, USA), 0.5 µL of each primer, 4 µL of sterile milliQ water, and a selected bacterial colony. Prior to PCR, heat treatment was performed to lyse the cells using TaKaRa PCR Thermal Cycler Dice ® Touch, with the following conditions: denaturation at 95-96 • C for 8 min. Amplification of the bacterial 16S rRNA gene was carried out with the following conditions: initially 95-96 • C (8 min), followed by 25 cycles of 95 • C (30 s), 55 • C (60 s), and 72 • C (120 s). After the final extension at 72 • C for 7 min, the PCR mixtures were kept at 4 • C. The PCR products were digested directly with HinfI (50 U/µL, TaKaRa Bio) by mixing 4.5 µL of PCR product with 1 µL of restriction enzyme. The reaction mixtures (5.5 µL) were incubated at 37 • C for 1 h and electrophoresed using 1.5% agarose LO3 (TaKaRa Bio) with TAE (Tris-Acetate-EDTA) buffer to visualize the restricted fragment (RF) pattern. Strains with an identical RF profile were grouped as an operational taxonomic unit (OTU) through manual assessment according to a previous report [36]. The relative abundance of each OTU was determined by dividing the OTU number by the total number of OTUs. The heatmap of OTU abundance in total was visualized based on log 2 transformations of the OTU table using the formula log 2 (10 × x + 1), where x is the relative abundance of each OTU.

Identification of Representative Strains of Each OTU Based on the 16S rRNA Gene
Genomic DNA was extracted from the bacterial cells using a phenol-chloroform extraction method [4]. Amplification of the 16S rRNA gene was performed using the 10F and 1541R primer pair (10 µM) [4]. The total volume of the PCR mixture was 50 µL, which contained 1 µL of template DNA (100 ng/µL), 2.5 µL of 10 µM primers, 5 µL of 10 × Ex Taq buffer (20 mM Mg 2+ plus), 4 µL of dNTP mixture (2.5 mM each), and TaKaRa Ex Taq polymerase (Takara Bio, Otsu, Japan). PCR reactions were carried out in a TaKaRa PCR Thermal Cycler Dice ® Touch with the following conditions: initially 95 • C (5 min), then followed by 25 cycles of 95 • C (30 s), 60 • C (1 min), and 72 • C (1 min). PCR mixtures were kept at 4 • C. To sequence the purified PCR amplicons, the reaction mixture (10 µL) containing 1.5 µL of ×5 sequence buffer (Thermo Fisher Scientific, Waltham, MA, USA), 1 µL of each primer, 1 µL of DNA template, 5.5 µL of sterile distilled milliQ water, and 1 µL of BigDye terminator™ v3.1 (Thermo Fisher Scientific) was performed in similar conditions as described above. After purification using the ethanol precipitation method, the amplicons were subjected to DNA sequencing using an ABI PRISM 3130xl Genetic Analyzer. Additionally, the shotgun metagenomic reads were mapped to the 16S rRNA gene sequences of representative isolates for the ARDRA-based OTUs, and then the relative abundances of the OTUs were calculated as gene copies per 1 million mapped reads (CPM), using the above-mentioned methods.

Measurement of the Potential Activity of Denitrification
The potential denitrification activities of each strain were determined using the acetylene-block method, as described previously [29,30]. Isolates were incubated in 5 mL of liquid DNB-NS medium in aerobic conditions for 3 weeks. After that, 30 µL of liquid culture (optical density measured at a wavelength of 600 nm was set at 0.7) was transferred into 20 mL vials containing 3 mL of liquid DNB-NS medium and covered by a rubber cap. The headspace air was replaced with dinitrogen (N 2 ) gas using a Replace Injector GR-16 (Sanshin, Yokohama, Japan) followed by 2 mL Ar-C 2 H 2 (90:10) gas injection into vials using a sterile 5 mL syringe. After incubation at 28 • C for 2 weeks, 1 mL of the headspace gas was transferred into a new 20 mL vial containing N 2 gas. Thereafter, about 1 mL of each diluted sample was analyzed for N 2 O concentration using a gas chromatography GC-14A machine (electron capture detector-ECD type, Shimadzu). The quantity of water dissolved N 2 O was calculated as Bunsen absorption coefficient, then strains reducing ≥20% of added NO 3 − to N 2 O were considered as denitrifying bacteria [37].

Phylogenetic Analysis
To assign the taxonomy of strains at the genus level, we filtered and trimmed the partial 16S rRNA gene sequences using seqtrace version 0.9.0 and searched the sequences using RDP seqmatch version 3 against RDP database (release11_6) [38]. All reference sequences and the related out sequences were aligned to construct a phylogenetic tree based on the maximum likelihood method using a Tamura 3-parameter model with 1000 bootstrap replicates within MEGA7 software (https://www.megasoftware.net, accessed on 24 November 2022) [39].

Accession Number
Metagenomic raw sequence data obtained in this study have been submitted to DDBJ/NCBI/EMBL under accession number DRA010653. The metagenome-assembled contigs, as well as the predicted CDSs with functional annotations and relative abundance have been deposited in the scientific repository Figshare with the DOI number 10.6084/m9.figshare.19328753. The partial 16S rRNA gene sequences obtained in the present study have been deposited in the DDBJ database with the following accession numbers: LC699492-LC699525.

Chemical Properties of Substrates in Site IG7
In the 2016 investigation of site IG7, a new organic soil layer (the A layer) approximately 10 cm in thickness was observed, and in the 2018 investigation, the A layer buried the volcanic deposits (the C1 and C2 layers) ( Figure 1D). The concentrations of total carbon and total nitrogen, as well as nitrate (NO 3 − ) and ammonium (NH 4 + ) in soil post-eruption, were approximately 3-to 40-fold higher than those in volcanic deposits that contained the least carbon and nitrogen sources (Table 1). In addition, these samples showed high carbon/nitrogen (C/N) ratios of 40.5 in 2017 and 19.8 in 2018, respectively ( Table 1). The total carbon content in 2017 and 2018 had increased, which appeared to be related to the development of Miscanthus. Although there was no remarkable difference in NO 3 − concentration on the surface at the IG7 site between March 2016 and September 2018, the concentration of NH 4 + in September 2018 apparently increased compared with that in 2016. The least water content in volcanic deposits was also found, and all substrate samples were exposed to acidic conditions (Table 1). The measurement was performed three times and the average value was showed. A, C1, and C2 layers indicate the soil post-eruption, latter volcanic deposits, and former volcanic, respectively. n.a., not available; n.d., No data; N.D., not detected.

Metagenome Assembly and Nitrogen Cycle-Related Gene Abundance
A total of 24.7 million high-quality paired-end reads (2 × 150 bp) was assembled into 1,194,015 contigs with a total length, maximum length, and N 50 value of 710,571,951 bp, 264,351 bp, and 569 bp, respectively. According to the taxonomic identification by Kraken, 541,517 bacterial, 4890 eukaryotic, 623 archaeal, and 381 viral contigs were identified for the Miscanthus rhizosphere microbiome, whereas the rest of the contigs were not successfully assigned to any known taxonomies. Furthermore, 1,499,158 CDSs were predicted from the metagenome assembly, and 71% of the predicted CDSs were assigned to COG categories.
DNA shotgun sequence analysis was carried out to identify genes that are directly influenced by root activities, preserved as the primary site for plant-microbe interaction, and tolerant of low nutrient conditions in the Miscanthus rhizosphere. We identified genes encoding the enzymes/proteins involved in chemotaxis and secretion systems, stress response sigma factors, transporters, and a mediator of the stringent response that coordinates a variety of cellular activities in response to changes in nutritional abundance ( Table 2). On the other hand, nitrogen cycle-related genes encoding the fixation, reduction, and oxidation enzymes based on the KEGG nitrogen metabolic pathway comprised 0.11% (1662 CDSs) of the total predicted genes. Functional annotation of nitrogen cyclerelated genes revealed highly consistent numbers for denitrification (nar, nap, nir, nor, and nos), nitrogen fixation (nif ), respiratory DNRA (nrfA), and fermentative DNRA (nirB and nirD) (Figure 2A). Further analysis of the denitrifying bacterial community showed that nirK (68 CDSs) and nirS (6 CDSs) were mainly conserved in the genomes of members in the classes Alpha-, Beta-and Gamma-proteobacteria ( Figure 2B). The specific key gene nrfA (19 CDSs) of the respiratory DNRA bacteria was mainly close to those possessed by Gamma-proteobacteria, whereas the gene nirB (547 CDSs) from fermentative DNRA bacteria was associated with a more diverse taxa including Acidobacteria, Planctomycete, Alpha-, Beta-, and Gamma-proteobacteria ( Figure 2B). In addition, 19 CDSs encoding NifH (total abundance: 26.5 copies per 1 million reads) were detected and mainly identified in the class Gamma-proteobacteria ( Figure 2B). Other nitrogen-cycling-related genes such as nitrification, assimilatory nitrate reduction to ammonium, anaerobic ammonium oxidation (anammox), and complete ammonia oxidation, were rarely found.   (68 CDSs) and nirS (6 CDSs) were mainly conserved in the genomes of members in the classes Alpha-, Beta-and Gamma-proteobacteria ( Figure 2B). The specific key gene nrfA (19 CDSs) of the respiratory DNRA bacteria was mainly close to those possessed by Gamma-proteobacteria, whereas the gene nirB (547 CDSs) from fermentative DNRA bacteria was associated with a more diverse taxa including Acidobacteria, Planctomycete, Alpha-, Beta-, and Gamma-proteobacteria ( Figure 2B). In addition, 19 CDSs encoding NifH (total abundance: 26.5 copies per 1 million reads) were detected and mainly identified in the class Gamma-proteobacteria ( Figure 2B). Other nitrogen-cycling-related genes such as nitrification, assimilatory nitrate reduction to ammonium, anaerobic ammonium oxidation (anammox), and complete ammonia oxidation, were rarely found. The key enzymes involved in the denitrification process (Nar, Nap, Nir, Nor, and Nos); nitrogen fixation (Nif); dissimilatory nitrate reduction to ammonium (DNRA; NrfA, NirB, and NirD); nitrification process (Hao, and Amo); assimilatory nitrate reduction (Nas); anaerobic ammonium oxidation (anammox; Hzo, and Hzs) are selected to show the nitrogen-cycling process. (B), heatmap shows taxonomic profile of the key genes involved in the nitrogen-cycling process.

Diversity and Taxonomy of Cultured Bacteria with Potential Nitrogen Transformation in the Miscanthus Rhizosphere
The diversity of the culturable bacteria was determined by T-RFLP profiling using four different restriction enzymes. According to the T-RF profiles, the α-diversity index of each sample was not significantly different, suggesting less impact from different culture conditions as well as root types (Table 3). However, by adding an antibiotic, nalidixic acid in anaerobic conditions, each diversity index of culturable bacteria was slightly reduced compared with those indices in anaerobic conditions without the antibiotics (Table 3). In addition, the culturable bacterial communities in different culture conditions were grouped into two clusters (Supplementary Figure S2). One cluster contained those isolated from the rhizosphere of the main and lateral roots in aerobic conditions and from those of the main root in anaerobic condition with antibiotics, whereas the communities incubated in the other anaerobic conditions were clustered together. A total of 304 isolates were obtained in six culture conditions and assigned to 34 OTUs according to the RF patterns of each strain generated by ARDRA, and a representative strain was randomly selected from each OTU to identify the taxonomy (Figure 3 and Table 4). Nucleotide sequences of 16S rRNA genes were determined and designated as "IG7-n", where n is the OTU number (Table 4). In addition, the metagenomic reads were mapped to the 16S rRNA gene sequences to deduce the relative abundance of OTUs and/or strains in the rhizosphere community ( Figure 3). These results showed that the isolated OTUs were affiliated to four phyla (Proteobacteria, Actinobacteria, Firmicutes, and Bacteroidetes), six classes (Alpha-proteobacteria, Beta-proteobacteria, Gamma-proteobacteria, Actinobacteria, Bacilli, and Flavobacteria), and eighteen genera (Agrobacterium, Rhizobium, Starkeya, Bosea, Duganella, Achromobacter, Alcaligenes, Paraburkholderia, Cupriavidus, Mitsuaria, Rahnella, Klebsiella, Stenotrophomonas, Pseudomonas, Mycobacterium, Leifsonia, Bacillus, and Flavobacterium) ( Table 4 and Figure 3).

Denitrifying Activity
In total, 34 representative strains of the ARDRA-based OTUs were determined for their potential denitrifying activity using the acetylene block method [29]. Five representative strains showed denitrifying activity that reduced more than 20% of the supplied nitrate to nitrous oxide ( Figure 4A), and these strains were characterized as the denitrifying bacteria in accordance with Tiedje (1994) [37]. Strains IG7-14 and IG7-17 identified as the genus Cupriavidus within the class Beta-proteobacteria ( Figure S3A), which were isolated from the main roots in anaerobic conditions and reduced 57.3% and 87.4% of the supplied nitrate to nitrous oxide, respectively (Figure 4). Strain IG7-9, identified as the genus Agrobacterium within the class Alpha-proteobacteria ( Figure S3B), and strains IG7-30 and IG7-31 affiliated with the genus Flavobacterium within the class Flavobacteria ( Figure S3C), which were obtained in aerobic conditions, were found to have reduced 20-30% of the suppled nitrate to nitrous oxide ( Figure 4A). Representative strains affiliated with the classes of Gammaproteobacteria, Actinobacteria, and Bacilli were able to reduce less than 9.4% of the suppled nitrite to nitrous oxide (Figure 4). shows strains that reduced more than 5% of the supplied NO3 − to N2O. Black circles and triangles indicate strains that produced N2O from the supplied NO3 − at levels of <20% and >20%, respectively. (B) shows strains that reduced less than 5% of the supplied NO3 − to N2O. The x-axis and y-axis show the denitrifying activity (%) and the class of OTUs, respectively.

Discussion
Previous studies have revealed the role of Miscanthus as pioneer plants in disturbed and acidic soil as well as their ability to form root-symbiotic relationships with diazotrophic bacteria and arbuscular mycorrhizal fungi to promote their growth in unfavorable environments [10,21,40,41]. Guo et al. (2021) demonstrated that a dynamic succession occurred in the rhizosphere microbiome of M. condensatus across the developmental phases on recent Miyake-jima volcanic deposits [8]. Moreover, a unique community structure was determined for the Miscanthus rhizosphere forming on the Miyake-jima volcanic deposits [8]. This community structure was obviously different from the early soil microbial communities in the recent Miyake-jima volcanic deposits, as well as from the rhizosphere communities of the Miscanthus grass colonizing various extreme habitats in Taiwan. The former was dominated by Beta-proteobacteria [6], whereas the latter mainly consisted of the phyla Acidobacteria and Chloroflexi [42]. Considering the extremely low nitrogen content of the volcanic substrate (TN ≤ 0.2 g/kg dry sample in 2016-2018), the present

Discussion
Previous studies have revealed the role of Miscanthus as pioneer plants in disturbed and acidic soil as well as their ability to form root-symbiotic relationships with diazotrophic bacteria and arbuscular mycorrhizal fungi to promote their growth in unfavorable environments [10,21,40,41]. Guo et al. (2021) demonstrated that a dynamic succession occurred in the rhizosphere microbiome of M. condensatus across the developmental phases on recent Miyake-jima volcanic deposits [8]. Moreover, a unique community structure was determined for the Miscanthus rhizosphere forming on the Miyake-jima volcanic deposits [8]. This community structure was obviously different from the early soil microbial communities in the recent Miyake-jima volcanic deposits, as well as from the rhizosphere communities of the Miscanthus grass colonizing various extreme habitats in Taiwan. The former was dominated by Beta-proteobacteria [6], whereas the latter mainly consisted of the phyla Acidobacteria and Chloroflexi [42]. Considering the extremely low nitrogen content of the volcanic substrate (TN ≤ 0.2 g/kg dry sample in 2016-2018), the present study focused on the nitrogen-cycle driven by root-associated microorganisms when this pioneer plant developed on the volcanic substrate.
Total carbon and nitrogen contents in the soil increased three to seven times from 2016 to 2018 (Table 1), appearing to be simultaneous with the rapid development of M. condensatus in the same period. Rapid increases in soil organic matter may promote nitrate reductive processes conducted by heterotrophic microorganisms [43,44], and we found that ammonium content was two-to six-fold higher than that of NO 3 − in all substrates in 2018, which was contrary to the tendency in 2016 (i.e., NO 3 − > NH 4 + , Table 1). These results implied that nitrate reductive processes, i.e., the respiratory denitrification and DNRA, might occupy a major part of nitrogen cycling in the M. condensatus-colonized volcanic substrate, particularly where present in the rhizosphere.
In the metagenomic analysis, we identified genes encoding chemotaxis and Type III/IV/VI secretion systems, RpoN that mediates gene transcription involved in nitrogen metabolism and other cellular processes [45,46], bacterial ammonium transporters, and the spoT/relA homologue gene encoding eubacterial ppGpp. The latter is a mediator of the stringent response that coordinates a variety of cellular activities in response to changes in nutritional abundance (Table 2). These findings suggested that the rhizosphere is the narrow region of soil directly influenced by root activities, preserved as the primary site for plantmicrobe interactions, and adapted to low nutrient conditions [47]. In addition, the plant roots actively secreted specific compounds that altered the bacterial community around the rhizosphere [48]. Thus, a selective association between plants and the rhizosphere microbiome will determine their ability to adapt and survive in harsh conditions [49].
Respiratory denitrification is characterized by the formation of N 2 and/or nitrous oxide (N 2 O) gas as a product of NO 3 − or NO 2 − reduction, followed by a greater increase in strain growth compared with if nitrogen oxide served as a terminal electron acceptor [37]. Our metagenomic results also showed that nirK and nirS were mainly affiliated with the classes Alpha-, Beta-and Gamma-proteobacteria ( Figure 2). Anaerobically isolated strains IG7-14 (i.e., OTU14) and IG7-17 (i.e., OTU17), which belonged to the genus Cupriavidus (Beta-proteobacteria), showed higher in vitro denitrifying activity than other cultured strains ( Figure 4); however, their 16S rRNA gene abundances were low in the rhizosphere (Figure 3). On the contrary, the other strains with low in vitro denitrifying activity seemed to dominate in the rhizosphere, which included members in the classes Flavobacteria, Alpha-proteobacteria, and Gamma-proteobacteria (Figures 3 and 4).
DNRA activity is mainly performed through an anaerobic microbial pathway, but Morley and Baggs (2010) demonstrated that it occurs even when O 2 is as high as 21% (v/v) [50]. This process retains the nitrogen source in the form of NH 4 + [51,52], which serves as N-nutrition for primary producers and heterotrophic microbes [12]. In the present study, a large number of nirB and nirD genes encoding the large and small subunits of cytoplasmic NADH-dependent nitrite reductase, respectively, were detected in the Miscanthus rhizosphere (Figure 2), suggesting that fermentative DNRA bacteria occupied the nitrogen-reducing community associated with the rhizosphere. In addition to the fermentative DNRA-related genes, a considerable proportion of nrfA genes responsible for respiratory DNRA, which is closely related to those possessed by Gamma-proteobacteria, was also found in the rhizosphere (Figure 2).
Several studies have compared the contribution of DNRA to the total NO 3 − consumption in biomes such as temperate and subtropical forests as well as temperate grasslands, which range from negligible to dominant, indicating the regulation of rates and the importance of DNRA by various environmental factors such as redox potential, C/N ratio, quality of carbon, NO 2 − /NO 3 − ratio, sulfur, and iron concentrations [12,53]. Recently, DNRA activity in a soil ecosystem was proposed to be increased or may outcompete denitrification by high carbon availability, and NO 3 − limitation is often expressed as the molar carbon/NO 3 − ratio above 12 [54,55], and such conditions occurred very often in the rhizosphere [56]. Moreover, DNRA rates were affected by the type of plant root exudate (mainly consisting of organic acid) following the order oxalate > citrate > glucose > acetate > malate [57]. In addition, genomic analyses have shown that some Beta-proteobacteria denitrifiers also possess the nirBDC-cysC gene cluster of fermentative DNRA in their genome [58,59], which are thought to function in ammonification and detoxification of NO 2 − when a huge amount of NO 3 − is respirated to NO 2 − by Nar [60]. Thus, the specific adaptive features of plants, as well as the active bacterial role in DNRA and denitrification, might be affected by the nutrient source type and availability around the Miscanthus rhizosphere [61,62].
A meta-analysis identified the classes of Gamma-proteobacteria and Alpha-proteobacteria consisting of the orders of Rhodobacterales, Rhodospirillales, Aeromonadales, and Alteromonadales to be the most important bacterium performing DNRA [12]. Our results confirmed that a large amount of some members belonging to these taxa inhabited the Miscanthus rhizosphere. In respect to the poor soil fertility and the heavy amount of annual precipitation on Miyake-jima, the aggressive association with DNRA bacteria in the rhizosphere may contribute toward the growth of Miscanthus on barren volcanic deposits by retaining nitrogen in the rhizosphere. For this reason, it could be suggested that the grass M. condensatus on the Miyake-jima volcanic deposits showed a specific pattern to shape a beneficial rhizobacterial community during their colonization.
Conversely, the negligible numbers of nitrification-and anammox-related genes in the current study might be explained by two reasons: (1) nitrification in acidic soil might be primarily performed by the fungal community [63] and other unsolved heterotrophic nitrifiers [64]; (2) the contribution of anammox in terrestrial ecosystems is usually low because of their specific ecological requirements [65][66][67][68]. Indeed, a few of the fungal contigs associated with the genera of Ilyonectria, Parastagonospora, Aspergillus, Nectria, Diplodia, Cordyceps, Acremonium, Kwoniella, and Rhodotorula have been detected in the metagenome dataset of the Miscanthus rhizosphere (data not shown).

Conclusions
In conclusion, the present study investigated the nitrogen-nutrient state of recent Miyake-jima volcanic deposits that were fully vegetated by the pioneer plant M. condensatus and evaluated the nitrogen-cycling microorganisms inhabiting the rhizosphere using a combination of shotgun metagenomics and culture-based analyses. Our results showed that nitrogen nutrients in the volcanic deposit increased with the development of the Miscanthusdominated vegetation, and a unique nitrate-transforming bacterial community consisting of denitrifiers and DNRA-bacteria inhabiting the rhizosphere. Additionally, diverse bacteria encompassing three phyla (Proteobacteria, Actinobacteria, and Bacteroidetes) were isolated from the Miscanthus-rhizosphere through stimulation using nitrate and succinate. Those affiliated with the genus Cupriavidus within the class Beta-proteobacteria showed the highest in vitro denitrifying activity; however, the shotgun metagenomic results indicated that the DNRA-related rhizobacteria appeared to take a major role in the nitrogen cycle in the rhizosphere of Miscanthus pioneering on recent Miyake-jima volcanic deposits.
Metagenomics analysis results provide important evidence in the initiation of ecosystem development in Miyake-jima primary volcanic deposits. Moreover, this study increases understanding of the association between the Miscanthus rhizosphere and a nitrate-transforming bacterial community on the newly placed volcanic deposits. Further studies through the metatranscriptomic and/or metaproteomic approaches as well as metabolic stable isotope labelled nitrogen are needed to determine the exact nitrogen-cycle processes that are actively present during Miscanthus colonization.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/microorganisms11020260/s1, Figure S1: Pictures shown the roots of Miscanthus collected from Miyake-jima in 2018, Figure S2: Dendogram based on T-RFLP profiles shows the similarity among the culturable communities in different culture conditions, Figure