The study was done at two locations in Trinidad (10° 3′N 60° 55′W and 10° 50′N 61° 55′W) on forests owned by the Ministry of Agriculture (Figure S1). Investigations were conducted on Teak plantations (T) and adjacent native forests (control) in the Tamana site (hereafter referred to as the north location) and the Brickfield, Tabaquite site (hereafter referred to as south location). Briefly, Teak is a deciduous tree and belongs to the Lamiaceae family. The native north and south stands were both dominated by Spondias monbin (Hog Plum), Brownea latifolia (Cooper Hoop), Pentaclethra macroloba (Bois multare), Carapa guianensis (Crappo) and Trichilia oblanceolata (Acurel). Additional tree species in the North were Tabebuia serratifolia (Yellow pouis), Swartzia pinnata (Bois pois), Eschweilera subglandulosa (Guatecare), Manilkara balata (Balata) and Sterculia caribaea (Mahoe). In the south, additional species were: Pachira insignis (Wild chataigne), Trichilia trinitensis (Obi), Pisonia eggersiana (Jiggerwood) and Isertia parviflora (Bois Fer).

Teak stands of three age groups were analyzed. One group was the oldest Teak trees that were established in 1943 and located only at the south site. A second set was middle-aged plantations that were established between 1956 and 1968 at both the north and south locations. The third set was a stand of young trees that was established in 2009 at the north site only. Hereafter, the Teak stands are abbreviated as from the North (N) or South (S) locations preceded by the year of establishment: 1957 (57N), 1966 (66N), 2009, (09N), 1943 (43S), 1956 (56S), 1968 (68S). Native forests from the north and south locations are designated as NatN and NatS, respectively.

Within each selected stand, soil samples and associated data were taken on a randomly determined line segment. For each stand, the length of the line transect was 200 m and five sampling points were determined systematically. Five replicate composited samples from each stand were obtained at a distance 40 m apart. At each of the five sampling point, a soil auger was used to collect five soil samples from the top 0 to 20-cm (A horizon) thus giving 25 soil samples. Hence, each of the five sample points consisted of a composite of five soil samples which were kept separate for all further analyses. Samples were transported to the laboratory where they were homogenized by 2-mm sieving, and fine roots were removed. Samples were then partitioned for analysis of enzymatic activity, PLFA extraction, DNA extraction, and edaphic characterization. Enzyme analysis and pH measure was done on all five of the samples. The remainder of the analyses was done on three samples randomly chosen from the set of five.

Soil pH was determined by the slurry method (1:5, w:v; sample:distilled deionized water) measured with an Eijkelkamp pH/mV/EC/Salinity/T/02m (Agrisearch Equipment ZG Giesbeek, the Netherlands). Total carbon (TC), organic carbon (OC), and inorganic carbon were determined by dry combustion using a LECO CNS-2000 analyzer. Total N (TN) was determined by combustion and by the Kjeldahl digestion method (TKN). Nitrate-N was measured in soil water extracts colorimetrically after reaction with phenoldisulphonic acid, and ammonium-N was determined by flow injection analysis. Major and minor trace elements were determined by inductively-coupled mass spectrometry. Cation exchange capacity (CEC) was determined by ammonium acetate extraction proceeded by cation quantification via atomic absorption spectrometry. Soil physicochemical characteristics are displayed in Table 1 and Table S1.

EAA was conducted according to Saiya-Corka et al. (2002). Soil (1 g wet weight) was added to 125 mL sodium acetate buffer (50 mM, pH 5) and homogenized in a Waring 7009G blender (Waring Laboratory Science, Torrington CT) run at top speed (22,000 RPM) for 1 min. The soil slurry was then decanted into a flat bottom tray containing a stir bar for rapid mixing to maintain the soil suspension, and 12 aliquots (200 μl each) were transferred simultaneously to a 96-well microplate (Thermo Fischer Scientific, Waltham, MA) by using a multichannel pipet. Activities of seven EEA were assayed (Table 2), which represented key reactions in carbon, nitrogen and phosphorus cycles; all EEA were measured by using substrates that generated a fluorescent product, either 4-methylumbellyferylone (MUB) or 7-amino-4-methylcoumarin (AMC). All substrates were obtained from Sigma Aldrich Chemicals (St. Louis, MO). Stocks (20 mM) of all substrates were made in dimethyl sulfoxide (DMSO, Sigma Aldrich) except 4-MUB-β-D-xylopyranoside, which was dissolved in pyridine (Sigma Aldrich). The working solutions were created by diluting all stocks to 200 μM (in DMSO or pyridine as appropriate), 50 μl of which was added to each assay giving a final substrate concentration of 40 μM. Fluorescence was measured at time zero, and after 3 h incubation in the dark at 30°C by using a Synergy plate reader (BioTek Instruments, Winooski, VT) with excitation = 365 nm, emission = 450 nm, and the chamber held at 30°C. Product formation was quantified by reference to six point standard curves (100 μM to 1 μM) of MUB or AMC in the assay buffer. The MUB and AMC standard curves were also established in 200 μl aliquots of the soil suspensions as quench controls. Blanks included: assay buffer alone, substrates in buffer alone, and soil suspension in assay buffer alone. For each soil sample, the 12 replicate fluorescence values obtained for each substrate were averaged, corrected for quenching and background, and rates of product formation calculated that were normalized to oven dry weight of soil (nmol/h/g soil).

