Root-Associated Microbial Communities of Abies nordmanniana: Insights Into Interactions of Microbial Communities With Antioxidative Enzymes and Plant Growth

Abies nordmanniana is a major Christmas tree species in Europe, but their uneven and prolonged growth slows down their production. By a 16S and 18S rRNA gene amplicon sequencing approach, we performed a characterization of root-associated bacterial and fungal communities for three-year-old A. nordmanniana plants collected from two nurseries in Denmark and Germany and displaying different growth patterns (small versus tall plants). Proteobacteria had the highest relative abundance at both sampling sites and plant sizes, and Ascomycota was the most abundant fungal phylum. At the order level, Acidobacteriales, Actinomycetales, Burkholderiales, Rhizobiales, and Xanthomonadales represented the bacterial core microbiome of A. nordmanniana, independently of the sampling site or plant size, while the fungal core microbiome included members of the Agaricales, Hypocreales, and Pezizales. Principal Coordinate Analysis indicated that both bacterial and fungal communities clustered according to the sampling site pointing to the significance of soil characteristics and climatic conditions for the composition of root-associated microbial communities. Major differences between communities from tall and small plants were a dominance of the potential pathogen Fusarium (Hypocreales) in the small plants from Germany, while Agaricales, that includes reported beneficial ectomycorrhizal fungi, dominated in the tall plants. An evaluation of plant root antioxidative enzyme profiles showed higher levels of the antioxidative enzymes ascorbate peroxidase, peroxidase, and superoxide dismutase in small plants compared to tall plants. We suggest that the higher antioxidative enzyme activities combined with the growth arrest phenotype indicate higher oxidative stress levels in the small plants. Additionally, the correlations between the relative abundances of specific taxa of the microbiome with the plant antioxidative enzyme profiles were established. The main result was that many more bacterial taxa correlated positively than negatively with one or more antioxidative enzyme activity. This may suggest that the ability of bacteria to increase plant antioxidative enzyme defenses is widespread.


INTRODUCTION
The soil ecosystem contains the highest reported microbial diversity on earth (Torsvik et al., 2002) that serves as a resource wherefrom plant roots can recruit specific microbial communities to the rhizosphere (Mendes et al., 2013). The microorganisms in the rhizosphere can be beneficial, harmful or neutral for the growth and health of the plant. For example, plant growth promoting bacteria and fungi may enhance protection against pathogens, promote plant growth (Hardoim et al., 2015), and facilitate plant nutrition . In contrast, plant pathogenic bacteria and fungi represent a threat toward plant health (Plett and Martin, 2018). Plants and microbes interact, e.g., by signaling via root exudates. The composition of root exudates varies between plant species (Zhalnina et al., 2018), and this variability plays an important role for the establishment of plant-rhizospheric microbial communities (Bais et al., 2006;Chaparro et al., 2013).
Abies nordmanniana (Steven) Spach is a major Christmas tree species grown throughout Europe, where the production exceeds 30 million trees annually (Bräuner Nielsen et al., 2011). However, the natural growth of this species is slow, and trees only reach commercial use after several years. The ability of A. nordmanniana to select specific rhizosphere microbial communities has so far not been addressed; nor have these communities been characterized, although they represent a possible source for plant growth and health promoting microorganisms. Many conifer species develop associations with beneficial microorganisms; in particular ectomycorrhizal fungi (EM), which enhance plant nutrition and growth (Marx, 1970;Ważny, 2014;Rudawska et al., 2016). Similarly, the study of Zulueta-Rodriguez et al. (2015) demonstrated positive effects of seed inoculation with known plant growth promoting rhizobacteria (Bacillus subtilis, Pseudomonas fluorescens, and P. putida), previously isolated from other plant species. These inoculants enhance germination and seedling growth of A. hickelii and A. religiosa; two endangered Mexican pine species. On the other hand, several fungal pathogens, e.g., from the genera Armillaria, Fusarium, and Rhizoctonia can negatively influence the health of Abies species (Oliva et al., 2009;Gordon et al., 2015), while EM in symbiosis with Abies may even offer protection against pathogens in seedlings (Ważny, 2014).
In plants, environmental stressors such as drought and high salinity lead to the accumulation of reactive oxygen species (ROS), which can cause severe damage to the cells (Caverzan et al., 2016). To protect themselves, plants launch an antioxidative defense, where antioxidant (ROS scavenging) enzymes play a major role (Gill and Tuteja, 2010). Hence, activities of enzymes such as catalase (CAT), superoxide dismutase (SOD) ascorbate peroxidase (APX) and peroxidase (POX) often increase under abiotic stress (Mhadhbi et al., 2004;Jebara et al., 2005;Kohler et al., 2008;Mandal et al., 2008;Bharti et al., 2016;Caverzan et al., 2016;Sarkar et al., 2018). In addition to their role in coping with ROS formation caused by environmental stresses, antioxidant enzymes are essential for the maintenance of cell redox homeostasis. By regulating ROS deriving from central processes in different cellular compartments such as photosynthesis and respiration, they can contribute to the regulation of plant growth processes (Das et al., 2015).
Microorganisms introduced to the rhizosphere can affect the activity of antioxidative defense pathways in plants. Hence, Bharti et al. (2016) reported that the expression of several defense enzymes was enhanced by an introduced rhizobacterium in wheat under salt stress, and thereby mediated salinity tolerance. However, there are even examples of introduced strains mediating drought tolerance that decrease the levels of antioxidative enzymes in plants (Sandhya et al., 2010;Armada et al., 2016). Moreover, it has been suggested that fungal endophytes might contribute to redox regulation, improving the antioxidant capacity of the plant host by regulating the host genes. Hence, Brotman et al. (2013) reported that Trichoderma induced expression of genes coding for antioxidant enzymes in Arabidopsis thaliana. While the effects of specific, introduced microorganisms on antioxidative plant enzymes have been assessed in inoculated plants and in vitro studies, these interactions, to the best of our knowledge are poorly understood in field-grown plants interacting with natural root-associated microbial assemblages. The different results, however, point at different mechanisms whereby microbes mediate beneficial effects to plants, ranging from stress reduction and avoidance to improved stress tolerance.
A deeper insight into the composition of microbial communities associated with A. nordmanniana and their interactions with plant antioxidative enzymes, and ultimately plant growth phenotype, is highly relevant to facilitate strategies for microbial mediated growth promotion of A. nordmanniana. In this study, we characterized and compared the composition of root-associated bacterial and fungal communities of A. nordmanniana showing different growth habits when grown at two field sites. In parallel, we determined signatures of A. nordmanniana antioxidative enzyme activities in roots. Finally, we established correlations between the relative abundances of specific microbial taxa and the plant antioxidative enzyme profile in the context of the tree growth phenotype.

