Fungi in soil and understory have coupled distribution patterns

Study site

Field sampling consisted of two separate forays in August 2017 (Aneityum) and in December 2017 (Tanna) in the province of Tafea, Vanuatu as previously described (Boraks et al., 2021). The two islands are separated by 86 km of open ocean and share a similar tropical climate. Aneityum receives 2,322 mm of rain annually and has an average annual temperature of 19.7 °C, although there is large seasonal and interannual variation (Vanuatu Meteorology and Geo-hazard Department, Australian Bureau of Meteorology, CSIRO, 2015). Our sampling locations are classified as low to mid-elevation rain forest (Mueller-Dombois & Fosberg, 2013) dominated by Syzygium, Dysoxylum, and Hedycarya species (Boraks et al., 2021).

Six forest transects (10 m by 40 m) were sampled, three on each island (Fig. S1). Within each transect we sampled 36 locations in grid formation. At each sampling location, we collected the fungal community from both the soil and the leaf surface (phylloplane). Fungal communities were harvested using flocked sterile swabs and preserved in CTAB buffer. Soil samples were collected by brushing away loose plant litter and inserting swabs into the organic layer of topsoil. Phylloplane samples were collected by swabbing leaf surfaces. Leaf sampling height varied by availability, but generally were located between 1 and 3 m above the forest floor. The total area sampled for each phylloplane sample was constrained to a surface area of 400 cm², or roughly the surface area of two hands. Plant taxonomy influences the phyllopsphere mycobiome to a greater (Kembel & Mueller, 2014), or lesser (O’Rorke et al., 2015), degree and we attempted to reduce this variance by swabbing leaves from multiple plant species using the same swab. We aimed to include at least three different plant species per sampling site.

Soil and phylloplane samples were taken in pairs from the same geographic location, such that the phylloplane-sampling grid and the soil-sampling grid were coupled. This paired-habitat and spatially explicit sampling scheme totaled 432 fungal community samples (36 sampling sites × 2 habitats × 6 transects) (Boraks et al., 2021).

DNA extraction, PCR and sequencing

Sequencing library preparation was previously described in Boraks et al. (2021). DNA was extracted from the swabs using the Qiagen DNeasy PowerSoil DNA Isolation Kit (Qiagen, Maryland). We slightly modified the extraction protocol by extracting directly from swabs. Amplicon libraries were prepared in a single PCR reaction using Illumina-barcoded fungal-specific primers as in Amend et al. (2019). The primers ITS1F and ITS2 were used to target the hypervariable nuclear ribosomal ITS1 region which is flanked by the 18S and 5.8S rDNA regions (Schoch et al., 2012). Sequence library was purified and normalized (Just-a-plate; Charm Biotech), quantified by qPCR and sequenced on Illumina MiSeq platform using V3 chemistry (Illumina Inc., San Diego, CA) at the Institute for Integrative Genome Biology (University of California–Riverside). Additional details of sequence library prep can be found in Boraks et al. (2021). Sequences are deposited in the NCBI Sequence Read Archive (BioProject: PRJNA634909).

ITS1 sequences were extracted from the flanking ribosomal subunit genes using ITSxpress (Rivers et al., 2018), then filtered by quality scores using the FASTX-Toolkit (Gordon & Hannon, 2010). Reverse sequence reads were discarded because of lower quality reads. Chimeras were detected and removed using vsearch (Rognes et al., 2016). Sequences were clustered at 97% identity and fungal taxonomy was assigned using the Python package constax (Gdanetz et al., 2017) and the Unite 8.0 database. Putative contaminants were identified based on their prevalence in extraction and PCR negative controls and removed using the R package decontam (Davis et al., 2017). OTUs that could not be assigned to a fungal phylum were discarded. The sequence library was normalized by variance stabilizing transformation (VST) (McMurdie & Holmes, 2014) in DESeq2 (Love, Huber & Anders, 2014) within R (R Core Team, 2019). Singletons were removed for all analyses except those that measured community dissimilarity. Sequence read abundance data, taxonomic assignments, sample metadata, and ancillary collection data were compiled with the R package phyloseq (McMurdie & Holmes, 2013). Using sequence read abundance as a metric for species abundance is contentious due to the compositional nature of high-throughput sequencing datasets (Gloor et al., 2017). Statistical analyses for this study rely on both presence-absence and sequence abundance data. Caution should be exercised when interpreting analyses that use relative sequence read abundance as metric for OTU abundance.

Statistical analyses

Area of occupancy