Phospholipid fatty acids (PLFA) were extracted from soil by using a Bligh and Dyer technique (Bligh and Dyer, 1959) as modified as described by (Moeskops et al., 2010). Briefly, each soil sample (2 g) was extracted with a chloroform–methanol–citrate buffer mixture (1:2:0.8 v/v/v), and then fractionated into neutral, glycol, and phospholipids by passage over a silicic acid column (SPE-SI; Bond Elut 3CC, Varian, Inc.). Phospholipids were then subjected to mild alkaline methanolysis for esterification of fatty acids. The extracts were analyzed by gas chromatography-mass spectrometry by using a Hewlett–Packard 6890 fitted with a DB-5 ms column (Agilent, Santa Clara, CA) interfaced to a Hewlett–Packard 5973 mass selective detector. Peak identification was done referencing retention times and mass spectra; quantification was done by using 19:0 as an internal standard. A total of 86 PLFA were identified, and up to 51 of these with carbon chains ≤C20 were used for biomass quantification. Selected PLFA were used to quantify specific guilds. Gram-positive (GP) bacteria were quantified by the sum of 15:0 ISO, 15:0 ANTEISO, 17:0 ISO, and 17:0 ANTEISO. Gram negative (GN) bacteria were quantified by the sum of PLFA 16:1ω7c, 17:0 CYCLO, 18:1 ω5c, 18:1 ω7c, 19:0 CYCLO and 19:0 CYCLO 11-12 2OH. The sum of 16:0 10 Methyl and 18:0 10 Methyl was used to quantify Actinomycetes. Fungal quantification was done with the sum of 16:1 ω5c, 18:1 ω9c, and 18:2 ω6,9c with the former used to quantify the subset of arbuscular mycorrhizal fungi (AMF). The fungi:bacteria ratio calculation was based on the sum of all fungal PLFA listed above, divided by the sum of the above-listed PLFA used for quantification of GN- and GP-bacteria (Joergensen and Wichern, 2008).

Three DNA extractions were prepared from each soil in batches of 0.25 g, and a 1 uL aliquot from each of the 0.75 g-equivalent extracts was used for generation of bacterial and fungal amplicon libraries. Bacterial libraries were created with universal prokaryote primers 515F and 806R that targeted the V4–V5 region of the 16S rRNA gene (Caporaso et al., 2010), amplicons were generated following PCR protocols described in the work referenced as the source of primers. Fungal libraries were created with primers ITS1F (Gardes and Bruns, 1993) and ITS4 (White et al., 1990) following the PCR protocol described by Manter and Vivanco Manter and Vivanco (2007). Region specific primers were modified to add Illumina adapter overhang nucleotide sequences to the amplicons. Following initial amplification, library size was verified on an Agilent DNA1000 chip, cleaned using a 1X volume of Mag PCR clean-up beads (Axygen Biosciences, Union City, CA) and then Illumina dual indexes and sequencing adapters were added by PCR (de Gannes et al., 2015). Following PCR, samples were cleaned and normalized by using a Sequal Prep Normalization Plate (Life Technologies, Carlsbad, CA). Quality and quantity of the libraries were assessed by using an Agilent DNA1000 chip and Qubit® dsDNA HS Assay Kit, respectively, and were standardized to 2 nM prior to pooling and sequencing. Sequencing was done with an Illumina MiSeq system (Illumina, San Diego, CA) by using Miseq reagent kit v. 3 (Illumina) to generate 2 × 250 bp paired end reads. Images were analyzed using the standard Illumina Pipeline, version 1.8.2.

