Airborne DNA reveals predictable spatial and seasonal dynamics of fungi.

Fungi are among the most diverse and ecologically important kingdoms in life. However, the distributional ranges of fungi remain largely unknown as do the ecological mechanisms that shape their distributions1,2. To provide an integrated view of the spatial and seasonal dynamics of fungi, we implemented a globally distributed standardized aerial sampling of fungal spores3. The vast majority of operational taxonomic units were detected within only one climatic zone, and the spatiotemporal patterns of species richness and community composition were mostly explained by annual mean air temperature. Tropical regions hosted the highest fungal diversity except for lichenized, ericoid mycorrhizal and ectomycorrhizal fungi, which reached their peak diversity in temperate regions. The sensitivity in climatic responses was associated with phylogenetic relatedness, suggesting that large-scale distributions of some fungal groups are partially constrained by their ancestral niche. There was a strong phylogenetic signal in seasonal sensitivity, suggesting that some groups of fungi have retained their ancestral trait of sporulating for only a short period. Overall, our results show that the hyperdiverse kingdom of fungi follows globally highly predictable spatial and temporal dynamics, with seasonality in both species richness and community composition increasing with latitude. Our study reports patterns resembling those described for other major groups of organisms, thus making a major contribution to the long-standing debate on whether organisms with a microbial lifestyle follow the global biodiversity paradigms known for macroorganisms4,5.

Experimental evidence that root-associated fungi improve plant growth at high altitude.

Unravelling how species communities change along environmental gradients requires a dual understanding: the direct responses of the species to their abiotic surroundings and the indirect variation of these responses through biotic interactions. Here, we focus on the interactive relationships between plants and their symbiotic root-associated fungi (RAF) along stressful abiotic gradients. We investigate whether variations in RAF community composition along altitudinal gradients influence plant growth at high altitudes, where both plants and fungi face harsher abiotic conditions. We established a translocation experiment between pairs of Bistorta vivipara populations across altitudinal gradients. To separate the impact of shifting fungal communities from the overall influence of changing abiotic conditions, we used a root barrier to prevent new colonization by RAF following translocation. To characterize the RAF communities, we applied DNA barcoding to the root samples. Through the utilization of joint species distribution modelling, we assessed the relationship between changes in plant functional traits resulting from experimental treatments and the corresponding changes in the RAF communities. Our findings indicate that RAF communities influence plant responses to stressful abiotic conditions. Plants translocated from low to high altitudes grew more when they were able to associate with the resident high-altitude RAF compared to those plants that were not allowed to associate with the resident RAF. We conclude that interactions with RAF impact how plants respond to stressful abiotic conditions. Our results provide experimental support that interactions with RAF improve plant stress tolerance to altitudinal stressors such as colder temperatures and less nutrient availability.

Symbiotic status alters fungal eco-evolutionary offspring trajectories.

Despite host-fungal symbiotic interactions being ubiquitous in all ecosystems, understanding how symbiosis has shaped the ecology and evolution of fungal spores that are involved in dispersal and colonization of their hosts has been ignored in life-history studies. We assembled a spore morphology database covering over 26,000 species of free-living to symbiotic fungi of plants, insects and humans and found more than eight orders of variation in spore size. Evolutionary transitions in symbiotic status correlated with shifts in spore size, but the strength of this effect varied widely among phyla. Symbiotic status explained more variation than climatic variables in the current distribution of spore sizes of plant-associated fungi at a global scale while the dispersal potential of their spores is more restricted compared to free-living fungi. Our work advances life-history theory by highlighting how the interaction between symbiosis and offspring morphology shapes the reproductive and dispersal strategies among living forms.

Elevation and plant species identity jointly shape a diverse arbuscular mycorrhizal fungal community in the High Arctic.

Knowledge about the distribution and local diversity patterns of arbuscular mycorrhizal (AM) fungi are limited for extreme environments such as the Arctic, where most studies have focused on spore morphology or root colonization. We here studied the joint effects of plant species identity and elevation on AM fungal distribution and diversity. We sampled roots of 19 plant species in 18 locations in Northeast Greenland, using next generation sequencing to identify AM fungi. We studied the joint effect of plant species, elevation and selected abiotic conditions on AM fungal presence, richness and composition. We identified 29 AM fungal virtual taxa (VT), of which six represent putatively new VT. Arbuscular mycorrhizal fungal presence increased with elevation, and as vegetation cover and the active soil layer decreased. Arbuscular mycorrhizal fungal composition was shaped jointly by elevation and plant species identity. We demonstrate that the Arctic harbours a relatively species-rich and nonrandomly distributed diversity of AM fungi. Given the high diversity and general lack of knowledge exposed herein, we encourage further research into the diversity, drivers and functional role of AM fungi in the Arctic. Such insight is urgently needed for an area with some of the globally highest rates of climate change.