To compare range sizes between habitats (soil or phylloplane) and lifestyle (specialist or generalist), we first categorized OTUs by habitat and lifestyle. Habitat specialization was assigned if an OTU was found more than once and exclusively in a single habitat type, whereas habitat generalists occurred more than once in each habitat from the same transect (Fig. S2). Occupancy was calculated by the number of samples in which an OTU occurred within a transect. This analysis was constrained within transects so the maximum occupancy was limited to the number of sampling points within a transect (36 locations, or roughly 400 m²). A Wilcoxon rank sum test, corrected for multiple hypothesis testing (Benjamini & Hochberg, 1995), was used to compare median range size.

Range size per substrate and specialization

To measure the range size of individual OTUs, we identified those that were spatially autocorrelated using a Mantel test with Bray–Curtis distances based on OTU relative sequence-read abundance and then built empirical semi-variogram models for each spatially autocorrelated OTU. Semi-variogram models are useful for analyzing the range of variation between geospatial phenomena. In this study we used semi-variograms to measure the geospatial decay of OTU relative sequence read abundance. Semi-variogram models were developed for each spatially autocorrelated OTU and geographic-range sizes were estimated from the range of each variogram model. The nugget (zero) and model type (exponential) were fixed for all variograms to increase comparability of range size between models. Range sizes were aggregated by habitat and specialization and then compared by pairwise Wilcoxon tests corrected for multiple hypothesis testing. Kriged surface maps were paired for the generalist OTU that occurred in both habitats of the same transect.

Occupancy-abundance relationships

To compare OTU rarity, we generated occupancy-abundance (OA) regressions for all OTUs of each habitat and then compared linear regressions (loess fit) for specialists and generalists of the soil and phylloplane. OTUs were first categorized by habitat and specialization, and then the number of occurrences per transect was tallied for each OTU. This type of OA regressions is known as interspecific occupancy-abundance relationships because the regressions compare multiple OTUs. The maximum occupancy of any single OTU is 36, because there were 36 sampling locations per transect. Abundance was measured using normalized sequence reads.

Are fungal communities different between habitats?

When all transects were combined, 4949 OTUs were identified in the soil and 5362 in the phylloplane. Ten OTUs were previously identified as suspect contaminats and were removed from the dataset. Transects contained an average of 4415 OTUs (standard error of mean (SE) ± 340; n = 6). On average, 55% (SE ± 7%) of the OTUs that were present in a transect could be found in the soil habitat. This was similar for phylloplane, in which 65% (SE ± 5%) of the OTUs that were present in a transect could be found in the phylloplane. Permanova indicated significant differences (p <  0.001) in community composition between soil and phylloplane (habitat R² = 0.11, p < 0.001) and island (R² = 0.03, p < 0.001; Fig. 1A).

Inter-habitat spatial autocorrelation

Composition of fungal communities of the soil and the phylloplane differ, yet 22% of OTUs co-occur in both habitats (generalists). Co-occurrence between habitats was spatially autocorrelated at a large scale (110 km), phylloplane communities were more similar to proximate soil communities than distant soil communities (Table 1). At a smaller scale (<40 m) the spatial autocorrelation of fungal communities from different habitats is less supported. Inter-habitat spatial autocorrelation was statistically significant for three of the six transects (Table 1).

Area of occupancy

An OTU was categorized as a habitat specialist if it was found more than once and exclusively among a single habitat type. We identified 2,187 soil specialist OTUs and 2,212 phylloplane specialist OTUs (Fig. 1B). We defined OTU generalists as those occurring more than once in each habitat. A total of 1307 unique OTUs were classified as generalists. Thirty-five percent of the generalist OTUs occurred in more than one transect. Generalist fungi had higher occupancy rates than specialist fungi (Fig. 2). Median occupancy rates were greater for generalist fungi (phylloplane generalist mdn = 6, iqr = 12; soil generalist mdn = 4, iqr = 6) than for specialist fungi (phylloplane specialist mdn = 3, iqr = 3; soil specialist mdn = 3, iqr = 4) (Fig. 2).

Semi-variogram modeling

Spatial autocorrelation was detected among 7.7% of the phylloplane specialists and 10.4% of the soil specialists. Spatial autocorrelation was detected among 10.5% of the phylloplane generalists and 18.8% of soil generalists. Range sizes were modeled using experimental variograms that were based on the number of sequence reads as a function of spatial distance. A total of 1,060 variogram models were generated. Variogram models were manually inspected, and we culled models (515/1060) that reported unrealistically small ranges (<1 m) for our given sampling scheme. In addition, 6 OTUs were identified as outliers and, despite robust variogram models, had exceptionally large ranges and were removed from downstream analysis (Table S1). The six outliers with large range sizes (mean = 737 m) included five OTUs within the Ascomycota and one OTU within Glomeromycotina. The removal of unreliable models and outliers resulted in 539 variograms remaining. Range size of specialist soil fungi (2.78 m; iqr 1.77) were significantly greater than phylloplane generalist fungi (2.11 m; iqr 1.53) and phylloplane specialist fungi (2.02 m; iqr 1.82) (p < 0.05). Soil generalist fungi (2.54 m; iqr 1.66) did not significantly differ in range size from soil specialists or phylloplane generalists (Fig. 3). Inclusion of the culled outliers OTU did not affect the significance of range-size comparisons. Spatially autocorrelated OTU sequence read abundance generally become insignificant around at distances greater than 5 m, although several ranges were estimated to be hundreds or thousands of meters.

