Soil bacterial community structures in relation to different oil palm management practices
Site description and soil sampling
The experiment was established as part of the EFForTS project (Ecological and socioeconomic Functions of tropical lowland rainForest Transformation Systems) in the Jambi province, located in Sumatra, Indonesia8.
The experimental sites are located in the state-owned oil palm plantation PTPNVI, which was planted in 2002 (Fig. 1). All planted palms were derived from Tenera seedlings, which are a crossing between Dura and Psifera palms, supplied by Marihat (Medan, Indonesia). Four different locations (referred to as OM1-4) harbor four treatments, which were established in November 2016. In each of these 16 plots (50 × 50 m), five subplots were randomly established, resulting in 80 samples total.
Fertilizer treatment was conducted in two intensities: for one application the conventional treatment usually used in the entire plantation with 130 kg nitrogen, 25 kg phosphorus and 110 kg potassium ha−1 and reduced fertilization with 68 kg nitrogen, 8.5 kg phosphorous and 93.5 kg potassium ha−1. Additionally, liming was conducted in all plots with equal amounts (213 kg dolomite and 71 kg micromag (micronutrients) ha−1). Fertilizer application and liming was done twice per year. The herbicide treatment used 375 cm3 glyphosate ha−1 sprayed within the palm circle four times per year and 375 cm3 glyphosate ha−1 in inter-rows applied twice per year15. The last application before sampling was done in April 2017. Mechanical weeding was done by cutting vegetation four times per year within the palm circle and two times per year in interrows with a brush cutter. The combination of these applications resulted in four different treatments: conventional fertilization with herbicide spraying (ch), conventional fertilization with mechanical weeding (cw), reduced fertilization with herbicide spraying (rh) and reduced fertilization with mechanical weeding (rw) (Table 1).
Topsoil was sampled in May 2017 with a soil corer from the upper seven centimeters in each subplot with a diameter of five cm. A soil corer was used to take three cores in each subplot with a distance of 1 m to each other and at least 1 m distance to trees. The three bulk soil samples per subplot were homogenized and coarse roots and stones were removed. To prevent nucleic acids, especially RNA, from degradation RNAprotect Bacteria Reagent (Qiagen, Hilden, Germany) was applied in a ratio of 1:1. For measurements of soil parameters, we collected an additional sample, which was not supplemented with RNAprotect solution. All samples were transported in cooling boxes and stored at −80 °C until further use.
Nucleic acid extraction
Frozen samples were thawed on ice. RNAprotect was removed from all samples by centrifuging for 20 min at 804.96 g and 4 °C and discarding the resulting supernatant. DNA and RNA were co-extracted from 1 g of soil by using the Qiagen RNeasy PowerSoil Total RNA kit and the RNeasy PowerSoil DNA Elution kit as recommended by the manufacturer (Qiagen), except that RNA was eluted with 50 µl elution buffer instead of 100 µl. DNA contamination was removed from RNA preparations by using the TurboDNAfree kit (Applied Biosystems, Darmstadt, Germany). For this purpose, 0.1 volume DNAse buffer and 1 µl DNAse were added and incubated for 30 min at 37 °C. Subsequently, a second digestion cycle was performed with 0.5 µl DNAse at 37 °C for 15 min. RNA was then purified with the RNeasy MiniElute Cleanup kit (Qiagen). In order to verify complete DNA removal, a control amplification of the 16 S rRNA gene was performed as described below for 16 S rRNA gene amplification. Purified RNA was then reverse-transcribed into cDNA with the Superscript IV reverse transcriptase and a specific primer (5′-CCGTCAATTCMTTTGAGT-′3) as recommended by the manufacturer (Thermo Fisher Scientific, Schwerte, Germany). After cDNA synthesis, we removed residual RNA by adding 1 µl RNase H (New England Biolabs, Frankfurt am Main, Germany) to each reaction and incubation for 20 min at 37 °C. Obtained DNA and cDNA were stored at −20 °C until further use.
16 S rRNA gene amplification and sequencing
For amplification of 16 S rRNA sequences, we used 16 S rRNA gene primers targeting the V3-V4 region (forward primer: S-D-Bact-0341-b-S-17 5′-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG-CCTACGGGNGGCWGCAG-3′, reverse primer: S-D-Bact-0785-a-A-21 5′-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG-GACTACHVGGGTATCTAATCC-3′) as described by Klindworth22 and Herlemann23 and added adapters for MiSeq sequencing (underlined). PCR reactions were performed in a total volume 50 µl containing 10 µl of 5-fold Phusion GC buffer, 0.2 µl 50 mM MgCl2 solution, 2.5 µl DMSO, 200 µM of each of the four deoxynucleoside triphosphates and 1 U of Phusion High-Fidelity DNA Polymerase (Thermo Fisher Scientific). We used 20 to 30 ng of DNA and 1 µl cDNA per reaction. The PCR reaction was started by an initial denaturation at 98 °C for 1 min, followed by 25 cycles of denaturation at 98 °C for 45 s, annealing at 60 °C for 45 s and elongation at 72 °C for 30 s. The final elongation was at 72 °C for 5 minutes. Amplicons were then purified by using MagSi-NGS PREP Plus magnetic beads following the procedure recommended by the manufacturer (Steinbrenner Laborsysteme GmbH, Wiesenbach, Germany) with the Janus Automated Workstation from Perkin Elmer (Perkin Elmer, Waltham Massachusetts, USA). Illumina MiSeq sequencing adapters were attached to the purified amplicons with the Nextera XT Index kit (Illumina, San Diego, USA). The Index PCR was done by using 5 µl of template PCR product, 2.5 µl of each index primer, 12.5 µl of 2x KAPA HiFi HotStart ReadyMix and 2.5 µl PCR grade water. Thermal cycling scheme was as follows: 95 °C for 3 min, 8 cycles of 30 s at 95 °C, 30 s at 55 °C and 30 s at 72 °C and a final extension at 72 °C for 5 min. The indexed products were purified as described before. Products were quantified by using the Quant-iT dsDNA HS assay kit and a Qubit fluorometer following the instructions of the manufacturer (Invitrogen GmbH, Karlsruhe, Germany). Purified amplicons were sequenced by the Göttingen Genomics Laboratory with a MiSeq instrument with a read length of 2 × 300 bp using dual indexing and reagent kit v3 (600 cycles) as recommended by the manufacturer (Illumina).
Sequence processing
We obtained 6,817,019 amplicon sequences with 5,183,993 remaining sequences after quality-filtering from DNA samples. At RNA level 6,412,838 raw sequences with 3,601,637 remaining sequences after quality-filtering were obtained24.
Obtained paired-end sequences were first quality-filtered with fastp version 0.2025 using a minimum phred score of 20, a minimum length of 50 bases, the default sliding window size (–cut_window_size = 4), read correction by overlap (option “correction”), adapter removal of the sequencing primers (option “adapter_fasta”), and the provided index sequences of Illumina. Quality-filtered paired-end reads were merged with PEAR version 0.9.11 and default settings26. Primer sequences were clipped with cutadapt version 2.5 and default settings27. All further steps, except mapping of sequences to ASVs (Amplicon Sequence Variant) were performed with functions implemented in vsearch version 2.1.4.128. Sequences were filtered by size with “sortbylength” with a set minimum length of 300 bp. Dereplication of identical sequences was done by “derep_fulllength”. Denoising and removal of low abundant sequences with less than eight replicates were done with the vsearch UNOISE3 module “cluster_unoise”. Chimeric sequences were removed by employing the UCHIME module of vsearch. This included a de novo chimera removal (“uchime3_denovo”) and a reference-based chimera removal (“uchime_ref”) against the SILVA SSU 138 NR database29. Sequences were mapped to ASVs by vsearch (“usearch_global”) with a set sequence identity threshold of 0.97. Taxonomy assignments were performed with BLASTN30 (version 2.9.0) against the SILVA SSU 138 NR database29 with an minimum identity threshold of 90%31. In addition to the taxonomy identity, we added the taxonomy id of the database, length of fragment, query percentage identity, query coverage and e-value in the taxonomy string of the table. We used identity (pident) and query coverage (qcovs) per ASV of the blast output to exclude uncertain blast hits. As recommended by the SILVA ribosomal RNA database project32, we removed the taxonomic assignment for blast hits if dividing the sum of percent identity and percent query coverage by 2 resulted in ≤93%. In total, 31,987 ASVs were used for downstream analysis.
Bacterial community analysis
The bacterial community composition was further analysed in R33 (version 3.6.1) and RStudio34 (version 1.1.463). ASV counts were normalized by using the Geometric Mean of Pairwise Ratios (GMPR) of the GMPR package version 0.1.335. Community compositions were then analysed by the ampvis2 package version 2.4.11 and “amp_heatmap” at genus level36. The fifteen most abundant genera were displayed as relative abundance and clustered at treatment level. Heat-trees were displayed by the metacoder37 package (version 0.3.2.9001).
For heat-tree calculation all counts were summed at order level and all taxa with a relative abundance of More