Distance decay 2.0 - A global synthesis of taxonomic and functional turnover in ecological communities.

Aim: Understanding the variation in community composition and species abundances (i.e., ß-diversity) is at the heart of community ecology. A common approach to examine ß-diversity is to evaluate directional variation in community composition by measuring the decay in the similarity among pairs of communities along spatial or environmental distance. We provide the first global synthesis of taxonomic and functional distance decay along spatial and environmental distance by analysing 148 datasets comprising different types of organisms and environments. Location: Global. Time period: 1990 to present. Major taxa studied: From diatoms to mammals. Method: We measured the strength of the decay using ranked Mantel tests (Mantel r) and the rate of distance decay as the slope of an exponential fit using generalized linear models. We used null models to test whether functional similarity decays faster or slower than expected given the taxonomic decay along the spatial and environmental distance. We also unveiled the factors driving the rate of decay across the datasets, including latitude, spatial extent, realm and organismal features. Results: Taxonomic distance decay was stronger than functional distance decay along both spatial and environmental distance. Functional distance decay was random given the taxonomic distance decay. The rate of taxonomic and functional spatial distance decay was fastest in the datasets from mid-latitudes. Overall, datasets covering larger spatial extents showed a lower rate of decay along spatial distance but a higher rate of decay along environmental distance. Marine ecosystems had the slowest rate of decay along environmental distances. Main conclusions: In general, taxonomic distance decay is a useful tool for biogeographical research because it reflects dispersal-related factors in addition to species responses to climatic and environmental variables. Moreover, functional distance decay might be a cost-effective option for investigating community changes in heterogeneous environments.

Microclimate structures communities, predation and herbivory in the High Arctic.

In a warming world, changes in climate may result in species-level responses as well as changes in community structure through knock-on effects on ecological interactions such as predation and herbivory. Yet, the links between these responses at different levels are still inadequately understood. Assessing how microclimatic conditions affect each of them at local scales provides information essential for understanding the consequences of macroclimatic changes projected in the future. Focusing on the rapidly changing High Arctic, we examine how a community based on a common resource species (avens, Dryas spp.), a specialist insect herbivore (Sympistis zetterstedtii) and natural enemies of lepidopteran herbivores (parasitoids) varies along a multidimensional microclimatic gradient. We ask (a) how parasitoid community composition varies with local abiotic conditions, (b) how the community-level response of parasitoids is linked to species-specific traits (koino- or idiobiont life cycle strategy and phenology) and (c) whether the effects of varying abiotic conditions extend to interaction outcomes (parasitism rates on the focal herbivore and realized herbivory rates). We recorded the local communities of parasitoids, herbivory rates on Dryas flowers and parasitism rates in Sympistis larvae at 20 sites along a mountain slope. For linking community-level responses to microclimatic conditions with parasitoid traits, we used joint species distribution modelling. We then assessed whether the same abiotic variables also affect parasitism and herbivory rates, by applying generalized linear and additive mixed models. We find that parasitism strategy and phenology explain local variation in parasitoid community structure. Parasitoids with a koinobiont strategy preferred high-elevation sites with higher summer temperatures or sites with earlier snowmelt and lower humidity. Species of earlier phenology occurred with higher incidence at sites with cooler summer temperatures or later snowmelt. Microclimatic effects also extend to parasitism and herbivory, with an increase in the parasitism rates of the main herbivore S. zetterstedtii with higher temperature and lower humidity, and a matching increase in herbivory rates. Our results show that microclimatic variation is a strong driver of local community structure, species interactions and interaction outcomes in Arctic ecosystems. In view of ongoing climate change, these results predict that macroclimatic changes will profoundly affect arctic communities.

Habitat fragmentation and species diversity in competitive communities.

