Missing checkerboards? An absence of competitive signal in Alnus-associated ectomycorrhizal fungal communities

A number of recent studies suggest that interspecific competition plays a key role in determining the structure of ectomycorrhizal (ECM) fungal communities. Despite this growing consensus, there has been limited study of ECM fungal community dynamics in abiotically stressful environments, which are often dominated by positive rather than antagonistic interactions. In this study, we examined the ECM fungal communities associated with the host genus Alnus, which live in soils high in both nitrate and acidity. The nature of ECM fungal species interactions (i.e., antagonistic, neutral, or positive) was assessed using taxon co-occurrence and DNA sequence abundance correlational analyses. ECM fungal communities were sampled from root tips or mesh in-growth bags in three monodominant A. rubra plots at a site in Oregon, USA and identified using Illumina-based amplification of the ITS1 gene region. We found a total of 175 ECM fungal taxa; 16 of which were closely related to known Alnus-associated ECM fungi. Contrary to previous studies of ECM fungal communities, taxon co-occurrence analyses on both the total and Alnus-associated ECM datasets indicated that the ECM fungal communities in this system were not structured by interspecific competition. Instead, the co-occurrence patterns were consistent with either random assembly or significant positive interactions. Pair-wise correlational analyses were also more consistent with neutral or positive interactions. Taken together, our results suggest that interspecific competition does not appear to determine the structure of all ECM fungal communities and that abiotic conditions may be important in determining the specific type of interaction occurring among ECM fungi.

buried at each point 5 cm below the soil surface. The bags were made of anti-static polyester fabric with 300 micron diameter pores. This pore size allowed fungal hyphae to grow into the bags, but prevented penetration of plant roots. We filled the bags with twice autoclaved #3 grade Monterey aquarium sand (Cemex, Marina, CA, USA). Aluminum tags on fluorescent string were added to facilitate bag recovery. The mesh bags at Plot 2 were buried on February 1, 2013 and February 22 at Plot 8. They were left undisturbed in the soil until May 31, when all were harvested. After removal from the soil, we placed the mesh bags into individual plastic bags and then onto ice for transport back to the laboratory. Soil cores and bags were stored at 4°C for <96 hours before further processing.