Sampling Site
Seedlings of A. nordmanniana were collected from two Christmas tree nurseries; sampling site 1: Baumschule Engler located in Hohenlockstedt, Germany (53 • 58 30.5 N 9 • 37 51.9 E), and sampling site 2: Primo Plant Ejendomme ApS located in Hadsund, Denmark (56 • 44 22.7 N 10 • 03 36.7 E). Sampling at both sites was performed during July 2016, the climatic conditions for both sampling sites were accessed using the German National meteorological Service 1 and the Danish National Meteorological Service 2 . During the sampling dates, the Danish site presented a higher maximum temperature of 28 • C, and the German site presented a maximum temperature of 21 • C (Supplementary Table S1). The German site presented the higher variation in rainfall during the 3 months prior to sampling, with values between 21 and 101 mm, and at the Danish site, the differences in rainfall ranged between 41 and 80 mm ( Supplementary Table S1).
Additionally, the climatic data for both sampling sites was also accessed for the 3 years of growth of the sampled A. nordmanniana plants. The yearly average temperature showed comparable maximum and minimum temperatures, and comparable yearly hours of sun for the German and Danish sampling sites, respectively (Supplementary Table S1). However, higher average rainfall values were observed for the Danish site (932 mm) during the 3-year period, compared with the German site (801 mm) (Supplementary Table S1).

Plant Material
Christmas tree production in each nursery derived from the same parental seed source: Berritzgaard F 665, which originates from a small area in Ambrolauri, Georgia. According to the seed providers Levinsen and Abies A/S, the second generation of Berritzgaard plants was raised at the Hildesvig forest in Denmark. Since 1991, the trees at the Berritzgaard plantation have been selected for uniformity in traits of importance for Christmas tree production.
Abies nordmanniana plants (3 years old) displaying different growth were collected at each sampling site, where five plants with an average height of 25 cm, hereafter referred to as "small" plants, and five plants with an average height of 60 cm, hereafter referred to as "tall" plants, were sampled (Supplementary Table S2). Each plant was considered as a biological replicate for plant size and sampling site. Plant heights were measured from the top of the apical bud to the end of the main root. The plants were collected in their entirety with the surrounding bulk soil and transferred into individual pots, in which they remained for 1.5 days at 18 • C in a greenhouse, at the Department of Plant and Environmental Sciences, University of Copenhagen, Frederiksberg, Denmark, before samples for DNA extraction, plant enzyme analysis and soil analysis were taken.

Soil Characterization
In order to perform soil characterization, three soil subsamples from each sampling site were randomly collected with a hand shovel and represented top soil down to 30 cm depth. Soil samples were placed in plastic bags and later combined per sampling site. For soil analysis, 500 g samples were used; soil physico-chemical analyses included the determination of concentrations of macro and micronutrients, soil organic matter and soil texture. The analyses were performed using ISO certified procedures at Eurofins Agro Testing, Denmark, and method specifications can be found at www.eurofins.dk. Elemental analyses employed Inductively Coupled Plasma -Optical Emission Spectrometry (ICP-OES), or Flow Injection Analysis-spectrometry (FIA) as specified in Supplementary  Table S3. Additionally, soil pH was measured by following the protocol of Hanlon (2012).

Rhizosphere and Bulk Soil DNA Extraction
Plants were removed from the pots and the bulk soil was removed from the roots with a brush. Subsequently, entire root systems with primary and secondary roots (∼5-10 mm diameter and ∼15-20 cm length) with firmly adhering soil were cut into short fragments with a sterile scalpel. The fragments were then macerated with liquid nitrogen in a mortar, to obtain finely ground and well-mixed samples representing plant tissue and rhizosphere soil. These root-associated samples are hereafter referred to as rhizosphere samples.
For each plant, 0.25 g of rhizosphere and bulk soil samples were used for DNA extraction, which was carried out using the PowerSoil R DNA Isolation Kit (Mo Bio Laboratories, Carlsbad, CA, United States) according to the manufacturer's instructions. The nucleic acid concentration of each sample was quantified using a NanoDrop 2000c spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, United States) (Desjardins and Conklin, 2010).