Distribution maps for generalist OTU

Kriged surface maps were produced for 11 OTUs that were identified as spatially-autocorrelated generalist fungi that occurred in the same transect. The OTUs were identified as OTU97_1160_Dothideomycetes sp.; OTU97_1992_Helotiales sp.; OTU97_2240_Agaricomycetes sp.; OTU97_228_Psathyrella sp.; OTU97_257_Polyporale sp.; OTU97_2608_Xylariales sp.; OTU97_2833_Agaricales sp.; OTU97_3832_Psathyrellaceae sp.; OTU97_3832_Psathyrellaceae sp.; OTU97_5320_Xylariales sp.; OTU97_6890_Ascomycota sp. In general, the spatial distribution of an OTU was synchronous between habitats (Fig. 4). Semi-variogram models and kriging standard errors are reported in the supplemental (Fig. S3)

Occupancy-abundance relationships

Phylloplane and soil occupancy-abundance relationships followed classical biogeographical expectations; the number of samples in which an OTU occurred, was positively correlated with the mean abundance of that OTU (Fig. 5). Among the phylloplane, specialists and generalists had similar occupancy-abundance regression slopes at low occupancy levels, but the regressions diverged at high occupancy levels (Fig. 5). Common phylloplane generalists were more abundant than common phylloplane specialists. Among the soil, low-occupancy specialist fungi were more abundant than low occupancy generalists, but specialist and generalist abundance converge at higher occupancy rates (Fig. 5).

Co-occurrence among habitat types

Here, we report a first attempt to couple distribution maps of fungal OTUs in the soil and understory. These map stacks provide a rich visualization of where generalist fungi persist simultaneously in either habitats across the forest landscape (Fig. 4). In addition to the visual evidence, inter-habitat spatial autocorrelation confirmed that fungal communities of the phylloplane are more likely to be similar to nearby soil communities than to distant soil communities (Table 1). The inter-habitat spatial autocorrelation observed here was especially pronounced at large spatial scales (<110 km) but became less evident at smaller spatial scales (<40 m). Stated differently, this result indicates the possibility of predicting the geo-location of a leaf surface based on the fungal community of soil samples (or vice versa), and that the accuracy of this prediction becomes less reliable at small spatial scales (<40 m).

Migration between habitats is important for community turnover of ephemeral habitats such as pools of water, phytotelma, and plant leaves. Leaf age is directly correlated with mycobiome composition. Young, newly emergent plant leaves are sparsely colonized by fungi (Wildman & Parkinson, 1979; Suryanarayanan & Thennarasan, 2004; Oono et al., 2015). The initial colonization is a stochastic processes via horizontally transmitted fungi from unknown environmental sources (Unterseher et al., 2018). As the leaf ages, fungal infection increases and beta-diversity decreases (Oono et al., 2015). The decrease in beta-diversity may be partially driven by host identity (Unterseher et al., 2013; Kembel & Mueller, 2014; O’Rorke et al., 2015). We attempted to compensate for the effect of host identity by sampling the leaves of multiple plant hosts from the same sampling location. After a leaf has senesced and fallen from the plant, it will carry a rich microbial community to the forest floor. The phylloplane and soil are physically disjunct habitats but spatially autocorrelated overlap in their mycobiomes (Table 1) is suggestive of migration between them.

It is well established that microbes from the forest floor migrate to the phylloplane and vice versa. Indeed, this is the main tenant of the latent saprotroph hypothesis which proposes that a saprotroph is able to colonize the living tissue of its host to gain competitive advantage when the host tissue begins to senesce (Porras-Alfaro & Bayman, 2011; Peršoh, 2013). Phytopathology literature contains examples of foliar pathogens that overwinter in ground litter, only to migrate back to the phylloplane when conditions are optimal. For example, the tree pathogen Venturia inaequlis is a hemibiotroph and must spend part of its lifecycle as a parasite, and part of its lifecycle on the dead tissue of the same host as a saprotroph (Agrios, 2005). In line with our results (59%, n = 3,519), a previous report examining the overlap of phyllosphere and leaf litter found 64% of leaf epiphytes could also be isolated from leaf litter (Osono, 2006). An important distinction in our study is that we sampled from the organic soil layer rather than leaf litter. Evidence from bacterial (Bai et al., 2015; Zarraonaindia et al., 2015; Wagner et al., 2016) and fungal systems (Amend et al., 2019; Bernard et al., 2021) have shown that soil may act as a reservoir of microbes that migrate to the phylloplane, an observation that is supported by our results (Table 1).