Habitat loss is one of the key drivers of the ongoing decline of biodiversity. However, ecologists still argue about how fragmentation of habitat (independent of habitat loss) affects species richness. The recently proposed habitat amount hypothesis posits that species richness only depends on the total amount of habitat in a local landscape. In contrast, empirical studies report contrasting patterns: some find positive and others negative effects of fragmentation per se on species richness. To explain this apparent disparity, we devise a stochastic, spatially explicit model of competitive species communities in heterogeneous habitats. The model shows that habitat loss and fragmentation have complex effects on species diversity in competitive communities. When the total amount of habitat is large, fragmentation per se tends to increase species diversity, but if the total amount of habitat is small, the situation is reversed: fragmentation per se decreases species diversity.

Accounting for environmental variation in co-occurrence modelling reveals the importance of positive interactions in root-associated fungal communities.

Understanding the role of interspecific interactions in shaping ecological communities is one of the central goals in community ecology. In fungal communities, measuring interspecific interactions directly is challenging because these communities are composed of large numbers of species, many of which are unculturable. An indirect way of assessing the role of interspecific interactions in determining community structure is to identify the species co-occurrences that are not constrained by environmental conditions. In this study, we investigated co-occurrences among root-associated fungi, asking whether fungi co-occur more or less strongly than expected based on the environmental conditions and the host plant species examined. We generated molecular data on root-associated fungi of five plant species evenly sampled along an elevational gradient at a high arctic site. We analysed the data using a joint species distribution modelling approach that allowed us to identify those co-occurrences that could be explained by the environmental conditions and the host plant species, as well as those co-occurrences that remained unexplained and thus more probably reflect interactive associations. Our results indicate that not only negative but also positive interactions play an important role in shaping microbial communities in arctic plant roots. In particular, we found that mycorrhizal fungi are especially prone to positively co-occur with other fungal species. Our results bring new understanding to the structure of arctic interaction networks by suggesting that interactions among root-associated fungi are predominantly positive.

Finding flies in the mushroom soup: Host specificity of fungus-associated communities revisited with a novel molecular method.

Fruiting bodies of fungi constitute an important resource for thousands of other taxa. The structure of these diverse assemblages has traditionally been studied with labour-intensive methods involving cultivation and morphology-based species identification, to which molecular information might offer convenient complements. To overcome challenges in DNA extraction and PCR associated with the complex chemical properties of fruiting bodies, we developed a pipeline applicable for extracting amplifiable total DNA from soft fungal samples of any size. Our protocol purifies DNA in two sequential steps: (a) initial salt-isopropanol extraction of all nucleic acids in the sample is followed by (b) an extra clean-up step using solid-phase reversible immobilization (SPRI) magnetic beads. The protocol proved highly efficient, with practically all of our samples-regardless of biomass or other properties-being successfully PCR-amplified using metabarcoding primers and subsequently sequenced. As a proof of concept, we apply our methods to address a topical ecological question: is host specificity a major characteristic of fungus-associated communities, that is, do different fungus species harbour different communities of associated organisms? Based on an analysis of 312 fungal fruiting bodies representing 10 species in five genera from three orders, we show that molecular methods are suitable for studying this rich natural microcosm. Comparing to previous knowledge based on rearing and morphology-based identifications, we find a species-rich assemblage characterized by a low degree of host specialization. Our method opens up new horizons for molecular analyses of fungus-associated interaction webs and communities. Fruiting bodies of fungi constitute an important resource for thousands of other taxa. The structure of these diverse assemblages has traditionally been studied with labour-intensive methods involving cultivation and morphology-based species identification, to which molecular information might offer convenient complements. To overcome challenges in DNA extraction and PCR associated with the complex chemical properties of fruiting bodies, we developed a pipeline applicable for extracting amplifiable total DNA from soft fungal samples of any size. Our protocol purifies DNA in two sequential steps: (a) initial salt-isopropanol extraction of all nucleic acids in the sample is followed by (b) an extra clean-up step using solid-phase reversible immobilization (SPRI) magnetic beads. The protocol proved highly efficient, with practically all of our samples-regardless of biomass or other properties-being successfully PCR-amplified using metabarcoding primers and subsequently sequenced. As a proof of concept, we apply our methods to address a topical ecological question: is host specificity a major characteristic of fungus-associated communities, that is, do different fungus species harbour different communities of associated organisms? Based on an analysis of 312 fungal fruiting bodies representing 10 species in five genera from three orders, we show that molecular methods are suitable for studying this rich natural microcosm. Comparing to previous knowledge based on rearing and morphology-based identifications, we find a species-rich assemblage characterized by a low degree of host specialization. Our method opens up new horizons for molecular analyses of fungus-associated interaction webs and communities.