Molecular Analyses
We processed the root tip samples by gently washing all roots away from the soil and removing all ECM colonized root tips from each core under a 10X dissecting scope (~10-50 root tips/core). All roots from each core were extracted using individual MoBio PowerSoil kits (Hercules, CA, USA), following manufacturer's instructions for maximum DNA yields. For the mesh bags, we followed the protocol outlined in Branco et al., (2013), which provided a cheaper and quicker protocol compared to direct DNA extraction from the sand within the mesh bags.
Briefly, each bag (including a negative control that was taken to the field, but not buried) was emptied into a sterile 50 ml centrifuge tube. We added 10 ml of sterile deionized water and vortexed each tube for two minutes, followed by a five minute settling period (hyphae have been previously observed to float to the water surface). We then transferred the top two ml top of water to a new 2 ml centrifuge tube and contents were pelleted via centrifugation. On the same day, we extracted total genomic DNA form the pellets using the Sigma REDExtract-N-Amp kit (Sigma-Aldrich, St, Louis, MO, USA) following manufacturer's instructions. Root tips and extracts were stored for one week at -20°C prior to PCR amplification. For the root tip samples, we combined equal quantity aliquots from all 97 DNA extractions (three cores contained no roots) into a single template for PCR. In contrast, for each mesh bag sample as well as extraction controls, we conducted individual PCR reactions. We processed these two types of samples differently because we were primarily interested in the spatial co-occurrence patterns in the soil ECM fungal communities and therefore only used the root tip samples to create a local sequence reference set of known Alnus-associated ECM taxa against which the mesh bag data could be compared. For all PCR reactions, we used the barcoded ITS1F and ITS2 primer set of Smith & Peay (2014), with each sample run in triplicate and pooled to minimize heterogeneity. Successful PCR products were determined by gel electrophoresis and magnetically cleaned using the Agencourt AMPure XP kit (Beckman Coulter, Brea, CA, USA) according to manufacturer's instructions. Final product concentrations were quantified using a Qubit dsDNA HS Fluorometer (Life Technologies, Carlsbad, CA, USA). Root tip and bag samples were run at different sequencing facilities under the same general conditions. For the root tips, the single PCR product was run at the University of Minnesota Genomics Center using 250 bp paired-end sequencing on the MiSeq Illumina platform. For the bags, we pooled the 192 successfully amplified bag samples at equimolar concentration and ran them on the same platform at the Stanford Functional Genomics Facility using 250 bp paired-end sequencing on the MiSeq Illumina platform. A spike of 20% and 30% PhiX was added to the runs to achieve sufficient sample heterogeneity, respectively. Raw sequence data and associated metadata from both the root tip and bag samples were deposited at MG-RAST (http://metagenomics.anl.gov/) under project #1080.

Bioinformatic Analyses
We used the software packages QIIME (Caporaso et al., 2010) and MOTHUR (Schloss et al., 2009) to process the sample sequences. Raw sequences were demultiplexed, quality filtered using Phred = 20, trimmed to 178 base pairs, and ends were paired, followed by filtering out of sequences that had any ambiguous bases or a homopolymer run of 9 bp. Following the guidelines discussed in Nguyen et al., (in press), we employed a multi-step operational taxonomic unit (OTU) picking strategy by first clustering with reference USEARCH (including de novo chimera checking) at 97% sequence similarity, followed by UCLUST at 97% sequence similarity. We used a 97% similarity threshold because it the most commonly employed in community-level ECM fungal studies, although some lineages, including Alnicola, may have greater sequence similarity among species (Tedersoo et al. 2009, Rochet et al. 2011. To assess the validity of the 97% threshold for sequences based on only ITS1 versus the full ITS region (i.e. ITS1, 5.8S, and ITS2), we examined seven known Alnus-associated Tomentella taxa (i.e. those present in Kennedy et al., 2011) and found that that threshold resulted in the same number of OTUs in both cases (data not shown). The UNITE database (Kõljalg et al., 2013) was used in both chimera checking and OTU clustering, with singleton OTUs discarded to minimize the effects of artifactual sequences (Tedersoo et al. 2010b). We assigned taxonomic data to each OTU with NCBI BLAST+ v2.2.29 (Altschul et al., 1990), using a custom fungal ITS database containing the curated UNITE SH database (v6) (http://unite.ut.ee/repository.php, Kõljalg et al., 2013) and more than 600 vouchered fungal specimens, including 46 representative sequences from Alnus forests at other HSC locations in Oregon (Kennedy & Hill 2010) and Mexico (Kennedy et al. 2011). Since sequences that had low subject length:query length matches were typically nonfungal, we further filtered out sequences with matches ≤ 90% to BLAST (i.e. at least 90% of the bases in the input sequence matches to another sequence in the database at some identity level).
Using the remaining sequence dataset, we rarefied all samples to 12946 sequences, which was the lowest number of sequences obtained across the 192 samples. Since there has recently been a question raised about the validity of rarefaction in next generation sequencing analyses  (Table S1), so present the data based on rarefied samples only. ECM OTUs within each sample were parsed out using a python script that searches for genera names from a list of 189 known ECM genera and their synonyms (Branco et al., 2013, appended from Tedersoo et al., 2010a. While this script provides a strong general filter for sorting the data by fungal lifestyle, some taxa belonging to clades that are polyphyletic for the ECM habit (e.g. Lyophyllum, Sebacinales) as well as taxa with low matches to Genbank (e.g. Uncultured Fungus) can be of questionable trophic status. For each of these groups, we carefully checked both the sequence matches and placement of our OTUs within phylogenetic trees of the clades to determine whether these taxa were properly classified at ECM. The resulting sample x OTU matrix contained 190 ECM taxa represented by at least one sequence per sample (min = 1, median = 34, mean = 1334, max = 209,187). We found that 15 of the 190 OTUs present were highly similar (>97% similar) to ECM fungi present in the dipterocarp rainforests of Malaysia, which were concurrently being studied in the Peay lab using the same next-generation sequencing approach (Fig. 1). Because these OTUs represented accidental contamination probably during library construction, they were eliminated from the final analyses. Although an additional 80 OTUs had greater than >97% similarity to taxa found in the Borneo study, because their closest BLAST match was not from Borneo, we conservatively considered these taxa as having cosmopolitan distributions and included them in the final analyses. The final OTU × sample matrix, including taxonomic matches and representative of sequences for each OTU, can be found in Table S2. Reviewing Manuscript samples). We utilized the C-score algorithm (Stone & Roberts 1990), which compares the number of checkerboard units (i.e. 1,0 x 0,1) between all pairs of species in the observed matrix (Cobserved) to that based in random permutations of the same matrix (Cexpected, i.e. the null models).

Statistical Analyses
Since randomized permutations of a matrix can be achieved in multiple ways (see Gotelli & Entsminger 2009 for details), we analyzed our datasets using both the 'fixed-fixed' and 'fixedequiprobable' options (which are recommended by the program guide and used in the previous ECM fungal co-occurrence analyses). In both options, the row (i.e. taxon) totals were fixed, so that the total abundances of each taxon in the observed and null matrices were identical. In the 'fixed-equiprobable' option, however, the column (i.e. sample) totals in the null matrices were no longer equivalent to those in the observed matrix. Instead, all samples in the null matrices had an equal probability of being colonized by any of the taxa in the observed matrix, which effectively eliminates differences in taxon richness among samples.
Of the ECM taxa present in the final root tip and bag datasets, over 90% (167/183) belonged to species never previously encountered with Alnus (Table S2, AlnusMatch = No).
Unlike other ECM host systems with large geographic ranges, the ECM fungal community associated with Alnus hosts is remarkably well characterized at local (Tedersoo et al., 2009, Kennedy et al., 2010, Walker et al., 2014), regional (Kennedy et al., 2011, Roy et al., 2013, and global scales (Polme et al., 2013. As such, it is highly likely the majority of the novel OTUs encountered were not part of the active ECM community in our plots, but rather present simply either as spores or additional lab contaminants. To account for this issue, we divided our checkerboard analyses into five different input matrices for the bag dataset (Plots 2 and 8). The first matrix included all 183 ECM fungal taxa (referred to as "All"). The second matrix included the 16 taxa that had >97% similarity matches to ECM samples from Alnus forests (referred to as 4 (referred to as AlnusRootOnly). To assess the robustness of the results generated using the larger Alnus matrix, the fourth matrix excluded the three most frequent and abundant species (Tomentella3, Alnicola1, Tomentella2) (referred to as AlnusMinusTop3). Finally, the fifth matrix included just the 10 taxa in the genus Tomentella (from the larger Alnus matrix) to look for evidence of species interactions among this subset of closely related taxa (referred to as AlnusTomentellaOnly). For all of the aforementioned C-score analyses, taxa present in less than 5 bag samples were removed, as low frequency taxa are generally considered non-informative (Koide et al., 2005). The observed input matrices were compared to 5000 null matrices.
Significant differences between the observed matrix C-score and that of the null matrices were determined along with standardized effect sizes (SES). Observed C-scores significantly higher than those generated from the null matrices are consistent with a community being structured by competitive interactions, whereas Cobserved significantly lower than the Cexpected is consistent with positive interactions.
To further assess the degree of association between known Alnus ECM fungal taxa, we also used an abundance-based approach (as opposed to the co-occurrence analyses, which are based on binary presence/absence data). Specifically, we calculated the pair-wise Spearman rank correlation coefficients among all pairs of the 16 Alnus-associated taxa using the cor function in R (R Development Core Team 2013). Coefficients >0.30 were tested for significance with the correlogram tests, we used the n.class=0 option, which uses Sturge's equation to determine the appropriate number of distance classes.

Results
We found 183 total ECM fungal taxa in the study (Table S2); 16 of which matched closely to known Alnus-associated ECM fungi. In the mesh bags, Alnus-associated ECM fungal taxa represented six of the ten most abundant OTUs present, including the dominant ECM fungal taxon, Tomentella3, which was present in all the bag samples in both plots and had sequence abundances nearly ten-fold higher than any other taxon (Fig. 2a,b). Two other Alnus-associated taxa, Alnicola1 and Tomentella2, were also present in all samples, whereas the remaining Alnusassociated ECM fungal taxa had frequencies varying from 2-96% (Plot 2 mean = 25%, Plot 8 mean = 31%) and lower sequence abundances. Eight of the 16 Alnus-associated ECM fungal taxa were present on both roots and in the bags, with abundances that were very similar (Fig. 1a). Of the eight ECM fungal taxa found on root tips, all were previously encountered on A. rubra root tips at other sites in Oregon, while the eight taxa found exclusively in bags had not been previously documented (Kennedy & Hill 2010).
ECM fungal taxon co-occurrence patterns were largely consistent between plots, but different between null models. Of the ten tests (i.e. 5 matrix types x 2 plots) using the 'fixedfixed' permutation option, nine indicated that the observed ECM fungal community did not differ significantly from random assembly (Table 1). In one case, Plot 2 All, the observed ECM fungal community had significantly more co-occurrence than expected by chance. In contrast, in the ten tests using the 'fixed-equiprobable' permutation option, three indicated that the observed ECM fungal community did not differ significantly from random assembly, while seven found that the observed ECM fungal community had significantly more co-occurrence than expected by chance.
Results remained the same for Alnus ECM fungal communities whether the top three taxa were removed or not. The Alnus and AlnusRootOnly analyses did differ under the 'fixed-equiprobable' option, with the former showing greater than expected co-occurrence and the latter having a pattern no different than one based on random assembly. Additionally, in the AlnusTomentellaOnly analysis, the ECM fungal community showed greater than expected cooccurrence in Plot 2 but not in Plot 8. In all of these cases, significant antagonistic patterns were not observed.
Spearman rank analyses revealed that pair-wise sequence abundances of some of the 16 Alnus ECM fungal taxa were significantly positively correlated ( Table 2) there was a single significant positive correlation between samples located 1-2 m apart (Fig. S2,   S3).

Discussion
We found that the ECM fungal communities in A. rubra forests displayed a different pattern of taxon co-occurrence compared to those seen for other ECM fungi. Unlike the consistent previous findings of less co-occurrence among species than expected by chance (Koide  et al., 2005;Pickles et al., 2012;Wadet et al., 2012), we observed no evidence of spatial patterns consistent with interspecific competition in Alnus-associated ECM fungal communities. In contrast, we consistently found co-occurrence patterns that were either no different from random assembly or consistent with positive interactions. Although we did not measure soil nitrate and acidity conditions in this study (see Martin et al., (2003) and Walker et al., (2014) for values from comparable age A. rubra forests at other sites in Oregon), Alnus soils are consistently characterized by abiotic conditions are generally considered stressful to ECM fungi. The results we obtained are thus consistent with the 'stress gradient hypothesis', which posits that species interactions shift from negative to positive as environmental conditions become harsher (Bertness & Callaway 1994). Although we emphasize that the patterns we found in this study are based solely on correlative inference, there is some experimental evidence that may support the stress gradient hypothesis for ECM fungal community dynamics. Koide et al., (2005) found a shift from significant negative co-occurrence patterns in their control plots to non-significant co-occurrence patterns in plots where either tannins or nitrogen were added experimentally. While they did not explicitly analyze these manipulations in terms of stress, both increased tannin and nitrogen levels have been shown to inhibit the growth of multiple ECM taxa (Koide et al., 1998;Cox et al., 2010). The direction of the response in the Koide et al., (2005) study is consistent with greater abiotic stress resulting in a decrease in antagonistic ECM fungal interactions. At the same time, it is plausible that resource limitation was eliminated with the addition of nitrogen, which could have allowed for greater spatial co-existence among ECM fungi. Since the Alnus system has naturally higher nitrogen availability than most ECM forests due to the co-presence of nitrogenfixing Frankia bacteria, it is also possible that greater resource abundance could drive the cooccurrence patterns we observed. Given the fact that the pattern could be explained by either increasing stress or resource availability, additional experimental tests are needed to distinguish among these explanations. One promising approach would be to examine the taxon co-occurrence patterns in younger and older Alnus forests, since soil nitrate and acidity concentrations increase in these forests over time (Daniere et al. 1986, Martin et al., 2003. If the stress gradient hypothesis were the most plausible explanation, then we would expect to see competitive and facilitative interactions to be dominant, respectively.
The presence of co-occurrence patterns consistent with significant negative species interactions was also missing in our analysis of more closely related ECM fungal taxa. For the ten Alnus-associated members of the genus Tomentella, co-occurrence patterns either did not differ Some positive spatial associations have been observed in other studies of ECM fungal communities (Agerer et al., 2002;Koide et al., 2005;Pickles et al., 2012), and have been suggested to be due to complementary resource acquisition abilities of among individual taxa (Jones et al., 2010). We speculate that in Alnus forests positive associations among ECM fungi could also reflect possible amelioration of local abiotic conditions. Huggins et al., (in press) found that Alnus-associated ECM fungi could more effectively buffer changes in local pH environments than non-Alnus ECM fungi, which may be key to persistence in the high acidity soils present in Alnus forests. While the exact buffering mechanism is not yet known, if it involves the release of molecules into the external environment, growing directly adjacent to another ECM fungus may result in greater buffering of local pH conditions than when growing in isolation. We believe it is important to note, however, that the patterning of positive associations were patchy and not consistent between plots, so it is hard to determine if local pH buffering is actually significant without local measurements of pH for each sample. Furthermore, sequence abundance of individual taxa has been shown not to correlate linearly with initial fungal tissue or DNA abundance in other studies using NGS techniques (Amend et al., 2010, Nguyen et al., 2014, so caution must be applied in using sequence abundance as an accurate ecological proxy. Like the co-occurrence and correlation-based patterns, we found that spatial autocorrelation patterns observed in Alnus ECM fungal communities were also anomalous relative to other studies. The specific distance of spatial autocorrelation appears to vary among systems, but there is typically strong spatial autocorrelation among community samples located less than 5 m apart (e.g. Lilleskov et al. 2004, Bahram et al. 2013. While the spatial extent of our study was very limited (the most distant samples within plots were only ~12 m apart), the found functionally identical results to those when those taxa were included (Table 2). A second difference between this and related studies was the sampling of ECM fungal hyphal communities in mesh bags. Previous studies assessing co-occurrence patterns have largely focused on ECM root tips, but Koide et al., (2005) found very similar taxon co-occurrence patterns for root-tip and soil-based analyses of ECM fungal communities in the same Pinus resinosa forest. Based on that result, and the fact that the sequence abundances of all the ECM fungi present on A. rubra root tips and the mesh bags showed highly similar patterns (Table S2), we do not believe assessing ECM hyphal communities was the source of our incongruous results either (however, in hindsight, a better experimental design would have been to sample the mesh bags and the ECM root tips directly around them within each plot). A third difference is the restricted taxonomic richness of Alnus ECM fungal communities. This explanation, however, also seems unwarranted, as Pickles et al., (2012) showed highly significant negative co-occurrence patterns in matrices of equivalent sizes. It is also possible that variation in soil nutrient availability could drive Alnus ECM fungal community structure and, because it was relatively homogenous in our small-sized plots, the resulting taxon distribution patterns were largely random. While we reiterate that we did not directly measure soil nutrient availability in this study, other studies of Alnus ECM fungi have shown some significant correlations between community structure and soil organic matter and nutrients such as K and Ca (Becerra et al., 2005, Tedersoo et al., 2009, Roy et al., 2013, Polme et al., 2013see Richard (1968) for a possible mechanism). In those studies, however, the percent of variance explained by soil nutrients was generally low, so we believe it is unlikely that variation in resource availability was the primary determinant of the distribution patterns observed. We recognize that additional differences likely exist, but feel confident that the cooccurrence results we observed are ecologically accurate and not generated by methodological or sampling artifact. NGS techniques clearly represent a powerful and efficient way to assess the richness and dynamics of fungal communities (Smith & Peay 2014), but we found that additional data quality control analyses beyond the standard sequence quality thresholds and chimera checking were needed to properly characterize ECM fungal community composition. Specifically, we found that a relatively high number of ECM fungal taxa present that appeared to be the result of PCR contamination. The PCR reactions of our extraction and PCR controls produced no bands indicating positive product, but the sensitivity of NGS techniques and the Illumina platform in particular makes the amplification of single DNA molecules highly probable (Tedersoo et al., 2010b;Peay et al., 2013). Fortunately, the atypical and well-described nature of Alnus ECM fungal communities made it relatively easy to identify the most obvious non-Alnus associated taxa and remove them prior to the final analyses. For taxa that belonged to ECM fungal lineages known to associate with Alnus hosts but which had not been previously documented, it was more difficult to determine their status (i.e. whether they represented PCR contaminants, were present in A. rubra soils as spores, or actually colonizing A. rubra root tips). In particular, the status of Thelephoraceae1, which had the third highest sequence abundance in the full dataset, was interesting because the closest BLAST match to Thelephoraceae1 was an ECM fungal root tip sample from Betula occidentalis in British Columbia, Canada. Bogar & Kennedy (2013) found that ECM fungal communities present on Alnus and Betula hosts can overlap, so it is possible this taxon was overlooked in previous surveys of Alnus ECM fungal communities that used less sensitive methods. However, the absence of this taxon from any the root tip samples in Plot 4 inclusion of spurious taxa is sufficiently high that we strongly recommend the sequencing of negative extraction and PCR controls to help try to account for any lab-based contamination (Nguyen et al., in press).
Taken together, our results suggest that while many ECM fungal communities appear to be strongly affected by competitive interactions, those present in Alnus forests are not. Although the reasons for this difference are not fully resolved in this study, the possibility of greater abiotic stress changing the way in which species interact in Alnus forests is likely an important factor.
The application of ecological theories such as the stress gradient hypothesis to better understand the factors driving ECM fungal community structure has grown rapidly in recent years (Peay et al., 2008;Koide et al., 2014) and new technologies such as next generation sequencing continue to make the study of ECM fungi increasingly tractable for ecologists. While we welcome this synergy, we stress the importance of a solid foundation in fungal biology as well as critical The top 20 ECM fungal taxa are color coded by whether there are known to be associated with Alnus hosts (black), of unknown host origin (grey), or laboratory contaminants (white).