Illumina datasets were de-multiplexed by using MiSeq Reporter v2.2.31. (Illumina) with a Q20 minimum value as a quality filter, and then reads trimmed of forward and reverse primers by using cutadapt (Martin, 2011). For bacterial libraries, paired sequences were merged into single reads by using FLASH (Magoč and Salzberg, 2011), and then length-filtered. The fungal ITS extractor was applied to fungal libraries to isolate ITS sequences by removing adjoining regions encoding ribosomal RNA (Nilsson et al., 2010). The QIIME (Quantitative Insights into Microbial Ecology) package (v. 1.8.0) was then utilized for generation of operational taxonomic units (OTU) generation by picking with against Greengenes (v. 2013_08), and supplemented by de novo clustering at 97% similarity. The sequence libraries generated by the forward and reverse ITS primers were maintained as separate data sets. The fungal community analyses was done using sequences from the forward primer, which covered ITS1 since there has been a larger use of that region in the literature (Chen et al., 2001; Narutaki et al., 2002; Hinrikson et al., 2005; Nilsson et al., 2008).

For bacterial libraries, taxonomy was assigned by BLAST against GreenGenes (v. 2013_08) and the libraries were screened for chimeric reads by using Chimeraslayer against GreenGenes (v. 2013_08). For fungal libraries, the UNITE database (alpha version 2012_11) was used for taxonomic assignment of OTU by BLAST, and for chimera screening. The QIIME package (v. 1.8.0) was used for OTU-based analyses including rarefaction and computation of diversity metrics (excluding singletons). For comparisons between samples, libraries were rarefied to a common number of sequences for the different taxonomic groups, which corresponded to the smallest library (Figure S1). Prism 6 (Graphpad, La Jolla, CA) was used to display the composition of each library by using average relative abundance for taxa believed to be best able to convey biologically relevant information.