Assessing changes in arthropod predator-prey interactions through DNA-based gut content analysis-variable environment, stable diet.

Analysing the structure and dynamics of biotic interaction networks and the processes shaping them is currently one of the key fields in ecology. In this paper, we develop a novel approach to gut content analysis, thereby deriving a new perspective on community interactions and their responses to environment. For this, we use an elevational gradient in the High Arctic, asking how the environment and species traits interact in shaping predator-prey interactions involving the wolf spider Pardosa glacialis. To characterize the community of potential prey available to this predator, we used pitfall trapping and vacuum sampling. To characterize the prey actually consumed, we applied molecular gut content analysis. Using joint species distribution models, we found elevation and vegetation mass to explain the most variance in the composition of the prey community locally available. However, such environmental variables had only a small effect on the prey community found in the spider's gut. These observations indicate that Pardosa exerts selective feeding on particular taxa irrespective of environmental constraints. By directly modelling the probability of predation based on gut content data, we found that neither trait matching in terms of predator and prey body size nor phylogenetic or environmental constraints modified interaction probability. Our results indicate that taxonomic identity may be more important for predator-prey interactions than environmental constraints or prey traits. The impact of environmental change on predator-prey interactions thus appears to be indirect and mediated by its imprint on the community of available prey.

Protax-fungi: a web-based tool for probabilistic taxonomic placement of fungal internal transcribed spacer sequences.

Incompleteness of reference sequence databases and unresolved taxonomic relationships complicates taxonomic placement of fungal sequences. We developed Protax-fungi, a general tool for taxonomic placement of fungal internal transcribed spacer (ITS) sequences, and implemented it into the PlutoF platform of the UNITE database for molecular identification of fungi. With empirical data on root- and wood-associated fungi, Protax-fungi reliably identified (with at least 90% identification probability) the majority of sequences to the order level but only around one-fifth of them to the species level, reflecting the current limited coverage of the databases. Protax-fungi outperformed the Sintax and Rdb classifiers in terms of increased accuracy and decreased calibration error when applied to data on mock communities representing species groups with poor sequence database coverage. We applied Protax-fungi to examine the internal consistencies of the Index Fungorum and UNITE databases. This revealed inconsistencies in the taxonomy database as well as mislabelling and sequence quality problems in the reference database. The according improvements were implemented in both databases. Protax-fungi provides a robust tool for performing statistically reliable identifications of fungi in spite of the incompleteness of extant reference sequence databases and unresolved taxonomic relationships.

How to make more out of community data? A conceptual framework and its implementation as models and software.

Community ecology aims to understand what factors determine the assembly and dynamics of species assemblages at different spatiotemporal scales. To facilitate the integration between conceptual and statistical approaches in community ecology, we propose Hierarchical Modelling of Species Communities (HMSC) as a general, flexible framework for modern analysis of community data. While non-manipulative data allow for only correlative and not causal inference, this framework facilitates the formulation of data-driven hypotheses regarding the processes that structure communities. We model environmental filtering by variation and covariation in the responses of individual species to the characteristics of their environment, with potential contingencies on species traits and phylogenetic relationships. We capture biotic assembly rules by species-to-species association matrices, which may be estimated at multiple spatial or temporal scales. We operationalise the HMSC framework as a hierarchical Bayesian joint species distribution model, and implement it as R- and Matlab-packages which enable computationally efficient analyses of large data sets. Armed with this tool, community ecologists can make sense of many types of data, including spatially explicit data and time-series data. We illustrate the use of this framework through a series of diverse ecological examples.

How are species interactions structured in species-rich communities? A new method for analysing time-series data.

Estimation of intra- and interspecific interactions from time-series on species-rich communities is challenging due to the high number of potentially interacting species pairs. The previously proposed sparse interactions model overcomes this challenge by assuming that most species pairs do not interact. We propose an alternative model that does not assume that any of the interactions are necessarily zero, but summarizes the influences of individual species by a small number of community-level drivers. The community-level drivers are defined as linear combinations of species abundances, and they may thus represent e.g. the total abundance of all species or the relative proportions of different functional groups. We show with simulated and real data how our approach can be used to compare different hypotheses on community structure. In an empirical example using aquatic microorganisms, the community-level drivers model clearly outperformed the sparse interactions model in predicting independent validation data.

Novel community data in ecology-properties and prospects.

