Community structure of total bacteria and methane emission-related prokaryotes in the rice fields applied with urea and biofertilizer

Rice fields are a source of methane emissions. Urea fertilization is considered to increase methane emission in rice field. Reduction in amount of urea applied with addition of biofertilizer, consisting of methanotrophic and N2O-reducing bacteria, is presumably to become an innovative fertilization technique to decrease methane emission from rice field. This current work aimed to investigate the community structure of total bacteria and methane emission-related prokaryotes in rice field soil treated with urea and biofertilizer at the vegetative and generative of rice stage. Two treatments were set up in the field experiment, i.e., 100% urea (250 kg/ha) without biofertilizer (B0) and 50% urea (125 kg/ha) with biofertilizer (B1). We used Illuminabased sequencing to investigate the soil microbial community in each treatment. Results showed that the soil bacterial community had minor changes in the two treatments throughout the rice growing period. Application of 50% urea with biofertilizer (B1) did not change the dominant bacterial phyla in rice field soil, i.e., Proteobacteria. However, there were differences in bacterial composition among the two treatments. Bacterial communities were partitioned into two clusters by the treatments (B0 and B1) rather than the rice growth phase. In addition, methanogens:methanotrophs ratio in the B1 treatment was lower than that of the B0 treatment.


Introduction
Global warming, which is currently the focus of world attention, is triggered by an increase in the concentration of greenhouse gases (GHG) in the atmosphere. Carbon dioxide (CO2), methane (CH4), and nitrous oxide (N2O) as the primary GHG have increased by 40%, 150%, and 20%, respectively, since 1750 to 2010 (IPCC, 2014). In 2010, GHG concentration in the atmosphere comprises 76% CO2, 16% CH4, 6.2% N2O, and 2.0% fluorinated gases (Fgases) (IPCC, 2014). Atmospheric CH4 concentration in 1750 was 700 ppb, then it increased to 1,762 -1,893 ppb in 2012 (USEPA, 2015). CH4 remains in the atmosphere for approximately 8.4 years. Moreover, CH4 possess a global warming potential (GWP) of 25 times stronger compared to CO2 for a 100-year timescale, in which small changes of CH4 concentration in the atmosphere significantly affect global warming (IPCC, 2007;Bridgham et al., 2013). CH4 emission is resulted from natural source and anthropogenic source with the percentage range of 21-47% and 45-80% to the total global emissions of 600 Tg of CH4 per year, respectively (Liu and Whitman, 2008). Rice (Oryza sativa L.) cultivation has been widely known as the main contributor to CH4 emission, which is the third-largest CH4 anthropogenic source subsequent to ruminants and energy generation (IPCC, 2007;Liu and Whitman, 2008). CH4 emission in rice cultivation ranges from 30 to 50 Tg per year or 10-20% of the total CH4 global (Philippot et al., 2009). Rice cultivation in Indonesia is a pivotal anthropogenic activity that produces GHG and reported as one of the top five contributing countries to GHG emission (USEPA, 2013). Indonesia is predicted to emit CH4 by 4 Tg per year or 6-7% of the total annual CH4 global from the rice cultivation sector (Husin et al., 1995). Rice cultivation in Indonesia is mainly carried out using conventional irrigation system. Memberamo is one of the most cultivated-rice varieties in Indonesia. Setyanto et al. (2004) reported that the average of 61.1 kg CH4/ha was produced by Memberamo cultivar under an irrigated field condition. CH4 flux of Memberamo cultivar reached ~100 mg/m 2 /d both at 35 days after transplanting (DAT) (the rice vegetative stage) and 67 DAT (the rice generative stage). Meanwhile, it achieved about 85 mg/m2/d at 83 DAT (the harvesting day). In relation to global warming, rice growth is facing to high temperature stress (Fahad et al., 2016a;Fahad et al., 2016b;Fahad et al., 2016c). High temperature stress decreased the pollen fertility, anther dehiscence, pollen retention, germination and metabolites synthesis in pollens of rice plants (Fahad et al., 2016a). High daytime and nighttime temperature significantly decreased the photosynthetic reaction of rice plant (Fahad et al., 2016b). Under high temperature at nighttime, a decrease in individual grain weight significantly decreases the production of rice grain per unit area (Fahad et al., 2016b). Temperature stress had been reported to affect the growth and yield performance of rice, which could reduce rice yield and the sustainability of its production (Fahad et al., 2016b;Fahad et al., 2016c). In addition, Ceccarelli et al. (2010) predicted that the temperature increase by up to 2 o C will decrease rice yield up to 41% in the end of 21 st century. Therefore, more efforts is necessary to meet rice demand. Increasing in rice cultivation area worldwide may be conducted to supply rice needs along with an increase in the human population. Consequently, rice demand improvement might be actively involved in CH4 emission. Atmospheric CH4 is mainly produced by biological activities, ranging from 70-80%. Complex microbial community inhabiting the rice field soil contributes to CH4 emissions. Methanogens, members of archaea, play significant roles in CH4 production under anaerobic conditions in rice field. Then, CH4 is used as a substrate for methanotroph, belonging to archaea and bacteria that are able to oxidize CH4 (Hanson and Hanson, 1996). CH4 fluxes are governed by a net balance between methanogen and methanotroph. Many factors influence microbial community in the rice field, such as crop rotation (Zhao et al., 2014), soil characteristics (Curd et al., 2018), water management (Breidenbach and Conrad, 2015), rice genotype (Shenton et al., 2016), rice age and growth phase (Edwards et al., 2018), fertilizer (Seo et al., 2014), and microbial inoculation (Bao et al., 2013). Microbes respond to environmental changes, including anthropogenic activities, then affect biogeochemical processes (Bissett et al., 2013). Manipulation of microbial communities in the soil might be a mitigation strategy in anthropogenic GHG emission (Singh et al., 2010). In previous studies, methanotrophic and N2O-reducing bacteria were used as biofertilizer to reduce CH4 and N2O emissions in the rice field. Application of the biofertilizer with the reduced amounts of synthetic fertilizer (NPK or urea) not only enhanced rice plant growth but also reduced CH4 and N2O emissions (Pingak et al., 2014;Sukmawati et al., 2015;Fatma, 2019a). However, the characteristic of the bacterial community in the treated rice filed soil has not been well documented. The information about the methanogens:methanotrophs ratio in relation to CH4 emission in the treated rice field is also still unclear. This current work aimed to investigate the community structure of total bacteria and methane emission-related prokaryotes in the Indonesian rice 3/12 Asian J Agric & Biol. 2021(2). field soil applied with urea and biofertilizer in relation to CH4 emissions at the rice vegetative and generative stage.
Understanding microbial community structure and their shift in the rice field is necessary in order to develop CH4 mitigation strategy. In this study, we analyze the microbial diversity and composition using an Illumina-based sequencing technique. This study provides basic information regarding the bacterial communities in the rice field under urea and biofertilizer application for Indonesian rice agroecosystems.