The R programming environment (www.r-project.org) was used for univariate analyses including analysis of variance (ANOVA), and pairwise comparison (Tukey's test). False discovery rates were determined by the method of Benjamini and Hochberg (1995) and significance was assessed by setting the resulting q-values ≥ 0.05. Multivariate analyses were done by using Primer-E v. 6 (PRIMER-E Ltd, Lutton, UK). For Illumina libraries, the numbers of reads assigned to an OTU were normalized to the total reads in a given library. Each dataset was square root transformed to down weight effects of highly abundant OTU and then a resemblance matrix constructed based on Bray-Curtis similarities (Clarke et al., 2008). Bray-Curtis similarity matrices were also constructed on PLFA data that was square root transformed to down weight importance of highly abundant individual PLFA. The Illumina and PLFA matrices were analyzed by non-metric multidimensional scaling (NMDS) ordination (25 restarts, minimum Kruskal stress = 0.01). Similarity contours overlaid on NMDS ordinations were based on groupings developed with the CLUSTER routine. The matrices were analyzed by two-way analysis of similarity (ANOSIM, tests the null hypothesis: the average of the ranks of within-group distances is greater than or equal to the average of the ranks of between-group distances) using both nested and crossed designs. Nested ANOSIM was first applied to determine if there were significant differences between Teak stands using the factors “stand” (the different teak stands according to age, selected for the study) and “location” (North or South). Crossed analysis examined the data for differences between the Teak plantations as a group compared to the native forest as a group, using the factors “forest type” (either the plantation teak forest or the natural forest) and “location.” Any factors determined to be significant by ANOSIM were then examined by canonical analysis of principal components (CAP) and ordination; CAP trace statistics were determined by 999 permutations. Relationships between beta-diversity patterns and environmental data were examined by using two non-parametric tests, the RELATE and the Bio-environmental Step (BEST) routines. RELATE is a Mantel-type test that generates a Spearman Rho test statistic that measures congruence between matrices of biotic or abiotic data. BEST identified subsets of physicochemical variables that gave rank order similarities (Euclidean distance) between soils that best matched the rank order Bray-Curtis similarities of microbial community composition (Clarke et al., 2008).

The Illumina sequence data reported here is deposited in the NCBI Sequence Read Archive (http://www.ncbi.nlm.nih.gov/sra) under accession number SRP058555.

Levels of β-xylosidase were similar in the north native forest and 57N samples, but were significantly (F = 5.29, p = 0.0163) lower than that of all other samples (Figure 1). At the north location, β-xylosidase activity was highest in the 66N sample, which was significantly greater (F = 4.68, p = 0.0113) than that of 09N (Figure 1). At the south location, β-xylosidase activity did not vary significantly among samples, and levels were on average similar to that of 66N. The north native forest and 57N samples were also similar in having comparatively high levels of phosphomonoesterase (acid phosphatase), which were significantly greater (F = 5.33, p = 0.0313) than those of 09N or 66N samples (Figure 1). Again, in the south samples, there were no significant differences in acid phosphatase activity, and levels were on averaged similar to those in the 09N and 66N samples. For the other four EAA assayed, there were no significant differences in activity in any of the samples.

Biomass levels based on PLFA content differed significantly (ANOVA F = 5.59, p = 0.0192) between the north stands (Table 3). The specific stand pairs from the north that showed significant differences (Tukey test p < 0.05) were: 57N vs. 09N, native vs. 57N and native vs. 09N. While these differences were detected, there were no apparent trends: soil biomass was not consistently higher or lower in Teak vs. native forest, and among Teak stands, biomass was not consistently higher or lower as a function of stand age. Microbial biomass levels in the south location were on average significantly greater than those of the north (Table 3; ANOVA F = 7.43, p = 0.021) but there were no significant differences between south stands in the amounts of microbial biomass (ANOVA F = 1.12, p = 0.511).

For PLFA-based community structure analysis, two broad indicators used were the ratio of fungi to bacteria (F:B) and the ratio of GP bacteria to GN bacteria (GP:GN). Neither of these ratios differed significantly among soils from native forest or Teak stands at either location. The main trend apparent from the data was the much greater dispersion of F:B values in the Teak stands relative to the native forest (Figure 2). For native stands, all F:B were in a relatively narrow range and all were <1. Teak stands had a much wider range with some >1. At the north, there were no significant differences among stands in the abundance of individual guilds assayed by PLFA (Table 3). However, at the south location, there were significant differences between soils in the abundance of Actinomycetes (F = 8.01, p = 0.0086), fungi (ANOVA F = 6.00, p = 0.0191), and GP bacteria (ANOVA F = 3.95, p = 0.0532). But, as with biomass, while differences were detected, there were no apparent trends as a function of forest type or of the age of Teak stands. The exception was the Actinomycetes, which, at the south location, were significantly less abundant in all Teak soils than in the soils of the native forest (ANOVA F = 6.97, p = 0.002).

Multivariate analysis of community structure based on PLFA data included unconstrained ordination by NMDS. The primary trend illustrated by this analysis was separation of microbial communities by location; there was no apparent clustering by stand type or by stand age in the case of the Teak samples (Figure 3). CLUSTER analysis revealed that the structure of microbial communities was 80% similar across all soils. Two-way nested ANOSIM examining differences among Teak stands showed that only location was significant (stand, p = 0.069; location, p = 0.01) while two-way crossed ANOSIM showed that both forest type and location were significant (forest type, p = 0.044; location, p = 0.01). CAP analysis using forest type as a factor gave strong canonical correlations (δ₁ = 0.89, δ₂ = 0.87) and by permutation yielded a highly significant first squared canonical correlation (δ12 = 0.789, p = 0.004) and a highly significant trace statistic (1.157, p = 0.001). In this constrained ordination, the first canonical axis (CAP1) separated soils by forest type, with native forest soils segregated to the left and Teak soils segregating to the right (Figure S2). The second canonical axis (CAP2) separated samples by location with the north and south samples segregating to the bottom and top of the axis, respectively (Figure S2).

Potential correlations of variation in PLFA-based community structure to edaphic properties were also examined by using the RELATE and BEST tests. The former analysis gave a Rho value of moderate strength (0.413) but which was strongly significant (p = 0.001). The BEST analysis gave a single correlation identified the parameters most closely associated (0.56) but which was strongly significant (p = 0.009) included six variables: TC, TN, NH₄-N, NO₃-N, CEC, and Mn.

Rarefaction curves showed that the sequencing work was relatively comprehensive to cover bacterial and fungal diversity as the rarefaction curves tended to approach saturation plateau (Figures S3, S4).

A total of 879 OTU were identified with six OTU assigned as Archaea and 872 OTU assigned as Bacteria (Table S2). For the bacterial community, 20 taxa accounted for the majority of sequences (Figure 4A). The Acidobacteria predominated across all soil samples (28–51%), followed by Proteobacteria (Gammaproteobacteria, 6–16%; Alphaproteobacteria, 9–13%; Deltaproteobacteria and Betaproteobacteria, (3–7%). Other top taxa were: Actinobacteria (4–15%), Bacteriodetes (4–10%), Verrucomicrobia (4–7%), and Plantomycetes (4–6%). Four phyla differed significantly between the native soils and those in the Teak plantations: those that were significantly more abundant in the native soils than in soils of Teak plantations were Cyanobacteria (ANOVA F = 6.54, p = 0.0180) and Nitrospirae (ANOVA F = 7.96, p = 0.0100) while those that were significantly more abundant in Teak soils than in soils of the native forest were candidate phylum TM7 (ANOVA F = 17.14, p = 0.0004) and candidate phylum WS6 (ANOVA F = 4.54, p = 0.0445).

For bacteria, the top two most abundant OTU were the same in all soils, being identified as Acidobacteria either to the order Ellin6513 (most abundant OTU) or to the family level as Koribacteraceae (second most abundant OTU; Table S2). Besides the top two OTU, the soils were further similar in that same collection of OTU comprised were most abundant in all samples although the specific order of abundance varied by soil. This group included OTU identified as (listed in order of decreasing average abundance across all soils): Acidobacteria order iii1-15, Bacteroidetes family Chitinophagaceae, Alphaproteobacteria family Rhodospirillaceae, Planctomycetes, order WD2101, Acidobacteria genus Candidatus Koribacter, Acidobacteria genus Candidatus Solibacter, Acidobacteria class Acidobacteria-5, Alphaproteobacteria genus Rhodoplanes and Verrucomicrobia genus DA101.

For the fungi, a large group of sequences (49–80%) were assigned only to the Fungi domain level (Table S3; Figure 4B). At the level of class, the Sordariomycetes predominated (3–16%) followed by the Agaricomycetes (2–21%). The Eurotiomycetes were prevalent in two of the south soils 43S, (12%) and 56S (8%). There were no significant differences in abundance of any of the fungal classes, except for one unidentified Glomeromycota class, which was significantly more abundant in Teak soils than in soils of the native forest (ANOVA F = 4.47, p = 0.047). A total of 793 OTU were identified, with one OTU (assigned only as Fungi) dominant across all soils comprising 29–79% of the total reads in a given library. As with the bacteria, the fungal libraries of all soils were composed of a few common OTU, which were identified as (listed in order of decreasing average abundance across all soils): Basidiomycota genus Hygrocybe, Ascomycota class Sordariomycetes, Basidiomycota genus Sparassis and as the phylum Ascomycota.

Alpha-diversity parameters (OTU richness and evenness) did not differ significantly (p = 0.748) between the bacterial communities as reflected in the Chao1 metric. Alpha-diversity characteristics were visualized with Whittaker plots, where line slopes reflected species (OTU) evenness and line lengths indicated OTU richness, which were similar for all soils (Figures 5A,B). The fungal communities also showed no significant differences (p = 0.627) in Chao1. Likewise, Whittaker plots displayed similar patterns of fungal diversity in all of the soils (Figure 5B). The domination of all fungal communities by a single OTU was readily apparent from the Whittaker plots (Figure 5B); this OTU was the same in all soils and was identified only to the kingdom level as Fungi.

NMDS ordination of the bacterial community data revealed a pattern similar to that of the PLFA data, with the prevailing trend being segregation by location and a lack of trends associated with forest type or age within teak stands (Figure 6). CLUSTER analysis showed bacterial communities were at least 80% similar across all soils (Figure 6). Two-way nested ANOSIM revealed no significant differences among Teak soils, while two-way crossed ANOSIM showed that both forest type and location were significant (p = 0.058 and p = 0.006, respectively). CAP analysis using forest type as a factor gave strong canonical correlations (δ₁ = 0.93, δ₂ = 0.86) and by permutation yielded a highly significant δ12 of 0.789 (p = 0.004) and a highly significant trace statistic of 1.157 (p = 0.001). Ordination of the CAP results showed Teak soils clustering near the center of CAP1, while CAP2 separated samples by location with the north and south soils segregating to the bottom and top, respectively (Figure S5). Phyla that characterized the northern samples (Teak and native forest) were Acidobacteria and candidate phylum AD3 (Figure S5). For the south location, Planctomycetes characterized the Teak soils, while Actinobacteria and candidate phylum OP3 were more strongly correlated with the native forest.

The NMDS ordination of the fungal community data was distinct from that of the PLFA data or of the bacterial communities, and showed no apparent trends in the data (Figure 7). CLUSTER analysis showed fungal communities were at least 60% similar across all soils. Two-way nested ANOSIM revealed no significant differences among Teak soils (p = 0.058), while two-way crossed ANOSIM showed that both forest type and location were significant (p = 0.006).

Potential correlations of bacterial and fungal community structure to edaphic properties were examined by using the RELATE and BEST tests. For the bacterial communities, RELATE gave a Rho value having relatively low strength and significance (0.224, Rho = 1 indicates a variation in soil properties exactly matches variation in community structure; p = 0.0588), whereas BEST analysis gave a no significant correlations. Fungal communities showed no correlation to soil properties: RELATE yielded a Rho of −0.206 (p = 0.938) and BEST analysis returned no significant correlations.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.