New technologies for monitoring biodiversity such as environmental (e)DNA, passive acoustic monitoring, and optical sensors promise to generate automated spatiotemporal community observations at unprecedented scales and resolutions. Here, we introduce 'novel community data' as an umbrella term for these data. We review the emerging field around novel community data, focusing on new ecological questions that could be addressed; the analytical tools available or needed to make best use of these data; and the potential implications of these developments for policy and conservation. We conclude that novel community data offer many opportunities to advance our understanding of fundamental ecological processes, including community assembly, biotic interactions, micro- and macroevolution, and overall ecosystem functioning.

Global Spore Sampling Project: A global, standardized dataset of airborne fungal DNA.

Novel methods for sampling and characterizing biodiversity hold great promise for re-evaluating patterns of life across the planet. The sampling of airborne spores with a cyclone sampler, and the sequencing of their DNA, have been suggested as an efficient and well-calibrated tool for surveying fungal diversity across various environments. Here we present data originating from the Global Spore Sampling Project, comprising 2,768 samples collected during two years at 47 outdoor locations across the world. Each sample represents fungal DNA extracted from 24 m3 of air. We applied a conservative bioinformatics pipeline that filtered out sequences that did not show strong evidence of representing a fungal species. The pipeline yielded 27,954 species-level operational taxonomic units (OTUs). Each OTU is accompanied by a probabilistic taxonomic classification, validated through comparison with expert evaluations. To examine the potential of the data for ecological analyses, we partitioned the variation in species distributions into spatial and seasonal components, showing a strong effect of the annual mean temperature on community composition.

Evaluating the predictive performance of presence-absence models: Why can the same model appear excellent or poor?