Study site description
This study was performed in the rice field in Tegal district, Central Java, Indonesia (7°02'44.8 "S, 109°08'16.6 "E). The rice soil had neutral pH (6.55) and classified as silt loam soil containing 23.31% sand, 75.98% silt, and 0.71% clay. The soil contained 3.30% of C (high category), 0.12% of N (low category), and 40.70 of C:N ratio (very high category).

Field experiment
Rice plants were grown in the rice field according to the existing practices of local farmers. The field experiment was designated to have two treatments, i.e., an initial amount of 100% urea (250 kg/ha) without biofertilizer (B0) and 50% urea (125 kg/ha) with biofertilizer (B1). Biofertilizer used in this study was a mixture of methanotrophic bacteria (Methylocystis rosea BGM1, Methylocystis parvus BGM3, Methylococcus capsulatus BGM9, and Methylobacter sp. SKM14) and N2Oreducing bacteria (Ochrobactrum anthropi BL2) with a density of 10 8 cells/mL (Rusmana and Akhdiya, 2010;Setyaningsih et al., 2010). The treatments were performed at three blocks, where each block divided into two plots designated for B0 and B1 with the total area were 1,500 m 2 , respectively. Rice seedling (Indonesian rice cultivar 'Memberamo') were transplanted correspond to each treatment, and marked as 0 day after transplanting (DAT). Biofertilizer in the B1 treatment was applied through root inoculation before seedlings transplantation. Urea application on both treatments was given at 15 DAT. Rice plants were grown under continuous flooding until the rice generative phase (approximately 15 days before harvesting), and dried until the harvesting day. Rice growth stage determination was carried out by observing the tiller number and plant height.

Soil samples collection and genomic DNA extraction
Soil samples were taken at three sampling times representing the rice growth phase, including 0 (before treatments), 36 (rice vegetative phase), and 69 DAT (rice generative phase). To investigate whether the bacterial application extends to rice field soil, we randomly collected three replicates of soil microcosms (0-15 cm in depth). Soil microcosms were taken from the site adjacent to the rice roots in each of two treatments (three replicate plots in each treatment) using soil core sampler (diameter 4.5 cm). The same amounts of soil samples from the three replicates of soil samples with the same treatment were pooled to obtain a homogeneous soil sample. A total of 6 pooled soil samples (2 treatments x 3 sampling times) were prepared for DNA extraction, i.e., 0B0, 0B1, 36B0, 36B1, 69B0, and 69B1, in which the number in front of the sample code denotes sampling time, B0 and B1 represent the treatments. Genomic DNA was separately extracted from each soil sample using ZymoBIOMICS DNA Mini Kit (Zymo Research, Irvine, CA, US) following the manufacturer's recommendation with few modifications. The DNA concentration and purity were assessed using NanoDrop (NanoDrop Technologies, Wilmington, US) and electrophoresis on 2% agarose gel. DNA was stored at -20°C for subsequent analysis.

16S rRNA gene amplicon generation and sequencing
Six soil DNA samples were prepared for amplicon sequencing. The V4-V5 region of bacterial 16S rRNA gene was amplified using a two-step PCR protocol recommended by the Illumina to obtain the barcoded amplicon library. In the first PCR step, the V4-V5 amplicon was generated using primer set 515F (5'-GTG CCA GCA GCC GCG GTA A-3') and 907R (5'-CCG TCA ATT CCT TTG AGT TT-3') (Caporaso et al., 2011;Armitage et al., 2012) added with the Illumina overhang adapter sequences. The initial PCR reaction mixture with a total volume of 25 µL consisted of 12.5 µL of Q5 High-Fidelity 2x Master Mix (New England Biolabs, Ipswich, MA, US), 1.25 µL of 10 µM forward primer, 1.25 µL of 10 µM reverse primer, 5 µL of template DNA, and 5 µL of nuclease-free water. The cycling parameter used in the first PCR reaction were 1 cycle at 98°C for 30 sec, followed with 30 cycles of 98°C for 10 sec, 72°C for 30 sec, 72°C for 30 sec, and 1 cycle at Asian J Agric & Biol. 2021(2).
72°C for 2 min. The PCR products from the same sample were pooled in equidensity ratios, then purified using Qiagen Gel Extraction Kit (Qiagen, Hilden, DE). The second PCR step was conducted to add the dual indices, namely Nextera XT i7 and Nextera XT i5 (Illumina Inc., San Diego, CA, US). The PCR reaction mixture was prepared in a total volume of 50 µL, containing 25 µL of Q5 High-Fidelity 2x Master Mix, 5 µL of Nextera XT i7 index, 5 µL of Nextera XT i5 index, 5 µL of PCR product generated in the first PCR step, and 10 µL of nuclease-free water. PCR reaction in the second step was as follows, a pre denaturation at 98°C for 3 min, 8 cycles of 98°C for 30 sec, 55°C for 30 sec, 72°C for 30 sec, and a final extension at 72°C for 5 min. DNA bands between 400-450 bp were chosen for further steps. The PCR products were determined for the quality and quantity using Agilent 4200 TapeStation System (Agilent Technologies, Inc., Santa Clara, CA, US), Helixyte TM Green dsDNA Quantifying reagent (AAT Bioquest Inc, CA, US), and NanoDrop (Thermo Scientific, Wilmington, DE, US). Sequencing of the PCR products obtained from the second PCR step was performed using an Illumina MiSeq platform by the Axil Scientific, Singapore. From the sequencing process, it is obtained 250 bp paired-end reads.

Analysis of sequence data
Obtained paired-end reads were set to the sample according to their specific barcode. The barcode and primer sequence was truncated from the paired-end reads, then the paired-end reads were merged and overlapped using FLASH V1.2.7 (Magoč and Salzberg, 2011). This software is able to merge some reads with the read from the opposite end of the same DNA fragment, generating the splicing sequences called the raw tags. To obtain the high-quality clean reads, the raw tags were filtered for the quality under the specific filtering condition according to QIIME V1.7.0 (Bokulich et al., 2013;Caporaso et al., 2010). For chimera detection, the reads obtained from the previous step were compared to the GOLD database using UCHIME algorithm (Edgar et al., 2011). The chimera sequences were removed to obtain the effective reads (Haas et al., 2011). The effective reads were analyzed using UPARSE V.7.0.1001 (Edgar, 2013) to cluster these sequences into the operational taxonomic unit (OTUs) with ≥97% similarity. Sequences frequently existed within each OTU was screened and chosen to obtain the representative sequence. The representative sequences were taxonomically annotated according to RDP classifier algorithm V2.2 (Wang et al., 2007) via the Greengene database at a confidence threshold of 0.8 (DeSantis et al., 2006). The information about OTUs abundance was normalized using a standard of sequence number, which corresponds to the sample with the lowest sequence number. Normalized data were used to examine alpha diversity and beta diversity. Alpha diversity showed bacterial diversity in a sample. The diversity was displayed as observed species, Chao1, ACE, Shannon index, and equitability index. Alpha diversity was determined using QIIME V1.7.0, then exhibited using R software V2.15.3. Chao1 and ACE were calculated to estimate the species richness in an assemblage. Shannon diversity index were used to identify the richness and evenness of microbial community. Rarefaction curve was used to assess the sequencing depth, showing that the data accurately reflect the microbial community within the sample. Each sample required a different amount of sample for a plateau. Meanwhile, beta diversity compares the similarity and dissimilarity of the bacterial communities between different samples. Beta diversity was assessed through a hierarchical clustering analysis using the Unweighted Pair-Group Method with Arithmetic Mean (UPGMA). This analysis was performed by QIIME V1.7.0, according to the Weighted UniFrac distance matrix.

Correlation analysis between bacterial taxa and plant age
The correlation between the relative abundance of bacteria and the plant age was analyzed using Pearson's correlation analysis. Correlation analysis was conducted for high abundance taxa in phylum and class level of taxonomy (> 1% relative abundance).

Analysis of methanogens:methanotrophs ratio
In this current work, we measured the relative abundance of OTUs classified as methanotrophic bacteria according to bacterial OTUs data of each sample. We separately collected the relative abundance of OTUs assigned to methanogens and anaerobic methanotrophic archaea from our related study which investigates the archaeal community structure from the same samples (Fatma, 2019b). The relative abundance of methanotrophic bacteria and methanotrophic archaea were calculated as

5/12
Asian J Agric & Biol. 2021(2). methanotroph group. Based on their relative abundance, we measured the methanogens:methanotrophs ratio in each sample. Pearson correlation analysis was carried out to determine the correlation between methanogens:methanotrophs ratio and CH4 fluxes that previously observed in our related study (Fatma, 2019a). Fatma (2019a) measured CH4 emission using a closed chamber method at the vegetative and generative of rice stage. Gas samples were collected from the closed chamber with a size of 55.5 x 55.5 x 89.9 cm 3 installed in rice field of each treatment. Gas sampling was conducted using 10 mL-syringe at 0 hours (t=0) and 3 hours after chamber installation (t=3). CH4 concentration was calculated using gas chromatography with flame ionization detector (GC-FID). Subsequently, CH4 emission was calculated according to International Atomic Energy Agency (IAEA, 1992).

Results
Bacterial diversity within the sample A total of 3,074,783 raw reads were obtained from 6 samples within bacterial sequencing reads. After lowquality reads, chimera, and singleton removal, we obtained 384,291 high-quality reads. Table 1 shows the sequences number and alpha diversity of bacterial sequences in all samples. The alpha diversity in B0 (0B0) was relatively similar to B1 (0B1), then steadily decreased from the rice vegetative to generative phase. During the rice growing period, bacterial diversity in B1 was slightly greater than that of B0. We observed the same trend in both treatments, where bacterial diversity decreased from 0 to 69 DAT. The rarefaction curves showed that the plateau was approached in each sample, which indicates that the bacterial sequencing reads were relatively sufficient in performing bacterial community in each sample (Figure 1).

Figure-1. The rarefaction curves of bacterial sequences from rice field in both treatments
during rice growing period. 0B0: 100% urea without biofertilizer at 0 DAT; 0B1: 50% urea with biofertilizer at 0 DAT; 36B0: 100% urea without biofertilizer at 36 DAT; 36B1: 50% urea with biofertilizer at 36 DAT; 69B0: 100% urea without biofertilizer at 69 DAT; and 69B1: 50% urea with biofertilizer at 69 DAT. Bacterial OTUs in the two treatments are presented in the Venn diagram ( Figure 2). We observed the bacterial OTUs shift in each treatment throughout the rice-growing period. Bacterial OTUs frequently present in all observation time was slightly higher in B1 treatment (37.1%) than those of B0 treatment (35.8%). Meanwhile, specific OTUs at each observation time (0, 36, or 69 DAT) in B0 treatment slightly differed from B1.

Bacterial composition in the rice field soil
Here, the result of the taxonomic assignment is shown in phylum and class level (Figure 3). We found minor changes in some bacterial phyla and class during rice growth, even though most of them are not significantly correlated (Table 2). In B0 treatment, plant age was positively correlated with Bacteroidetes and Verrucomicrobia (r > 0.8), and negatively correlated with Acidobacteria, Planctomycetes, Cloroflexi, and Actinobacteria (r < -0.8). The correlation analysis in B1 treatment indicated a positive correlation among plant age and the relative abundance of Proteobacteria and Nitrospirae (r > 0.8). Negative relationship in B1 treatment was observed between plant age and Acidobacteria, Bacteroidetes, and Verrucomicrobia were observed (r < -0.8) ( Table 2). a b

Figure-3. Relative abundance of the bacterial community showing the dynamic of top ten dominant bacterial community at phylum (a) and class level (b) in treated rice field throughout the rice growth period.
The 16S rRNA gene sequences were assigned to the phylum (a) and class (b) level of taxonomy based on the Greengenes database. 0B0: 100% urea without biofertilizer at 0 DAT; 0B1: 50% urea with biofertilizer at 0 DAT; 36B0: 100% urea without biofertilizer at 36 DAT; 36B1: 50% urea with biofertilizer at 36 DAT; 69B0: 100% urea without biofertilizer at 69 DAT; and 69B1: 50% urea with biofertilizer at 69 DAT.

Methanogens:methanotrophs ratio
The relative abundance of putative methanogen and methanotroph were obtained from our archaeal 16S rRNA gene data (Fatma, 2019b). For putative methanotroph, we collected the relative abundance of methanotrophic bacterial and archaeal assemblages. Within bacterial assemblages, we only detected Methylocystaceae in the same proportion in all samples accounting for 0.2%. We did not find methanotrophic bacteria belonging to Verrucomicrobia and NC10 phyla. From the archaeal community, we observed ANME-2D (anaerobic oxidation of methane) with varied proportion for each sample (Fatma, 2019b). We added the information about CH4 flux at the rice vegetative and generative phase (Fatma, 2019a). Table 3 shows the methanogens:methanotrophs ratio and CH4 fluxes in the two treatments. Data showed that methanogens:methanotrophs ratio was positively correlated with CH4 fluxes (r = 0.843, P = 0.157). Higher methanogens:methanotrophs ratio indicated higher CH4 fluxes, as shown in B0 treatment at both vegetative and generative phase of rice plants. Methanogens:methanotrophs ratio at the rice vegetative phase was higher than that of the rice generative phase.

Discussion
Illumina-based sequencing showed that obtained sequences had represented the whole bacterial diversity in each sample (Figure 1). Before fertilizer application, bacterial diversity was relatively similar between both treatments, then the bacterial diversity decreased over the rice growth phase (Table 1). The decrease in bacterial diversity from the rice vegetative phase to the generative phase was also observed in previous investigations (Pittol et al., 2018;Edwards et al., 2018). It was possibly due to root exudates and nitrogen availability were declined. We found that bacterial diversity was slightly higher in B1 in comparison with B0, indicating that B1 treatment led to minor effect the bacterial diversity in the rice field soil (Table 1). Additionally, bacterial OTUs commonly found at 0, 36, and 69 DAT was slightly higher in B1 compared to B0 treatment ( Figure 2). Inoculation of methanotrophic bacteria in B1 treatment might support the growth of certain bacterial OTU from 0 to 69 DAT. Bacterial inoculation may affect soil resident microbial growth (Bao et al., 2013). Soil from B0 treatment has less bacterial diversity, which was likely due to a higher urea application rate. Our result was consistent with a previous report which showed that nitrogen fertilizer declined soil bacterial diversity (Yang et al., 2019). Phylogenetic analysis revealed that Proteobacteria was the most common phylum in rice field soil of the two treatments. We suggest that Proteobacteria dominance is likely due to high diverse phylum members with a wide variety of metabolic activity, particularly chemoorganotrophs and chemolithotrophs. They may highly adapt in rice field soil with a variety of physicochemical characteristics caused by anthropogenic activity (Pittol et al., 2018). Proteobacteria, Acidobacteria, Planctomycetes, Chloroflexi, Actinobacteria, and Bacteroidetes made up approximately 80% of the total bacterial sequences within each sample ( Figure  3a), which is in accordance with prior observations (Lee et al., 2014;Pittol et al., 2018). Our findings suggest that the application of 50% urea and biofertilizer (B1) affected minor changes in the soil bacterial community at phylum and class level of taxonomy ( Figure 3, Table 2). The bacterial community at the phyla and class level responded to the different fertilizer treatment throughout the rice growing phase, even though most of them were not significant (Table 2, Figure 3). Consistent with our study, little effects on rice-associated bacteria have been reported in rice-seedling inoculated by Azospirillum sp. B510 (Bao et al., 2013). On the contrary, no significant impact on rhizosphere microbial communities was observed in rice inoculation by a commercial formulation consisting of Pseudomonas fluorescens and Azospirillum brasilense (Salamone et al., 2012). Different fertilizer treatments may influence bacterial community and their interaction with the environment. Zeng et al. (2016) reported that nitrogenous fertilizer directly decreased bacterial diversity by ammonium availability enhancement and indirectly decreased bacterial community composition by soil pH acidification. Declining soil pH may unfavorable to soil microbes which were due to most of them are more adaptive in neutral conditions. Soil microbes had the highest phylogenetic diversity in neutral pH . In this present study, the decrease in bacterial diversity caused by higher nitrogen application rate was shown in B0 treatment (Table 1). It might be caused by nitrophilous species enrichment, which outcompetes other species (Bobbink et al., 2010). The changes of some bacterial phylum were likely due to other soil characteristic influenced by nitrogen application (Figure 3, Table 2). Bacterial communities were partitioned into two clusters by fertilizer treatment rather than the rice growth phase (Figure 4). Nitrogen input seemed to cause the rice soil more oxidative and slightly changed the abundance of certain groups, including Proteobacteria . Additionally, urea application declined redox potential (Eh) around the rice root region and dissolved oxygen (DO) at the soil-water interface (Liu et al., 2019). The association between bacterial community, soil, and plant responds to some changes in soil characteristics by 9/12 Asian J Agric & Biol. 2021(2). fertilizer addition, including soil pH, labile N availability (Zeng et al., 2016;Yang et al., 2019), Eh, DO (Liu et al., 2019), soil organic carbon and cation exchange capacity (CEC) . Members of the same bacterial taxa in phylum and class level did not always respond to nitrogen addition in a similar way. Their response is controlled by taxa abundance at the lower level of taxonomy, yet some minor changes in the lower level were not always exhibited at the higher level (Zeng et al., 2016). The application of nitrogen fertilizer also affects microbial community related to CH4 emission. Urea addition sharply elevated methanogenic 16S rRNA gene and methanotrophic pmoA gene copies, increasing CH4 emission from rice field ecosystems (Liu et al., 2019). The abundance of methanogens was enriched in the soil with standard nitrogen input rather than lower nitrogen input. In contrast, low nitrogen input was significantly increased the population size of methanotrophic bacteria and CH4 oxidation activity on rice roots . Our experiment showed that the application of 100% urea (B0) causes the relative abundance of methanogens were higher in comparison to the relative abundance of methanotrophs (Table 3). Inorganic nitrogen could trigger the relative abundance and activity of methanotrophic bacteria, as reported by a previous investigation (Seo et al., 2014). In this study, lower nitrogen input in B1 treatment might stimulate methanotrophic activity, leading decreases in CH4 emission rate. On the contrary, inorganic nitrogen may inhibit methaneoxidizing bacteria (Mishra et al., 2018). However, the effect of nitrogen fertilizer application on methanogenic and methanotrophic community and CH4 emission remains in debate. In the recent study, methanotrophic bacteria belonging to Methylocystaceae was detected in low relative abundance of 0.2% in each sample (Fatma, 2019b) This number is lower than previous studies, which reports that 0.79-1.75% and approximately 1.1% of bacterial sequencing reads were classified as aerobic methanotroph (Lee et al., 2014;Vaksmaa et al., 2017). We suggest that Methylocystaceae commonly presents in all soil samples and plays essential roles in CH4 emission. Methylocystaceae, such as Methylosinus (type II methanotrophic bacteria), are major contributors to CH4 oxidation and N2 fixation in rice roots . Methanogens and methanotrophs communities were possibly correlated with CH4 emission. Lee et al. (2014) reported that pmoA:mcrA transcript:gene ratio plays prominent roles in CH4 production in the rice field, and could be the best predictor of actual CH4 emission rate. In our study, methanogen:methanotroph ratio was positively correlated with CH4 flux, even though it was not significant (Table  3). Higher methanogen:methanotroph ratio showed higher CH4 fluxes, as observed in both treatments at the rice vegetative phase compared to the generative phase (Table 3). A previous study reports that CH4 emission was higher at tillering phase when soil was flooded (Sheng et al., 2016). Methanogen was found in larger population size and activity at the vegetative phase because of anaerobic condition caused by flooding (Breidenbach and Conrad, 2015). Meanwhile, CH4 oxidation activity was limited by O2 and CH4 concentration, yet CH4 was actively oxidized in the oxic soil surface layer (Conrad and Rothfuss, 1991). Interestingly, inoculation of methanotrophic bacteria in B1 treatment led to lower methanogens:methanotrophs ratio and CH4 fluxes along the rice growth period (Table 3). This result showed that methanotrophic bacteria is highly possible to decrease CH4 emission in the rice field by reducing the methanogens:methanotrophs ratio. To improve our knowledge of the microbial community and function in relation to CH4 emission, further studies with more replicates would be needed. Additionally, metatranscriptomic analysis may improve our understanding of the changes in microbial community to evaluate the effect of the biofertilizer in the rice field soil related to CH4 mitigation.

Conclusion
The bacterial communities showed a minor dynamic during the rice-growing period under different fertilizer treatments. Proteobacteria was the most dominant bacterial phyla in both treatments along rice growing period. The application of 50% urea with biofertilizer seemed to have a minor effect on the bacterial diversity and composition in the rice field soil. Although the dominance of bacteria phyla did not differ between treatments, the bacterial communities were partitioned into two clusters by the treatments (B0 and B1). Within microbial assemblages, methanogens:methanotrophs ratio in B1 treatment was lower compared to B0 treatment. The methanogens:methanotrophs ratio was positively correlated with CH4 fluxes in the rice field. The biofertilizer combined with reduced amounts of urea would be an alternative candidate for CH4 mitigation.