Enzyme Activity Signatures of Plant Antioxidative Metabolism
After the recovery of plant material with firmly adhering rhizosphere soil for DNA extraction, the remaining root tissue was washed with tap water to remove any remaining soil particles and was used for the determination of plant physiological parameters. For these determinations, root material was snap frozen and ground in liquid nitrogen and stored at −80 • C until further processing.
Activities of the antioxidative enzymes APX, CAT, cell wall and cytoplasmic peroxidases (cwPOX, POX), glutathione reductase (GR) and SOD were determined photometrically by kinetic assays in a miniaturized 96-well plate format following the approach presented in Jammer et al. (2015). The assays were based on principles published by Yoshimura et al. (2000) for APX, Aebi (1984) for CAT, Polle et al. (1994) for (cwPOX) POX, Edwards et al. (1990) for GR and McCord and Fridovich (1969) for SOD.

16S rRNA Gene Amplicon Sequencing
To confirm the presence of bacterial DNA in the rhizosphere and bulk soil samples, the 16S rRNA gene region was amplified by PCR using the following protocol; 2 µl of extracted DNA was added to a PCR reaction mix containing 10× PCR buffer without salt (Sigma-Aldrich Co.), 25 mM of MgCl 2 (Sigma-Aldrich Co.), 10 mM of each dNTP, 5 U Taq DNA polymerase (Sigma-Aldrich Co.), 10 µM of 341 Forward primer and 10 µM of 806 Reverse primer [primer sequences are in Sundberg et al. (2013)] and molecular grade water (Sigma-Aldrich Co.) up to a final volume of 50 µl. Primers were obtained from TAG Copenhagen A/S, Copenhagen, Denmark. The PCR amplification was conducted in a Mastercycler Pro thermal cycler (Eppendorf, Hauppauge, NY, United States), as follows: an initial denaturation step consisting of 4 min at 94 • C; 35 cycles of 30 s at 94 • C, 1 min at 55 • C, 1 min at 72 • C, and a final extension of 7 min at 72 • C. The PCR products were visualized using a gel documentation system (Bio-Rad Gel Doc TM 2000 Bio-Rad Laboratories CO) after agarose gel electrophoresis and gel staining with GelRed TM Nucleic Acid Gel Stain (Biotium Inc., United States).
Sequencing data was imported to CLC Genomics Workbench Version 8.0 3 . The raw data was processed following the CLC Microbial Genomics Module pipeline. Briefly, sequences were paired, then, the primer sequences were trimmed off with the following parameters: removal of low-quality sequence (limit = 0.05), removal of ambiguous nucleotides: maximal two nucleotides allowed. The removal of adapter sequences was performed using the following adapters: Bakt_805R (GACTACHVGGGTATCTAATCC), strand = Plus, action = Remove adapter, score = [2: 3: 8: 8], and Bakt_341F (CCTACGGGNGGCWGCAG), strand = Plus, action = Remove adapter, score = [2: 3: 8: 8], removal of sequences on length: minimum length five nucleotides. Next, the paired read sequences were merged, and the length was fixed by trimming to an exact length of 402 bp.
Samples were filtered based on the number of reads, with a minimum number of reads set to 100, discarding the samples that did not fulfill this parameter. Subsequently, the filtered sequences with comparable length and coverage were clustered and assigned to operational taxonomic units (OTUs). A 97% of similarity was used when comparing with the Greengenes 16S rRNA gene database. A list of sequences that clustered with sequences in the Greengenes database was obtained and used to generate the abundance OTU table. Furthermore, chimeras were detected and removed in this step. Subsequently, an OTU relative abundance table was generated by summing up all OTU counts and dividing each count with the total sum.

18S rRNA Gene Amplification and Sequencing
Isolated DNA from rhizosphere samples was PCR amplified using one forward and two reverse primers for multiplex PCR targeting two universally conserved domains. The assay amplified a variable region (of ∼400 bp) from the small subunit of 18S rRNA gene of plant-associated microbial eukaryotes: basidiomycetes, ascomycetes, glomeromycetes, zygomycetes, oomycetes, plasmodiophorids and nematodes. Primers were modified after Stokholm et al. (2016) to match the requirements for amplicon sequencing using the Illumina MiSeq platform. The 18S rRNA gene reverse primers contain at least one mismatch relative to sequence accessions representing Pinaceae species. Primers were obtained from Eurofins Genomics, Ebersberg, Germany. The forward primers were fused to a sequencing linker overhang, followed by 10 bp different internal multiplex identifier (MID), tags for barcoding of different samples, 2 × NN, and a variable heterogeneity spacer from 0 to 7 bp to generate primers of different lengths according to Fadrosh et al. (2014; Supplementary Table S4). The reverse primers consisted of a linker overhang, followed by internal MID tags, 2 × NN and a 12 bp insert in front of the reverse primer.
Equimolar concentrations of the 18S rRNA gene amplicons, were mixed and sequenced in parallel at the Illumina MiSeq platform at the Danish Technological Institute, Department of Food Security, Taastrup, Denmark. Next a de-multiplexing of the sequences using sabre (a barcode de-multiplexing and trimming tool for FastQ files 4 ) was made in a virtual Linux environment, according to the internal MID tags. The resulting de-multiplexed paired-end FastQ files were then individually overlap-merged and quality-filtered through a series of sequence trimming and quality control (QC) steps. The quality filtered datasets were then subjected to the open source SCATA pipeline (Sequence Clustering and Analysis of Tagged Amplicons 5 ) at the Swedish University of Agricultural Sciences (Clemmensen et al., 2016), for single linkage clustering into OTUs, using a 98% similarity cut-off and parameters as described in Stokholm et al. (2016). The representative sequences of the OTUs were further analyzed and annotated down to the order level, or to the lowest possible taxonomic level, by comparison to similar sequence accessions in GenBank 6 , as well as by comparison, where possible, to sequences in the SILVA SSU database for 18S sequences 7 (Yilmaz et al., 2014).