Habitat breadth and range size

Range size and niche breadth are positively correlated among plants and animals (Slatyer, Hirst & Sexton, 2013; Bernard et al., 2021). Here we show that for fungi, niche breadth may or may not correlate with range size, and that this equivocal result depends on how range size is calculated (Figs. 2 and 3). Range size and habitat specificity are two of the major axes of rarity (Rabinowitz, 1981) and by including sequence read abundance we explore commonness and rarity of fungi from the phylloplane and soil.

There is no standard methodology for measuring species’ geographic range. The large suite of techniques available can capture fundamentally distinct features of species distribution (Gaston, 1996). We used two different approaches to measure fungal geographic range size. First, we considered area of occupancy per transect (Fig. 2). This is a tally of the number of sampling-sites at which an OTU is recorded within a transect. Calculating range size by areas occupied is advantageous because of its simplicity, it requires only presence-absence data and it is inclusive of rare species (Gaston, 2008), and was successfully used in a study of bacterial range size (Choudoir et al., 2018). The drawback to this approach is that it lacks distribution nuance; it is unable to differentiate between a sparse and clustered distribution of the same occupancy level. The second method we used to measure range size was variography which provided a measurement most similar to Gaston’s (1996) geometric circle (Fig. 3). In this study the radius a range was calculated from the range of an empirical semi-variogram. This technique was previously used in a study of ectomycorrhizal fungi (Pickles et al., 2010). Another way to think of this measurement is to consider the linear distance, from the geometric center of a range, at which variation in OTU sequence abundance is no longer spatially autocorrelated. The advantage of the variogram approach is that it provides a linear distance measurement which can be used to make visually illustrative maps (Fig. 4). When interpreting kriged variogram maps (Fig. 4) generated using high throughput sequencing data, it is important to consider their limitations. Variogram modeling is only possible on a small fraction of our metabarcoding dataset (Fig. S2) because of the large amounts of spatially autocorrelated data required and most OTU datasets are dominated by low-abundance taxa. Relative sequence read abundance is compositional and therefore comparisons between samples may be misleading (Gloor et al., 2017; Morton et al., 2019). Additionally, amplicon sequencing is indiscriminate and we are unsure whether we sequenced spores, hyphae, or exogenous DNA (Carini et al., 2016).

Brown (1984) affirmed that species which have a broad environmental range and are able to use a wide range of resources will be both locally abundant and widespread. For fungi, adjacency of steppingstones within a habitat matters little if you can just as easily step between habitats. In an analysis of soil bacteria, Luo et al. (2019) found that generalist bacteria occurred at greater abundances than specialist bacteria. Results from our area of occupancy analysis support theoretical expectations insofar as we found that generalist fungi occur more frequently than specialists (Fig. 2). Interestingly, greater transect occupancy levels did not translate into larger modeled range sizes. Variogram estimations indicated no significant differences in range sizes between habitat specialists and generalists (Fig. 3). Lack of difference in the variogram range sizes might be attributable to the insensitivity of variograms to low occupancy OTU. Differences in range size results from each technique highlights the complicated nature of spatial distribution and habitat specialization.

A positive interspecific abundance–occupancy relationship is one of the most widely observed patterns in macroecology. We tested for variation in abundance-occupancy regression slopes between specialists and generalists from two separate habitats. OTUs of this study exhibited a positive occupancy-abundance relationship which describes the greater-than-proportional increase in abundance with respect an increase in occupancy (Fig. 5). When transect occupancy is low, soil specialists are more abundant than soil generalists (Fig. 5). We believe this may be result of the comparative stability of soil. Without significant perturbation, soil microbial communities show limited turnover (Kivlin et al., 2014; Žifčáková et al., 2016) potentially due to the longevity of spore banks (Bruns et al., 2009; Nguyen, Hynson & Bruns, 2012). In less disturbed habitats, a strong competitor can persist in one location, cumulatively adding to the volume of spores within the spore bank of that location. This process might explain the pattern of locally abundant soil specialists towards the left side of Fig. 5B. The result of this study are similar to a previous report that belowground fungi had high local abundance but small species ranges, while aboveground fungi had lower local abundance but larger spatial ranges (Kivlin et al., 2014). Generalists occur more commonly and abundantly among the phylloplane, whereas in the soil specialist fungi were found at high abundance but low occupancy (Fig. 5). By contrasting the physical and temporal characteristics of these two habitats we were able to infer potential strategies employed by fungi to persist on the landscape.