When comparing multiple models of species distribution, models yielding higher predictive performance are clearly to be favored. A more difficult question is how to decide whether even the best model is "good enough". Here, we clarify key choices and metrics related to evaluating the predictive performance of presence-absence models. We use a hierarchical case study to evaluate how four metrics of predictive performance (AUC, Tjur's R 2, max-Kappa, and max-TSS) relate to each other, the random and fixed effects parts of the model, the spatial scale at which predictive performance is measured, and the cross-validation strategy chosen. We demonstrate that the very same metric can achieve different values for the very same model, even when similar cross-validation strategies are followed, depending on the spatial scale at which predictive performance is measured. Among metrics, Tjur's R 2 and max-Kappa generally increase with species' prevalence, whereas AUC and max-TSS are largely independent of prevalence. Thus, Tjur's R 2 and max-Kappa often reach lower values when measured at the smallest scales considered in the study, while AUC and max-TSS reaching similar values across the different spatial levels included in the study. However, they provide complementary insights on predictive performance. The very same model may appear excellent or poor not only due to the applied metric, but also how predictive performance is exactly calculated, calling for great caution on the interpretation of predictive performance. The most comprehensive evaluation of predictive performance can be obtained by evaluating predictive performance through the combination of measures providing complementary insights. Instead of following simple rules of thumb or focusing on absolute values, we recommend comparing the achieved predictive performance to the researcher's own a priori expectations on how easy it is to make predictions related to the same question that the model is used for.

The role of seasonality in shaping the interactions of honeybees with other taxa.

The Eltonian niche of a species is defined as its set of interactions with other taxa. How this set varies with biotic, abiotic and human influences is a core question of modern ecology. In seasonal environments, the realized Eltonian niche is likely to vary due to periodic changes in the occurrence and abundance of interaction partners and changes in species behavior and preferences. Also, human management decisions may leave strong imprints on species interactions. To compare the impact of seasonality to that of management effects, honeybees provide an excellent model system. Based on DNA traces of interaction partners archived in honey, we can infer honeybee interactions with floral resources and microbes in the surrounding habitats, their hives, and themselves. Here, we resolved seasonal and management-based impacts on honeybee interactions by sampling beehives repeatedly during the honey-storing period of honeybees in Finland. We then use a genome-skimming approach to identify the taxonomic contents of the DNA in the samples. To compare the effects of the season to the effects of location, management, and the colony itself in shaping honeybee interactions, we used joint species distribution modeling. We found that honeybee interactions with other taxa varied greatly among taxonomic and functional groups. Against a backdrop of wide variation in the interactions documented in the DNA content of honey from bees from different hives, regions, and beekeepers, the imprint of the season remained relatively small. Overall, a honey-based approach offers unique insights into seasonal variation in the identity and abundance of interaction partners among honeybees. During the summer, the availability and use of different interaction partners changed substantially, but hive- and taxon-specific patterns were largely idiosyncratic as modified by hive management. Thus, the beekeeper and colony identity are as important determinants of the honeybee's realized Eltonian niche as is seasonality.

DNA traces the origin of honey by identifying plants, bacteria and fungi.

The regional origin of a food product commonly affects its value. To this, DNA-based identification of tissue remains could offer fine resolution. For honey, this would allow the usage of not only pollen but all plant tissue, and also that of microbes in the product, for discerning the origin. Here we examined how plant, bacterial and fungal taxa identified by DNA metabarcoding and metagenomics differentiate between honey samples from three neighbouring countries. To establish how the taxonomic contents of honey reflect the country of origin, we used joint species distribution modelling. At the lowest taxonomic level by metabarcoding, with operational taxonomic units, the country of origin explained the majority of variation in the data (70-79%), with plant and fungal gene regions providing the clearest distinction between countries. At the taxonomic level of genera, plants provided the most separation between countries with both metabarcoding and metagenomics. The DNA-based methods distinguish the countries more than the morphological pollen identification and the removal of pollen has only a minor effect on taxonomic recovery by DNA. As we find good resolution among honeys from regions with similar biota, DNA-based methods hold great promise for resolving honey origins among more different regions.

Fungal Communities Are Important Determinants of Bacterial Community Composition in Deadwood.

Fungal-bacterial interactions play a key role in the functioning of many ecosystems. Thus, understanding their interactive dynamics is of central importance for gaining predictive knowledge on ecosystem functioning. However, it is challenging to disentangle the mechanisms behind species associations from observed co-occurrence patterns, and little is known about the directionality of such interactions. Here, we applied joint species distribution modeling to high-throughput sequencing data on co-occurring fungal and bacterial communities in deadwood to ask whether fungal and bacterial co-occurrences result from shared habitat use (i.e., deadwood's properties) or whether there are fungal-bacterial interactive associations after habitat characteristics are taken into account. Moreover, we tested the hypothesis that the interactions are mainly modulated through fungal communities influencing bacterial communities. For that, we quantified how much the predictive power of the joint species distribution models for bacterial and fungal community improved when accounting for the other community. Our results show that fungi and bacteria form tight association networks (i.e., some species pairs co-occur more frequently and other species pairs co-occur less frequently than expected by chance) in deadwood that include common (or opposite) responses to the environment as well as (potentially) biotic interactions. Additionally, we show that information about the fungal occurrences and abundances increased the power to predict the bacterial abundances substantially, whereas information about the bacterial occurrences and abundances increased the power to predict the fungal abundances much less. Our results suggest that fungal communities may mainly affect bacteria in deadwood.IMPORTANCE Understanding the interactive dynamics between fungal and bacterial communities is important to gain predictive knowledge on ecosystem functioning. However, little is known about the mechanisms behind fungal-bacterial associations and the directionality of species interactions. Applying joint species distribution modeling to high-throughput sequencing data on co-occurring fungal-bacterial communities in deadwood, we found evidence that nonrandom fungal-bacterial associations derive from shared habitat use as well as (potentially) biotic interactions. Importantly, the combination of cross-validations and conditional cross-validations helped us to answer the question about the directionality of the biotic interactions, providing evidence that suggests that fungal communities may mainly affect bacteria in deadwood. Our modeling approach may help gain insight into the directionality of interactions between different components of the microbiome in other environments.

Computationally efficient joint species distribution modeling of big spatial data.

The ongoing global change and the increased interest in macroecological processes call for the analysis of spatially extensive data on species communities to understand and forecast distributional changes of biodiversity. Recently developed joint species distribution models can deal with numerous species efficiently, while explicitly accounting for spatial structure in the data. However, their applicability is generally limited to relatively small spatial data sets because of their severe computational scaling as the number of spatial locations increases. In this work, we propose a practical alleviation of this scalability constraint for joint species modeling by exploiting two spatial-statistics techniques that facilitate the analysis of large spatial data sets: Gaussian predictive process and nearest-neighbor Gaussian process. We devised an efficient Gibbs posterior sampling algorithm for Bayesian model fitting that allows us to analyze community data sets consisting of hundreds of species sampled from up to hundreds of thousands of spatial units. The performance of these methods is demonstrated using an extensive plant data set of 30,955 spatial units as a case study. We provide an implementation of the presented methods as an extension to the hierarchical modeling of species communities framework.