Data Analysis
The bacterial and eukaryotic OTU relative abundance tables were used to estimate OTU diversity using PAST version (2.17). The α-diversity, represented by Shannon H and Chao-1 estimations, was calculated in PAST and compared by Student's t-test. The dissimilarity of the microbial communities from bulk soil (bacterial community only) and rhizosphere samples from the two sampling sites, and for small versus tall plants was calculated using the Bray-Curtis index.
The differences in taxa abundance at the phylum level between the sampling sites and plant sizes were compared by Student's t-test using PAST version (2.17). Dissimilarities of the relative abundance to determinate the microbial taxa at order and genus level contributing for the dissimilarity among sampling sites were calculated using the similarity percentage analysis (SIMPER) with the Bray-Curtis dissimilarity index in PAST version (2.17). With the SIMPER analysis, the average Bray-Curtis dissimilarities between all analyzed groups (sampling site and plant size) is decomposed into percentage of contributions of each OTU filtered at the order level, giving the list in decreasing order of contribution of dissimilarity (Clarke and Gorley, 2006). The SIMPER analysis should be considered carefully for taxa with a high variance, as it can fail to properly identify taxa with large between group effects (Warton et al., 2012). In the current analyses, the variances for the taxa contributing to variability between treatments were low; therefore, we have chosen to use SIMPER for identifying the taxa contributing to differences between sampling sites and plant sizes.
Principal Coordinate Analysis (PCoA), based on the Bray-Curtis dissimilarity matrix between rhizosphere samples were performed with the Phyloseq R package version (2.4-3). A permutational multivariate analysis of variance (Adonis), based on the Bray-Curtis distances with 9999 permutations, was used to access the statistical significance of the microbial community between sampling sites, bulk soil versus rhizosphere samples and tall versus small plants. Additionally, a Canonical Correspondence Analysis (CCA) was performed using PAST version (2.17), including the data from the antioxidative enzymatic activities of each root sample, with the bacterial and fungal OTU relative abundance tables filtered at the order level for bacterial communities and order and genus (for Hypocreales) level for fungal communities. Additionally, a oneway ANOVA was used to test the effect of plant size on the fungal relative abundance at order level. Further, the correlation statistic according to Spearman was used to test the significant correlations between the enzymatic variables and the OTU tables at the order or genus level. 6 http://www.ncbi.nlm.nih.gov/ 7 http://www.arb-silva.de/browser/ Venn diagrams were constructed using the software Venny (Oliveros, 2006), to determinate the unique and shared OTUs between each sampling site and plant size. For this, the filtered list of OTUs at the order level for both bacterial and fungal communities, were used. For each constructed Venn diagram, the mean relative abundance of the shared and unique OTUs at the order level was used to construct stacked bar plots.
The 16S and 18S rRNA gene amplicon sequences were deposited in NCBI's Sequence Read Archive under BioProject PRJNA515250.

Soil Characteristics
Soil from the two sampling sites were both sandy loams with 5-7% clay and 84-90% sand. Overall soil from both sampling sites presented comparable concentrations of nitrogen and organic matter (Supplementary Table S3). More than twofold differences in chemical composition were observed for manganese and boron content, which was higher in soil derived from the sampling site in Germany compared with the sampling site in Denmark. Danish bulk soil samples had higher amounts of potassium, magnesium, copper and iron (Supplementary Table  S3). Soil pH values were comparable, with values of 5.8 and 5.5 for Germany and Denmark, respectively.

Root Antioxidative Enzyme Activity Signatures
Abies nordmanniana plants originating from the German site showed higher activities of the determined antioxidative enzymes in roots (Figure 1), except for APX, which showed similar levels in plants from both sampling sites (Figure 1A). At both sites, root samples from small plants exhibited higher APX activities than the tall plants; the difference was, however, only statistically significant (P < 0.01) for the plants from the German site. Similarly, POX ( Figure 1B) and SOD ( Figure 1C) activities were higher in roots of small plants at both sites, with statistically significant differences (P < 0.05) for POX in plants from the Danish site, and for SOD in plants derived from the German site, respectively. The activities of GR (Figure 1D), CAT ( Figure 1E) and cwPOX ( Figure 1F) were higher by trend in the roots from small plants in Germany. In contrast, roots of the plants from the Danish site did not show consistent differences in these enzyme activities. Overall, the roots from small plants at the German site showed generally higher activity levels of antioxidative enzymes than roots from tall plants, while the differences between tall and small plant enzyme activities were lower at the Danish site.

Alpha Diversity of Bacterial and Fungal Communities
For bacterial diversity analysis, a total of 3,917,647 16S rRNA gene sequence reads were obtained, and after quality filtering and chimera removal, a total of 1,631,616 sequence reads were used for the analyses. The sequence reads clustered into 28,275 OTUs at 97% sequence similarity level. The rarefaction curves Values represent mean ± SEM derived from five biological replicates of three technical replicates each. * , * * indicate statistically significant differences at the P < 0.05 and P < 0.01 level, respectively, between small and tall plants based on two-sided, unpaired Students t-test.
indicated that the number of OTUs increased with the number of sequences read; hence, a deeper sequencing effort is needed to cover the complete diversity of the bacterial communities in the current environments (Supplementary Figure S1). When comparing the rhizosphere bacterial communities from the two sampling sites, a significantly (P < 0.05) higher evenness was found for the German versus the Danish rhizosphere samples, while the richness, diversity and evenness indices did not differ between plant sizes ( Table 1).
After sequencing of the 18S rRNA gene region, a total of 62,078 paired-end reads were obtained, of which 58,653 could be de-multiplexed using sabre according to the internal barcode tags, while 3,425 pairs of FastQ records did not match the barcodes. The individually qualityfiltered datasets using SCATA resulted in a range of 40-60% accepted reads per sample (average = 54%), with a total of 31,820 reads passing QC. The cluster analysis resulted in 516 OTUs (total reads = 22,864) and 8,956 singletons. The rarefaction curves showed that the sequencing effort was not deep enough to cover the complete diversity of fungal communities in the environments under study (Supplementary Figure S2). A comparison of the richness, diversity and evenness indices for fungal rhizosphere communities between the sampling sites, showed no significant differences. However, a significantly higher (P < 0.05) fungal community richness (Chao-1) was found for the small plants at the German site, and significantly higher (P < 0.05) fungal community diversity (Shannon)

Variation in the Structure of A. nordmanniana Rhizosphere Microbial Communities
The differences in microbial community compositions between rhizosphere samples were measured using a PCoA, based on the relative abundance of OTUs. The PCoA analysis indicated that both bacterial and fungal communities clustered according to the sampling site (Figures 2A,B). Hence, the microbial rhizosphere communities from Germany and Denmark were dissimilar. The non-parametric multivariate analysis performed using the Adonis test at the OTU level, also revealed significant differences in both bacterial and fungal rhizosphere community composition between the sampling sites (Supplementary Table S5). Neither the bacterial nor the fungal community compositions differed significantly between tall and small plants for any of the two sites according to the dissimilarity analysis (Supplementary Table S5).

Bacterial Community Composition
Across both sampling sites, an overall of 24 bacterial phyla were identified. Proteobacteria was significantly more abundant than any other phylum (P < 0.05), with relative abundances of 27-50%, but the relative abundance was not significantly different between sampling sites. Actinobacteria, Acidobacteria as well as Bacteriodetes were also abundant with relative abundances above 20% (Supplementary Figure S3). At sampling site level, the relative abundance of Actinobacteria was significantly higher (P < 0.05) in the Danish samples, while the relative abundance of Bacteriodetes, Firmicutes, and Verrucomicrobia were significantly higher (P < 0.05) in the German samples. At the bacterial order level, the relative abundances of Actinomycetales, Burkholderiales, and Rhizobiales, were higher in the rhizosphere of plants collected in Denmark as compared to Germany (P < 0.05). Nevertheless, they constituted the three most abundant orders among all samples. Furthermore, the relative abundances of Acidobacteriales, Rhodospirillales and Xanthomonadales, were significantly (P < 0.05) higher in rhizosphere from plants collected in Germany (Figure 3). A SIMPER analysis identified Acidobacteriales, Actinomycetales, Burkholderiales, Rhizobiales, and Xanthomonadales as the taxa that accounted for 46% of the dissimilarities between the Danish and German samples (Supplementary Table S6).
At the plant size level, Acidobacteriales and Xanthomonadales were particularly more abundant (P < 0.05) in the rhizosphere from tall plants in Danish samples, compared with the corresponding small plants. Additionally, the relative abundances of Burkholderiales and Rhizobiales were higher in the small plants from both sampling sites, however, the differences were not statistically significant.
The A. nordmanniana bacterial core microbiome constituted 49.1% of the identified OTUs at the order level, which were shared across the two sampling sites and across plant sizes (Figure 4); it was dominated by OTUs representing the Acidobacteriales, Actinomycetales, Burkholderiales, Rhizobiales, and Xanthomonadales (Supplementary Figure S4). The German samples presented a higher percentage of unique OTUs 28.3% compared with the Danish samples (8.7%, Figure 4). For the OTUs exclusively found in the microbiomes of small plants there was little overlap between German and Danish samples and the same was the case for OTUs exclusively found in tall plant microbiomes (Figure 4 and Supplementary Figure S5); overall, these unique taxa presented low relative abundances.
Correlations between activities of antioxidative enzymes and the relative abundances of the taxa at the order level, which contributed at least 3% to the variation between samples according to SIMPER (Supplementary Table S6), were determined by a CCA. The analyses showed correlations for several taxa to one or more enzyme activity. Hence, relative abundances of Actinomycetales, Rhizobiales, Sphingobacteriales, and Xanthomonadales were positively correlated with APX (P < 0.05); abundances of Acidobacteriales, Burkholderiales, Rhodospirillales, Saprospirales and Xanthomonadales were positively correlated with SOD (P < 0.05); and Acidobacteriales, Caulobacterales, Saprospirales, and Sphingomonadales were positively correlated with POX (P < 0.001). In contrast only Burkholderiales was negatively correlated with APX (P < 0.01) and Actinomycetales and Rhizobiales with SOD (P < 0.05). Finally, Rickettsiales did not show any significant correlation with the analyzed enzymes ( Figure 5A).

Fungal Community Composition
The 18S rRNA amplicon-targeted analysis identified three fungal phyla, Ascomycota, Basidiomycota and Murocomycota. OTUs that could be assigned to oomycetes, plasmodiophorids or nematodes represented less than 3% the total DNA amplified and were consequently left out of further analyses.
Overall, Ascomycota was more abundant in the rhizosphere of plants collected from Denmark, with a relative abundance of 75-90%, while in German samples, Ascomycota showed a relative abundance of 45-60%, but these differences were not statistically significant. Additionally, for the German samples, Basidiomycota was found in significantly (P < 0.05) higher relative abundance (30-50%) compared with the Danish samples (10-20%; Supplementary Figure S6). The relative abundance of Mucoromycota was also higher compared with the samples from Denmark but these differences were not statistically significant; nevertheless, this phylum had the lowest overall relative abundance at both sampling sites. When comparing the relative abundances of fungal phyla at plant size level, Ascomycota dominated in tall and small plants from Denmark and in the small plants from Germany. For the tall plants in Germany, Basidiomycota was the dominating phylum; however, this dominance was not significant (Supplementary Figure S6).
A comparison of fungal community composition at the order level between the sites showed that Agaricales was significantly (P < 0.05) more abundant in rhizosphere samples from Germany, while Hypocreales and Pezizales was significantly (both P < 0.05) more abundant in the Danish rhizosphere samples (Figure 6). SIMPER analysis identified Agaricales, Hypocreales and Pezizales as the taxa that each accounted for at least 20% of the dissimilarities between the Danish and German samples (Supplementary Table S7). Further, the order Hypocreales was dominated by the genus Fusarium. Hence, according to SIMPER analysis, Fusarium was the genus that contributed the most (21%) to the dissimilarity of fungal communities between the sampling sites (Supplementary Table  S8). Plant size had a significant effect on the relative abundance of Agaricales (P < 0.05), which was higher in the tall plants when comparing with the small plants at both sampling sites. Contrary the order Hypocreales had a significantly (P < 0.05) higher relative abundance in the Danish small plants than in Danish tall plants (Figure 6).
The core fungal microbiome constituted 87.5% of the OTUs identified (order level) (Figure 7) were shared across the sampling sites and plant sizes and was dominated by Agaricales, Hypocreales and Pezizales (Supplementary Figures S7A,B). We did not find OTUs that were unique for small or tall plant microbiomes (Supplementary Figure S7B).
A search for correlations between the relative abundance of fungi at the order level (within the Hypocreales; two genera Trichoderma and Fusarium could be analyzed independently), and activities of antioxidative enzymes by CCA, revealed that the relative abundance of Fusarium (Hypocreales) was significantly (P < 0.01) negatively correlated to APX and SOD. The relative abundance of Agaricales was significantly (P < 0.05) positively correlated to POX. Finally, the relative abundance of the genus Trichoderma (Hypocreales) was significantly (P < 0.01) positively correlated to SOD. The fungal orders Helotiales, Mortierellales, Trechisporales, Pezizales and Xilariales, did not show any significant correlations to the analyzed antioxidative enzymes (Figure 5B).

DISCUSSION
The current study of the bacterial and fungal communities naturally associated with the rhizosphere of A. nordmanniana is, to our knowledge, the first study documenting the rootassociated microbial diversity in relation to plant growth for this Christmas tree species. The most abundant rhizosphere bacteria phyla at both sampling sites (Acidobacteria, Actinobacteria, Bacteriodetes, and Proteobacteria), have been reported before as dominating bacterial communities associated with conifers (Yu et al., 2013;Naumova et al., 2015;Shi et al., 2015;Proença et al., 2016;Uroz et al., 2016). For the fungal communities, Ascomycota and Basidiomycota were the most abundant and common fungal phyla. Our findings are in agreement with previous studies that reported Ascomycota to be associated with conifers (Yu et al., 2013;Stenstrom et al., 2014), while Basidiomycota appears to be the major fungal phylum in the Pinus tebulaeformis rhizosphere (Yu et al., 2013).

Root-Associated Core Microbiome of A. nordmanniana
The A. nordmanniana fungal core microbiome comprised taxa within the orders of Agaricales (Basidiomycota), Hypocreales (Ascomycota) and Pezizales (Ascomycota). Ascomycetes might have roles in nutrient acquisition and transformation of organic nitrogen and phosphorus into their inorganic forms (Stenstrom et al., 2014), and therefore may contribute significantly to the nutrient supply needed for maintaining plant growth, but the order Hypocreales even contains important plant pathogens (Dean et al., 2012). Members of the order of Pezizales might both have biotrophic and decomposing potential (Vrålstad et al., 1998). Finally, the high abundance of Agaricales, which comprises many EM species (Tedersoo and Smith, 2013), points toward ectomycorrhization of A. nordmanninana in the studied field sites, and ectomycorrhization in Abies spp. has been previously reported (Taylor et al., 2000;Rudawska et al., 2016;Oros-Ortega et al., 2017). However, in our study, and for Agaricales, it was not possible to assign OTUs to lower taxa to confirm ectomycorrhization of our root samples.
The bacterial core microbiome of A. nordmanniana was dominated by the orders Burkholderiales and Rhizobiales, but even included Acidobacteriales, Actinomycetales and Xanthomonadales. Non-nodulating Rhizobium spp. and some members of Burkholderiales, as well as of Xanthomonadales, and Actinomycetales are dominantly detected on the ectomycorrhizal (EM) root tips of many conifer species (Khetmalas et al., 2002;Izumi et al., 2006;Burke et al., 2008;Nguyen and Bruns, 2015;Lladó et al., 2017). Therefore, the high abundance of these taxa associated with roots of A. nordmanniana could be influenced by their interactions with EM (Nguyen and Bruns, 2015). Additionally, Burkholderiales and Rhizobiales, includes nitrogen-fixing and mineral-weathering bacteria (Suna et al., 2014) and are part of the most abundant root-associated bacterial orders and core microbiomes for a wide range of plant hosts (Yeoh et al., 2017;Garrido-Oter et al., 2018).
Root exudates have been reported to play an important role in shaping the root-associated microbial communities (Chaparro et al., 2013). Hence, the ability of A. nordmanniana to maintain a core microbiome may be ascribed to specific nutrients and signals in these exudates. Moreover, the possible intraspecific plant genetic variation can alter root-associated microbial communities, thus the interaction of the plant genotype and the soil environment might contribute to the root-microbiome assembly and their complexity in natural environments (Lundberg et al., 2012;Wagner et al., 2016). Furthermore, even different soil characteristics and other factors like climate influence the rhizosphere microbiome associated with a particular plant species (Costa et al., 2006). For the current study including two sandy soils, we observed different relative abundances of microbial taxa between the two sites, as well as a site effect on the bacterial species evenness. These results support the notion that the soil type is an important factor shaping the composition of root-associated microbial communities as seen for other plant species (Mendes et al., 2013;Schreiter et al., 2014;Vanderkoornhuyse et al., 2015). Furthermore, the two sampling sites presented differences in the rainfall during the period prior to sampling. This could have an influence on the microbial communities associated with the sampled rhizospheres, as soil moisture is a key environmental component shaping soil and rhizosphere microbial communities. For example, periods of drought and rewetting could alter the microbial community structure (Cruz-Martínez et al., 2012;Barnard et al., 2013;Chodak et al., 2014).

A. nordmanniana Growth Phenotype
We analyzed the composition of rhizosphere microbial communities for small (growth-retarded) versus tall plants to look for biotic causes or consequences of the growth retardation. The microbiome of small, growth-retarded plants was generally characterized by a higher relative abundance of Fusarium spp., and a lower relative abundance of Agaricales compared to tall plants. The dissimilarities are even reflected in significant differences at plant size level in the fungal community richness, diversity and evenness. Fusarium is a diverse genus including plant pathogenic species (Dean et al., 2012), and many conifers such as Abies spp. are affected by Fusarium infections (Gordon et al., 2015), while some Fusarium species have been described as potential beneficial endophytes, with insecticide and nematocide properties (Vu et al., 2006;Paparu et al., 2008). Associations of Fusarium spp. with roots of conifers have been reported before (Lilja et al., 1997;Stenstrom et al.,  2014). A generally higher abundance of Fusarium spp. has previously been found in the rhizosphere of diseased clonal Picea mariana plants compared to healthy plants (Vujanovic et al., 2007), and these authors indicated that plant health was influenced by the composition of taxa among diseased and healthy P. mariana. The plants we sampled in the current study seemed healthy despite their size difference. However, the higher abundance of Fusarium in the small plants might indicate that some of these fungi could represent latent pathogens, with a potential to cause stress in the small plants. Nevertheless, the identified Fusarium OTUs could also represent part of the fungal biota in the root-attached soil, and not necessary represent latent pathogens inside the sampled roots. Therefore, future studies aiming to a deeper understanding of the interactions between Fusarium spp. and A. nordmanniana, should distinguish the fungal microbiota associated to root tissues, from the one in the rhizosphere soil. Additionally, such studies should employ sequence analysis targeting the ITS region and include a cultivation approach to obtain fungal isolates for detailed investigations.
Alternatively, the order Agaricales that showed a higher relative abundance in tall plants, contains several potential beneficial plant-growth promoting and EM-forming species (Tedersoo and Smith, 2013;Contreras-Cornejo et al., 2015) that may counteract biotic and abiotic stressors. Hence, the tall plants may experience less stress or better nutrient uptake than the small plants, due to beneficial effects by Agaricales spp. (Contreras-Cornejo et al., 2015;Floudas et al., 2015;Rudawska et al., 2016), thus allowing for more resources being allocated to growth. It has been suggested that some fungal species associated with healthy plants could have the potential to compete against root pathogens in the rhizosphere of Picea spp., thus minimizing their pathogenic effect or/and reducing their abundance in healthy plants (Vujanovic et al., 2007).
Hence, an increased relative abundance of potential latent pathogens included in Fusarium spp., as well as a decreased  relative abundance of potentially EM forming members of Agaricales, could be important for the observed plant size differences. Further, the relative abundance of the bacterial order Xanthomonadales was lower in small, growth-retarded plants. Xanthomonadales contain taxa that can cause diseases in conifers, but are even associated with EM (Khetmalas et al., 2002). Their lower abundance in small plants, may be related to the low abundance of Agaricales. We suggest that bacterial disease could not be the major driver for plant size differences.

Root Antioxidative Enzyme Activity Profiles of A. nordmanniana
Results on antioxidative enzyme activities may contribute to the understanding of the plant growth phenotype, due to their essential involvement in regulating growth-related processes and stress responses (Apel and Hirt, 2004;Das et al., 2015). It has previously been reported that an increase in the antioxidative enzyme levels in plants occurs as a mechanism to reduce oxidative stress caused by biotic or abiotic factors (Bharti et al., 2016;Caverzan et al., 2016;Sarkar et al., 2018), and that plants lower their growth in response to stress (Zhu, 2001).
Interestingly, the plants in Germany showed higher antioxidative enzyme activities than plants at the Danish site. Activity profiles of antioxidative enzymes integrate the plant response to abiotic and/or biotic stressors experienced over a longer time than transcriptome profiles (Rasmussen et al., 2013). Hence, the different antioxidative enzymes levels for plants from the two sampling sites, points to differential exposure to biotic and/or abiotic factors prior harvesting. Considering the highly variable precipitation recorded for the German sampling site, we propose that transient exposure to drought, which increases ROS levels in plants (Mittler, 2002), could contribute to the high antioxidative enzyme activities. Moreover, the generally higher enzyme activities observed in small plants, particularly at the German site, could be indicative for increased ROS levels. Subsequently, an increase in ROS levels may demand more plant resources, to cope with the stress and thus limiting plant growth (Zhu, 2001;Caverzan et al., 2016). With the current approach, however, the specific stress source(s) and detailed mechanisms whereby the plants cope with the stress and/or regulate growth, differentially in the tall and small plants could not be resolved.

Abies nordmanniana Root-Associated Microbiome Correlation With Root Antioxidative Enzyme Activity Profiles
We observed that many more bacterial taxa correlated positively than negatively with one or more antioxidative enzyme activity. The ability of bacteria to impact plant antioxidative enzyme activities has primarily been determined for introduced plant beneficial strains and primarily under controlled conditions, where salt or water stress could be imposed on the plants (Jebara et al., 2005;Bharti et al., 2016;Vurukonda et al., 2016). Under such conditions, the introduced bacteria may either increase (Helepciuc et al., 2014;Bharti et al., 2016;Chandra et al., 2018) or decrease the enzyme activities (Vardharajula et al., 2011). The different responses are probably caused by the different mechanisms whereby bacteria can affect plant stress physiology. Some bacteria can increase plant antioxidative enzyme activities directly by interfering with plant stress signaling (Jebara et al., 2005;Vurukonda et al., 2016), while others decrease the activities due to their ability to ameliorate stress, as seen for drought stress by exopolysaccharide production (Vardharajula et al., 2011).
For the fungal taxa, Fusarium correlated negatively with APX and SOD, while Agaricales correlated positively to POX. In contrast to the current finding, an increase in POX activity has been observed at the initial stages of the root necrosis infection with Fusarium spp. in Picea abies and Pinus silvestris (Asiegbu et al., 1999). Studies for other plant species like chickpea, have even documented an increase in antioxidative enzymes during confrontation with Fusarium, as well as a time-dependent enzymatic response reliant on compatibility between host and pathogen (Garcia-Limones et al., 2009;Dalvi et al., 2017). The relation between Agaricales and antioxidative enzyme activities has not previously been addressed for Abies or other conifers. However, our results suggest that the establishment of associations with plant beneficial members of this fungal group, such as EM species (Tedersoo and Smith, 2013), might affect these enzyme activities. For comparison, arbuscular mycorrhizal fungi can both, decrease and increase antioxidative enzyme activities e.g., under drought stress (Armada et al., 2016).
In summary, the current data show a relationship between the composition of rhizosphere microbiota, the plant antioxidative enzymes, and the plant growth habit. They indicate that the microbiota have a strong influence on the cellular metabolism in the roots, and that their ability to increase plant antioxidative enzyme defenses is widespread, or that the plant root recruits beneficial microbes under stressful conditions. Nevertheless, from this field study, we cannot disentangle the causes from the consequences regarding these correlations. However, we suggest that future studies may provide novel knowledge about the possible mechanisms that different microbial groups employ during association with stressed plants. Finally, many of the bacterial and fungal taxa that differ in relative abundance between sites or between tall and small are important for nutrient turnover in soil and rhizosphere. Hence, future studies should focus on relations between the relative abundance of these taxa, soil or plant nutrient concentrations, and the activities of enzymes involved in soil nutrient transformations.

DATA AVAILABILITY
The datasets generated for this study can be found in the PRJNA515250.

AUTHOR CONTRIBUTIONS
AG-L contributed to the experimental design, the sampling, the acquisition and analysis of the data and the interpretation, and the writing of the manuscript. DG contributed to the data acquisition and analysis of plant antioxidative enzymes, and the writing of the manuscript. MS contributed to the eukaryotic gene amplification and amplicon sequencing, the processing of 18S rDNA sequencing data, and the writing of the manuscript. OL contributed to the processing and interpretation of 18S rDNA sequencing data and the revising of the manuscript. MN contributed to experimental design and data analysis, and the revising of the manuscript. TR contributed to the data analysis, and the revising of the manuscript. BV contributed to the experimental design, the sampling, and the revising of the manuscript. ON contributed to the experimental design, the sampling, the analysis of the data and interpretation, and the writing of the manuscript.

ACKNOWLEDGMENTS
We thank Rebecca Andrea Dölker, for her arduous work of processing the root samples for the enzyme profile analysis.