Evaluating Alternative Metacommunity Hypotheses for Diatoms in the McMurdo Dry Valleys Using Simulations and Remote Sensing Data

Diatoms are diverse and widespread freshwater Eukaryotes that make excellent microbial subjects for addressing questions in metacommunity ecology. In the McMurdo Dry Valleys of Antarctica, the simple trophic structure of glacier-fed streams provides an ideal outdoor laboratory where well-described diatom assemblages are found within two cyanobacterial mat types, which occupy different habitats and vary in coverage within and among streams. Specifically, black mats of Nostoc spp. occur in marginal wetted habitats, and orange mats (Oscillatoria spp. and Phormidium spp.) occur in areas of consistent stream flow. Despite their importance as bioindicators for changing environmental conditions, the role of dispersal in structuring dry valley diatom metacommunities remains unclear. Here, we use MCSim, a spatially explicit metacommunity simulation package for R, to test alternative hypotheses about the roles of dispersal and species sorting in maintaining the biodiversity of diatom assemblages residing in black and orange mats. The spatial distribution and patchiness of cyanobacterial mat habitats was characterized by remote imagery of the Lake Fryxell sub-catchment in Taylor Valley. The available species pool for diatom metacommunity simulation scenarios was informed by the Antarctic Freshwater Diatoms Database, maintained by the McMurdo Dry Valleys Long Term Ecological Research program. We used simulation outcomes to test the plausibility of alternative community assembly hypotheses to explain empirically observed patterns of freshwater diatom biodiversity in the long-term record. The most plausible simulation scenarios suggest species sorting by environmental filters, alone, was not sufficient to maintain biodiversity in the Fryxell Basin diatom metacommunity. The most plausible scenarios included either (1) neutral models with different immigration rates for diatoms in orange and black mats or (2) species sorting by a relatively weak environmental filter, such that dispersal dynamics also influenced diatom community assembly, but there was not such a strong disparity in immigration rates between mat types. The results point to the importance of dispersal for understanding current and future biodiversity patterns for diatoms in this ecosystem, and more generally, provide further evidence that metacommunity theory is a useful framework for testing hypotheses about microbial community assembly.

Diatoms are diverse and widespread freshwater Eukaryotes that make excellent microbial subjects for addressing questions in metacommunity ecology. In the McMurdo Dry Valleys of Antarctica, the simple trophic structure of glacier-fed streams provides an ideal outdoor laboratory where well-described diatom assemblages are found within two cyanobacterial mat types, which occupy different habitats and vary in coverage within and among streams. Specifically, black mats of Nostoc spp. occur in marginal wetted habitats, and orange mats (Oscillatoria spp. and Phormidium spp.) occur in areas of consistent stream flow. Despite their importance as bioindicators for changing environmental conditions, the role of dispersal in structuring dry valley diatom metacommunities remains unclear. Here, we use MCSim, a spatially explicit metacommunity simulation package for R, to test alternative hypotheses about the roles of dispersal and species sorting in maintaining the biodiversity of diatom assemblages residing in black and orange mats. The spatial distribution and patchiness of cyanobacterial mat habitats was characterized by remote imagery of the Lake Fryxell sub-catchment in Taylor Valley. The available species pool for diatom metacommunity simulation scenarios was informed by the Antarctic Freshwater Diatoms Database, maintained by the McMurdo Dry Valleys Long Term Ecological Research program. We used simulation outcomes to test the plausibility of alternative community assembly hypotheses to explain empirically observed patterns of freshwater diatom biodiversity in the long-term record. The most plausible simulation scenarios suggest species sorting by environmental filters, alone, was not sufficient to maintain biodiversity in the Fryxell Basin diatom metacommunity. The most plausible scenarios included either (1) neutral models with different immigration rates for diatoms in orange and black mats or (2) species sorting by a relatively weak environmental filter, such that dispersal dynamics also influenced diatom community assembly, but there was not such a strong disparity

INTRODUCTION
Metacommunity theory provides a framework to integrate dispersal dynamics with local community assembly processes, such as environmental filtering, over a patchy and heterogeneous landscape (Leibold et al., 2004;Leibold and Chase, 2017). In general, different levels of dispersal, habitat heterogeneity, and species sorting are predicted to result in different biodiversity outcomes (Leibold et al., 2004;Sokol et al., 2017). However, dispersal can be difficult to assess, especially for microorganisms, which have been long-assumed to have nearly unlimited, widespread dispersal (Baas Becking, 1934). The dogma that all microbes are ubiquitous has been challenged over the past two decades, and evidence is mounting that metacommunity theory applies to microbial life (Martiny et al., 2006;Nemergut et al., 2013;Langenheder and Lindström, 2019). Simulations can be used to overcome some of the challenges of understanding microbial community assembly by characterizing the likely roles of environmental filters (e.g., species sorting) and dispersal (e.g., source-sink dynamics) in structuring empirically observed metacommunities (Sokol et al., 2015) to compliment field-based experimental observations.
The ice-free desert oases of Antarctica represent ideal locations to test the application of metacommunity theory to microbial communities because these ecosystems are characterized by geographical isolation, relatively low anthropogenic impacts, and low trophic complexity. Here, freshwater habitats often harbor dense cyanobacterial mats, which are the major primary producers in the landscape . Because these ecosystems lack vascular plants, data products derived from satellite imagery (e.g., NDVI maps) can be used to map the distribution of such mats throughout the landscape . In the McMurdo Dry Valleys (MDV) of East Antarctica, "orange" and "black" mat types are common in and around glacial meltwater streams. These mats differ based on their appearance, dominant taxa, and location in the stream channel (Alger, 1997;McKnight and Tate, 1997;McKnight et al., 1998;Kohler et al., 2015;Van Horn et al., 2016). Orange mats are dominated by the filamentous cyanobacterial genera Oscillatoria and Phormidium, and typically occur in areas with consistent stream flow. In contrast, black mats are dominated by Nostoc spp. and occur in intermittently wet habitats, such as seeps and along stream margins. Furthermore, the coverage of mat growth varies dramatically within and among streams, where surveys of stream reaches have reported cyanobacterial mat percent coverages ranging from near 0 to near 100% (Alger, 1997;Kohler et al., 2015).
These cyanobacterial mats provide the primary habitat for diatom (Bacillariophyceae) assemblages in the ice-free landscape of the MDV (Esposito et al., 2006(Esposito et al., , 2008, where diatoms typically account for up to 10% of total mat biomass (Alger, 1997;McKnight et al., 1998;Kohler et al., 2015). Different diatom taxa are often associated with specific mat types (Darling et al., 2017;Andriuzzi et al., 2018;Kopalová et al., 2019). Thus, the heterogeneous distribution of mat types within and across MDV streams provides a patchy mosaic in which spatial habitat complexity can mediate the dispersal dynamics of these benthic diatoms, integrating local community assembly processes across the landscape, and creating a unique opportunity for study.
Here we focus on the resident diatom assemblages in the cyanobacterial mats because diatoms are excellent organisms for ecological monitoring and examination of biogeography due to their species-specific morphology, excellent preservation of their siliceous cell walls, and individual environmental tolerances (Smol and Stoermer, 2010;Spaulding et al., 2010). In the MDV, diatoms are an important sentinel taxonomic group because they are the most diverse and broadly distributed Eukaryote in the region (Adams et al., 2006;Spaulding et al., 2010;Kociolek et al., 2017). Thus, it is particularly important to understand the dynamics that control diatom biodiversity in the MDV.
Past work suggests that metacommunity theory can provide insight into the processes organizing MDV microbial biodiversity in general (Sokol et al., 2013) and diatom biodiversity in particular (Sakaeva et al., 2016). In MDV streams, several aspects of the flow regime appear to be dominant environmental filters for diatom communities. Specifically, the occurrence of flood pulses and periods without flow during the austral summer, as well as the interannual frequency of flow, have been linked to diatom assemblage composition (Esposito et al., 2006;Stanish et al., 2011Stanish et al., , 2012. Laboratory and field experiments indicate that optimal growth temperatures (Darling et al., 2017) and nutrient concentrations  may also influence species composition, but these factors are also strongly influenced by stream hydrology (Wlostowski et al., 2018). Thus, there is strong support that hydrology is a master environmental variable for understanding diatom assemblage structure in the MDV (Esposito et al., 2006;Stanish et al., 2012). In addition, dispersal dynamics can also influence benthic diatom community composition (Soininen et al., 2004;Verleyen et al., 2009;Heino et al., 2010;Sakaeva et al., 2016;Chen et al., 2019). Recently, Sakaeva et al. (2016) provided evidence that diatom assemblages residing in cyanobacterial mats at pond margins showed trends in composition across both spatial and environmental gradients in the MDV, indicating that dispersal dynamics may modify how environmental filters structure diatom assemblages across the landscape, aligning with other similar studies from the region (Michaud et al., 2012;Sokol et al., 2013).
Both wind and water are potential vectors for diatom dispersal in the MDV. Wind-driven transport has been documented for organic matter and biota in Taylor Valley Diaz et al., 2018), where the most extreme, and likely most consequential, foehn wind events occur during the austral winters (Nylen et al., 2004). Many MDV diatom genera are considered to be aerophilic (e.g., Luticola, Hantzschia, Humidophila, and Muelleria), and thus adapted to survive the challenges associated with both aeolian transport (Diaz et al., 2018) and the desiccation-rehydration and freeze-thaw cycles associated with life in MDV cyanobacterial mat habitats (McKnight et al., 1999). Streams have also been shown to be an important vector for the movement of organic matter, including diatoms, downstream toward the terminal lakes (Cullis et al., 2014;Kohler et al., 2018). Black mats are probably more easily mobilized than orange mats (e.g., by high wind or high flow events) and appear to be more sensitive to changes in the hydrologic regime (Kohler et al., 2015(Kohler et al., , 2018Gooseff et al., 2017). Thus, there is the potential for differences in the dispersal of diatoms among patches of cyanobacterial mats in the MDV landscape to mediate how local (i.e., "alpha-scale" or "patchscale") community assembly processes scale to maintain the biodiversity of the diatom metacommunity (i.e., "gamma-scale" or "regional-scale").
It is to be expected that diatom biodiversity in the MDV will be influenced by physical controls over the hydrologic habitat template (Poff and Ward, 1990), largely because hydrology controls the distribution of the cyanobacterial mats that provide habitat for diatom assemblages. Indeed, the past three decades of research in the MDV provide strong support for this general hypothesis. However, the ways in which an altered landscape translates to a shift in biodiversity can vary dramatically given different metacommunity dynamics (e.g., Sokol et al., 2015Sokol et al., , 2017. Thus, if we assume incorrect underlying community assembly mechanisms, our predictions about how future climatic conditions will affect MDV diatom biodiversity will be incorrect. Furthermore, as warming is expected in the coming decades (Shindell and Schmidt, 2004), a mechanistic understanding of the underlying diatom community assembly dynamics will be crucial for interpreting trends in the composition and diversity of these indicator taxa in a non-stationary world (Milly et al., 2008;Fitzpatrick et al., 2018).
Given that cyanobacterial mats provide a patchy mosaic of habitats for the MDV diatom metacommunity, our alternative hypotheses are: H1: Neutral model (NM) hypothesis: all taxa are functionally equivalent and respond similarly to environmental filters (sensu Hubbell, 2001), but orange and black mats represent patches characterized by different dispersal dynamics (i.e., one habitat type allows for a high colonization success rate for immigrants).
H2: Moderate species sorting (MSS) hypothesis: different mat types represent patches with both different dispersal dynamics and different local environmental filters that moderately affect species sorting.
H3: Strong species sorting (SSS) hypothesis: differences in environmental filters between mat types exert a very strong influence over species sorting, largely overwhelming any spatial patterns that would otherwise emerge from dispersal dynamics.
All three of these hypotheses provide different mechanisms for the maintenance of diatom biodiversity in the MDV, and different processes by which an altered landscape might drive changes in diatom biodiversity in future climate scenarios. H1 provides a mechanism based on NM metacommunity dynamics where the stochastic processes that determine emergent diatom biodiversity patterns are governed by the spatial arrangement of habitable patches in the MDV landscape and the movement of diatoms among these patches (i.e., cyanobacterial mats). Alternatively, H3 provides a deterministic mechanism based on SSS dynamics where environmental filters sort diatom assemblages by functional traits at the patch (i.e., a cyanobacterial mat) scale. Following H3, emergent biodiversity patterns are simply predicted by the relative dominance of different mat types, and thus different types of environmental filters, in the MDV landscape. H2 (MSS) represents a hybrid hypothesis allowing for stochastic dispersal and colonization dynamics (H1) to play out in a non-neutral landscape (H3). Given existing evidence, we expect both species sorting and dispersal dynamics to play roles in organizing the diatom metacommunities residing in cyanobacterial mats in and around MDV streams. Modeling provides an opportunity to characterize what those stochastic (dispersal) and deterministic (species sorting) dynamics might look like.
Here, our objectives are (1) to test the plausibility of the three alternative metacommunity hypotheses described above to explain observed patterns in diatom biodiversity in cyanobacterial mats in stream channels in the Fryxell Basin of Taylor Valley, MDV and (2) to characterize the stochastic (dispersal) and deterministic (environmental filtering) dynamics of the most plausible scenarios. To evaluate the plausibility of the alternative hypotheses, biodiversity must be characterized at both the patch (alpha-diversity) and metacommunity (gammadiversity) scales in order to understand its drivers. Among-patch variation in community composition (beta-diversity) provides a measure of scaling, i.e., the contribution of each local patch to the regional species pool. We apply a spatially explicit numerical model, MCSim (Sokol et al., 2015Sokol, 2019), to predict diatom biodiversity outcomes under different metacommunity scenarios, using satellite imagery (i.e., a derived map of cyanobacterial mats) and ground-based observations of diatom assemblages to constrain simulation initial conditions (Figure 1). We use simulation outcomes to evaluate the plausibility of each of the three hypotheses outlined above.

Site Description
The MDV are approximately 4,500 km 2 in area, making them the largest ice-free area on the continent (Levy, 2013). The MDV FIGURE 1 | Workflow to incorporate source data from the McMurdo Long Term Ecological Research (MCM LTER) program and remotely sensed imagery to constrain spatially explicit metacommunity simulations based on the MCSim framework to assess the plausibility of different metacommunity scenarios. Simulations that best maintained initial biodiversity levels were considered the most plausible. provide an excellent model ecosystem for integrating community assembly theory and models with patterns of diatom distribution observed in the field because the diatom metacommunity presents a tractable number (∼50) of taxa (Esposito et al., 2008) whose distribution and ecological preferences have been rigorously described in streams since 1994 by the McMurdo Dry Valleys Long Term Ecological Research (MCM LTER) program (Esposito et al., 2006;Stanish et al., 2011Stanish et al., , 2012Kohler et al., 2016;Darling et al., 2017).
During the austral summer, the MDV streams are fed by glacial meltwater and are underlain by permafrost, receiving no groundwater inflow, and flow into perennially ice-covered lakes (Conovitz et al., 1998). The active layer thaws to a depth of about 60 cm and water in the stream channel exchanges with water in saturated sediments underlying and adjacent to the streambed, referred to as the hyporheic zone (Conovitz et al., 2006). This melt period typically lasts for 4-10 weeks, and provides a short growing season when cyanobacterial mats become biologically active (McKnight and Tate, 1997). Due to the influence of the adiabatic lapse rate, the elevation of the source glacier or snowfield exerts a strong influence on the flow regime for a given stream. During cold cloudy summers meltwater generation for the higher elevation glaciers is limited and storage of meltwater in the hyporheic zone further limits flow in the downstream reaches of the longer streams (Conovitz et al., 1998). During warm sunny summers, high flows can occur in most streams and flood pulses can act to scour the mats (Cullis et al., 2014). Overall, these differences in flow regime among streams influence the abundance and coverage of the different mat types (Kohler et al., 2015).
The water chemistry in MDV streams has been monitored by the MCM LTER program for over three decades and is generally dilute compared to temperate streams. Higher concentrations of major ions and silica occur in the longer streams driven by the greater extent of weathering reactions occurring in the hyporheic zone as the water flows from the glacier to the lake (Gooseff et al., 2002). The primary sources of nutrients for the cyanobacterial mats are nitrate deposited on the glacier surfaces associated with auroral activity, weathering of apatite releasing phosphorus in hyporheic sediments, and recycling of nutrients from mat biomass that has become entrained in hyporheic sediments (McKnight et al., 2004;Kohler et al., 2018). In streams with abundant mats, nutrient concentrations are low ( Table 1) compared with streams where geomorphic conditions, such as unstable sandy substrates, do not support mat growth (McKnight et al., 2004).
In this study, our aim was to model the diatom metacommunity dynamics among patches of cyanobacterial mats located in the sub-catchments of eight streams that drain meltwater from alpine glaciers into Lake Fryxell, which is one of three terminal lakes located in Taylor Valley, in the MDV (Table 1 and Figure 2). The Lake Fryxell drainage basin was selected as our focal system because the hydrological and biological characteristics of its streams have been well documented by the MCM LTER (e.g., Kohler et al., 2015;Wlostowski et al., 2016;Gooseff et al., 2017). Four of the streams are relatively short and have less variable hydrology (Kohler et al., 2015;Wlostowski et al., 2016). These include Canada Stream and Huey, which flow into Lake Fryxell from the north, and Bowles and Green Creek from the west. Canada Glacier and Lake Fryxell provide  geographic barriers between the two pairs of streams. The other four streams (Delta, Harnish, Crescent, and Von Guerard) have much longer flow paths, which result in more variable hydrographs, and drain meltwater from alpine glaciers from the south into Lake Fryxell (Kohler et al., 2015;Wlostowski et al., 2016). Streams were grouped such that cyanobacterial mat habitats located in the same stream group would be located in catchments with similar hydrologic characteristics and are not separated by any significant geographic barriers (i.e., Lake Fryxell or Canada Glacier). Stream group 1 (SG1) includes Canada Stream and Huey; stream group 2 (SG2) includes Bowles and Green; and stream group 3 (SG3) includes Delta, Harnish, Crescent, and Von Guerard. These groupings were used in our data sampling procedure (detailed below) to create a more realistic initial metacommunity to serve as a common simulation starting point.

MCSim Metacommunity Modeling
We used v0.4.9 of the MCSim metacommunity simulation package for R (Sokol et al., 2015Sokol, 2019) to model scenarios that represent H1, H2, and H3. The simulations required three general components to run (Figure 1, MCSim inputs): (1) a landscape map, (2) an initial metacommunity, and (3) parameters that determine the dispersal and recruitment dynamics that occur with each time step of the simulation. We used remotely sensed imagery from the WorldView-2 satellite (Figure 2)  and data from the MCM LTER monitoring record (Figure 3 and Supplementary Data Sheet 1) to create a simulation map ( Figure 4A and Supplementary Data Sheet 2) and an initial contrived diatom metacommunity (Supplementary Data Sheet 3) as a common starting point for all simulation scenarios. We created 10,000 independent simulations to explore how well different dispersal dynamics crossed with NM, MSS, and SSS metacommunity dynamics could maintain initial metacommunity composition.

The Simulation Landscape: Based on Remote Imagery
Satellite imagery (Figure 2) was used to derive a two-dimensional map ( Figure 4A) highlighting the distribution of photosynthetic mat communities throughout the watershed of Lake Fryxell (i.e., Fryxell Basin). Recent studies have demonstrated that satellite data and the use of common spectral parameters like the Normalized Difference Vegetation Index (NDVI) can effectively identify the distribution of these cyanobacterial mats in the MDV (Salvatore, 2015;Salvatore et al., 2020). NDVI is a commonly used spectral parameter in remote sensing that is used to identify the distribution and relative photosynthetic capacity of vegetation.
Most photosynthetic vegetation contains dark pigments that absorb solar radiation at visible wavelengths, while the lack of absorption in the near-infrared results in a significant increase in reflectance at wavelengths just beyond the visible. These contrasting spectral signatures are unique to photosynthetic species relative to nearly all other natural materials. More specifically, NDVI is a ratio between the difference in reflectance at near-infrared and visible wavelengths normalized to the sum of the reflectance at these wavelengths ([NIR-Red]/[NIR + Red], where red and NIR wavelength bands are centered at 0.659 and 0.831 µm, respectively). The product of this ratio is a value between −1 and +1, with higher values representing greater mat density, photosynthetic activity, or the presence of otherwise spectrally dominant vegetation. While NDVI values are nonunique in most terrestrial ecosystems due to the heterogeneity in photosynthetic species, higher NDVI values are far more closely linked to the presence and areal coverage of cyanobacterial mats in the MDV due to the relative simplicity of these ecosystems . Here, we specifically use NDVI to identify likely patches of cyanobacterial mat communities in the Fryxell Basin to provide a proxy for the spatial distribution of diatom habitats. High-resolution WorldView-2 data were processed in accordance with the methods described in Salvatore (2015) and Salvatore et al. (2020). Data were calibrated to surface reflectance using the Dark Object Subtraction and Regression (DOS-R) methods before NDVI was derived. Surface reflectance represents a fundamental property of the surface itself and, in theory, is independent of other external influences. These calibration steps ensured that differences in the surface composition and coverage were the only properties being measured and not, for example, relative changes in atmospheric properties or illumination conditions. The derived NDVI raster image ( Figure 2C) was then imported into R and processed to only include pixels located in Fryxell Basin itself, minimizing any influences due to bedrock composition, illumination geometry, or topographic effects (R code available in Sokol et al., 2020). We then selected pixels with the 500 greatest NDVI values (summarized in Table 2) to represent 500 patches in the simulated metacommunity landscape ( Figure 4A). Thus, patch locations in the simulated landscape were based on the pixel coordinates in the source NDVI raster image. Further, we assumed that cyanobacterial mat areal coverage within a pixel was positively correlated with the pixel's NDVI score , and that the size of a diatom assemblage within a pixel (J L ) scaled monotonically with the available habitat (i.e., mat coverage) in the pixel. Following these assumptions, total diatom assemblage size (J L , i.e., the count of extant individuals in an assemblage at a location) in each of the 500 patches included in a simulated metacommunity were a function of observed NDVI score, such that values for J L followed a log-normal distribution with minimum value of 50 and maximum value of 10,000. These bounds for J L were set based on a prior sensitivity analysis of how simulated biodiversity metrics were affected by both J L and overall metacommunity size (J M ) , and were chosen to allow for reasonable compute times, but still produce metacommunities that were sufficiently large such that simulation behavior was consistent with previously published studies (Sokol et al., 2015. While the metacommunity map, including patch location and patch carrying capacity (J L ), were determined from remote imagery, we were not able to derive other patch attributes in the same way. Specifically, we do not currently have a workflow to derive mat type (orange or black) or local environmental characteristics from available remote imagery. Thus, to create simulation scenarios in which patch attributes reflected realistic distributions of mat type and local environmental characteristics, we used a stratified sampling procedure to assign attributes from the MCM LTER stream monitoring records (described below) to patches in the simulated metacommunities.

Initial Conditions: Based on Empirical Observations
Simulation initial conditions were derived from the long-term stream monitoring data maintained by the MCM LTER program (see Stanish et al., 2011Stanish et al., , 2012. The source dataset, assembled from data available on the Antarctic Freshwater Diatoms Database 1 and the MCM LTER data portal 2 , spans 1994-2013 and includes 189 records of data from samples (mat cores) collected from orange and black cyanobacterial mats in and near stream channels in Fryxell Basin (summarized in Table 3, raw data available in Supplementary Data Sheet 1). Each record includes mat type (orange or black), a measurement of cyanobacterial mat chlorophyll-a standing biomass, and estimates of the relative abundances of observed diatom taxa.
To construct the patch attributes (see Supplementary Data Sheet 2) and diatom assemblages for the initial contrived metacommunity (see Supplementary Data Sheet 3) to serve as a starting point for the simulations, we randomly sampled, with replacement, 500 records from the source dataset derived from the MCM LTER monitoring program (see Supplementary Data FIGURE 3 | Summary of relative abundances (mean ± sd) of diatom taxa in the MCM LTER record for samples collected from orange and black cyanobacterial mats for each "stream group". Morphospecies designations in the source database were preserved. Sheet 1), described above. The data from each random sample were assigned to a patch in the simulation map ( Figure 4A), thus, defining mat type, local environmental characteristics, and the relative abundances of the starting diatom assemblage for each patch in the simulated metacommunity. The sampling procedure was stratified by stream group, such that only empirically observed records from SG1 streams were assigned to the patches in SG1 on the simulation map, only records from SG2 streams were assigned to the SG2 patches on the simulation map, and only records from SG3 streams were assigned to the SG3 patches on the simulation map. This sampling procedure maintained a relative dominance by mat type, a distribution of chlorophyll-a standing biomass values, and diatom community compositions similar to what was observed to occur in each group of subcatchments (Figure 3).
Using the ade4 package for R (Chessel and Dufour, 2006), we conducted an Outlying Mean Index (OMI) analysis (Dolédec et al., 2000) using the 189 records from the MCM LTER FIGURE 4 | All metacommunity simulations used the same map, where J L indicates community size at each patch, x-and y-axes indicate patch coordinates in pixels, which are the spatial units used in the simulation and correspond to a ∼2 m × 2 m area in the original raster map (A), each simulation scenario was randomly assigned one of the dispersal kernels (see Eq. 2) plotted in (B), and included species habitat preferences and niche breadths resulting in different types of metacommunity dynamics: (C-i), where C niche = 1, strong species sorting (SSS); (C-ii), where C niche = 4, moderate species sorting (MSS); or (C-iii), where C niche = 20, neutral model (NM). The "rug" (upward ticks along the x-axis) in (B) indicates all the pairwise distances between patches in the simulation landscape, and the "rug" in (C) is the distribution of patches along the environmental gradient used in the simulation.
Frontiers in Ecology and Evolution | www.frontiersin.org 8 September 2020 | Volume 8 | Article 521668 source data to construct a univariate environmental gradient (E) and diatom trait matrix. Specifically, black mat occurrence and log(x + 1) transformed chlorophyll-a standing biomass measurements were scaled and used in a PCA to create a composite variable (E) to serve as an environmental state variable. We used these measures of chlorophyll-a standing biomass in cyanobacterial mat cores from the MCM LTER record as a proxy for local hydrologic characteristics (Kohler et al., 2015). We included mat type (i.e., a binary dummy variable indicating if a mat was "black" or "orange") to provide a mechanism by which mat color could act as an environmental filter on species composition (i.e., an alternative to mat type affecting immigration rate, m). We chose not to include water chemistry data, as previous studies have demonstrated such data are not major predictors of diatom community composition in this system (Esposito et al., 2006;Stanish et al., 2012). The E scores calculated for each MCM LTER record were mapped onto their respective patches in the simulation (see Supplementary Data Sheet 2). E scores were then used in an OMI analysis along with Hellinger transformed (Legendre and Gallagher, 2001)

Metacommunity Model Details
The MCSim metacommunity model is explained in detail elsewhere (Sokol et al., 2015, but briefly, it is a spatially explicit, zero-sum lottery model based on Hubbell's neutral model (Hubbell, 2001), modified following Gravel et al. (2006) to be spatially explicit, include a dispersal kernel term, W(r), and the influence of an environmental filter, λ(E). Similar to the neutral model, MCSim parameters include m and ν, where m represents a relativized immigration rate at a patch in the metacommunity and ν is the probability a non-extant taxon will be recruited in a particular lottery event. The value of m for a patch with J L individuals is related to the number of immigrants, I, such that The dispersal kernel is defined as: where W(r) is the probability of a successful dispersal event, r is the distance between the source and sink patch in the landscape, and w is the slope of the dispersal kernel. Larger values for w create steeper dispersal kernels and increase the effect of dispersal limitation (e.g., Figure 4B). MCSim uses W(r) to weight the contributions of the diatom communities from each of the neighboring patches to the immigrant pool for any location in the landscape. Then m and ν (each take a value between 0 and 1) represent the proportion of individuals that are recruited from neighboring patches within the metacommunity, and from a pool of potential invaders from outside the modeled metacommunity (e.g., propagules transported by wind from other lake basins located up-valley), respectively. Given the relative contributions of the novel invader pool (defined by ν), the immigrant pool (defined by m), and the extant resident community to the recruitment pool, the environmental filter, λ(E), further weights the recruitment probabilities of each taxon, following where λ(E) is a taxon's re-weighted effective relative abundance in the recruitment pool, E represents the value of a local environmental state variable, µ is the taxon's optimal value for E, σ is the taxon's niche-breadth, and C niche is a constant that can modify the strength of an environmental filter. When C niche is increased for all taxa in a metacommunity, it weakens the influence of the environmental filter uniformly across all taxa, while preserving cross-taxa heterogeneity in niche breaths (e.g., Figure 4C). Given a set of coordinates describing the location and attributes of patches in a landscape (e.g., Figure 4A and Supplementary Data Sheet 2) and a table of initial taxon abundances at each location (see Supplementary Data Sheet 3), MCSim allows metacommunity dynamics to play out over a series of generations. For a metacommunity with J M total individuals, J M lottery recruitment events occur with each time step. Each recruitment event is stochastic, but the probability of each taxon being recruited is weighted by its relative abundance in the effective recruitment pool as determined by the interaction of all the terms described above (also see Sokol et al., 2017).

Simulation Experiment
We simulated 10,000 independent metacommunities using random combinations of the parameter values outlined in  Invasion rate (ν), probability of novel taxa invading the metacommunity 10 −6 , 10 −5 , 10 −4 , 10 −3 , 10 −2 , 10 −1 See Supplementary Data Sheet 4 for a complete design matrix for how parameter values were applied across the 10,000 simulations used in this study. *see Figure 3C for more information on how different niche breadth coefficients were implemented.
design matrix, Sokol et al., 2020 for R code to run the simulations). Each simulation used the same map ( Figure 4A and Supplementary Data Sheet 2) and initial metacommunity composition (Figure 3 and Supplementary Data Sheet 3) and ran for a duration of 100 time-steps. This approach provided a relatively balanced exploration of how simulation outcomes vary across the metacommunity parameter space represented in Table 4. Choices for simulation duration and parameter ranges were based off of previous work (Sokol et al., 2015, where we found simulations with a duration of >30-50 time steps usually reached a steady state with respect to measures of alpha, beta, and gamma diversity. Preliminary model runs were used to determine that the values chosen for C niche would create SSS, MSS, and NM metacommunity dynamics ( Figure 4C). An important modification to how MCSim was applied in this study is that immigration rate (m) was based on mat type. Within each simulation all orange mats had immigration rates of m orange and all black mats had immigration rates of m black . The values of m orange and m black were varied across simulation scenarios (Table 4 and Supplementary Data Sheet 4), and immigration disparity between mat types in each simulation was summarized as log(m black /m orange ). Thus, log-ratio immigration disparity provides a measure of complexity in dispersal dynamics that can be compared across simulations, where negative values indicate higher immigration rates in orange mats, positive values indicate higher immigration rates in black mats, and values near zero indicate uniform immigration (less complex dispersal dynamics) across all mat types.

Evaluating Simulation Plausibility
The MDV ecosystem is considered "metastable" (Gooseff et al., 2017) and sediment cores from Lake Fryxell indicate the regional diatom species pool has been relatively stable throughout the Holocene (Whittaker et al., 2008;Konfirst et al., 2011), but diatom assemblage composition has been shown to vary through time at the patch scale (Esposito et al., 2006;Stanish et al., 2011). Thus, we considered a suitable criterion for a "plausible" simulation to be a scenario that maintains richness and evenness at both the patch (alpha) and metacommunity (gamma) scales but can allow species composition to turnover as time progresses. To this end, we used multiplicative diversity partitioning (Jost, 2007) to calculate mean alpha-scale diversity (D α ), gammascale diversity (D γ , calculated using all patches in a simulated metacommunity), and compositional turnover among patches (D β ) in each simulated metacommunity, where D β = D γ /D α . We used order q = 1 species equivalents because this version of the metric incorporates both richness and evenness. Diversity partitions were calculated for the first and last time steps of each simulation with the vegetarian package for R (Charney and Record, 2012), using equal weights for all local communities because original source data available in the MCM LTER database were effectively relative abundance data. To quantify deviation from the reference condition for each of the simulations, we used a modification of a X 2 statistic (Sokol et al., 2015): where χ 2 sim provides a multiscale measure of change in biodiversity. Specifically, this statistic provides a measure of how much D α , D β , and D γ at the final time step (t = 100) deviate from initial values for a given simulation, sim. Squared differences are each standardized by the standard deviation of diversity measured at that scale at the final time step across all simulations, so that deviation at each scale contributes equally to the χ 2 sim metric. Thus, simulations with the smallest values represent those in which D α , D β , and D γ changed the least during the course of the simulation.
Among the ∼3333 simulations for each metacommunity type, we considered those with χ 2 sim scores in the 5th percentile (i.e., the scenarios that best maintained initial diversity) to represent the most plausible metacommunity simulations. We used this approach to identify the simulations with the optimal dispersal kernel slope (w), immigration rates (m black and m orange ), and invasion rate (ν) for each metacommunity type. Thus, comparing the top 5% most plausible simulation outcomes (those with χ 2 sim scores in the 5th percentile) for each metacommunity type provides a more apples-to-apples comparison between metacommunity types because the dispersal and immigration parameters have been optimized independently for each metacommunity type.
To quantify the change in taxonomic composition during the course of each simulation, we calculated patch-scale temporal beta-diversity. For each patch in a simulated metacommunity, we calculated the Bray-Curtis dissimilarity between the extant diatom assemblages at the first (t = 1) and last (t = 100) time steps using Hellinger transformed abundances (Legendre and Gallagher, 2001). We took the median Bray-Curtis dissimilarity as a representative measurement of patch-scale temporal betadiversity for each simulation. We used the decostand and vegdist functions in the vegan package for R (Oksanen et al., 2017) to transform the data and calculate dissimilarity metrics, respectively.

RESULTS
The dataset that was simulated for this analysis includes 10,000 independent metacommunities, each with a duration of 100 time-steps (i.e., generations). The simulation output provides the information necessary to explore how local and metacommunity-scale trends in biodiversity vary over the region of metacommunity parameter space described in Table 4. Simulations were roughly evenly split across the neutral model (NM), moderate species sorting (MSS), and strong species sorting (SSS) metacommunity types ( Table 4). Simulations with both NM and MSS metacommunity dynamics produced outcomes with lower χ 2 sim scores than SSS simulations (Figure 5), indicating that NM and MSS type metacommunity dynamics were better able to maintain D α , D β , and D γ than SSS metacommunity dynamics within the metacommunity simulation framework that was used for this study.
In comparing the top 5% scenarios across metacommunity types, it is clear that NM and MSS metacommunity types, both with χ 2 sim scores < 1, are each capable of providing much more plausible dynamics than SSS scenarios (min χ 2 sim = 1.99) for diatoms in Fryxell Basin (Figure 5).
Overall, all three metacommunity types produced scenarios that maintained gamma diversity (D γ ∼ 22 species equivalents, Figure 6). Additionally, simulations with both MSS and NM metacommunity dynamics were capable of maintaining values of D α and D β similar to initial conditions, which were 10.6 and 2.07 species equivalents, respectively. However, none of the scenarios with SSS dynamics were able to maintain D α , and therefore produced metacommunities with inflated values of D β by the final time step.
Additionally, NM and MSS scenarios better maintained initial community composition than SSS scenarios (Figure 7). Patchscale temporal beta-diversity is a metric that is bounded between 0 and 1, and represents the median change in diatom assemblage taxonomic composition during the course of a simulation. Patch-scale temporal beta-diversity was much greater for SSS scenarios (>0.60) than NM and MSS scenarios (<0.50). However, even the most plausible scenarios experienced some taxonomic FIGURE 6 | Diversity outcomes for the final time step (t = 100) of the simulated metacommunities. Diversity partitions follow the Jost (2007) order q = 1 multiplicative method to calculate diversity at alpha (patch), beta (turnover among patches), and gamma (metacommunity) scales. Outcomes are displayed separately for scenarios with strong species sorting (SSS), moderate species sorting (MSS), and neutral model (NM) metacommunity dynamics. Dashed lines indicate initial diversity partitions at t = 1, which were the same for all 10,000 simulations. turnover, as the lowest observed patch-scale temporal betadiversity was still >0.30.
A comparison of the distribution of dispersal parameter values for the "most plausible" simulated metacommunities (i.e., with 5th percentile χ 2 sim scores) against the full population of ∼3333 simulations provided insight into the dispersal dynamics necessary to maintain diversity for each metacommunity type (Figure 8). For NM metacommunities, the most plausible scenarios had higher immigration rates in orange mats than black mats (i.e., m orange > m black ) and low rates of invasion of non-extant taxa (low values of ν), relative to the full population of all ∼3333 NM simulations. On the other hand, there was no difference in the distribution of dispersal kernel slopes (w) between the top 5% and the full population of NM simulations. For MSS simulations, the top 5% simulations had a bias toward higher immigration rates in black mats than orange mats (i.e., m black > m orange ), more moderate invasion rates (ν), and lower dispersal kernel slopes (w) than the full population of all ∼3333 MSS simulations. As described above, no SSS metacommunity simulations were capable of maintaining initial biodiversity, however, the top 5% SSS simulations that deviated the least from initial conditions had high values of ν and dispersal kernels with shallow slopes (small values for w) relative to the full population of ∼3333 SSS simulations.

DISCUSSION
Here we tested the plausibility of alternative metacommunity hypotheses to explain controls over diatom biodiversity in cyanobacterial mats residing in Fryxell Basin, MDV, using spatially explicit, process-based simulations (Figure 1). By using the NDVI metric derived from a remotely sensed image of Fryxell Basin to create the simulation map, our simulated landscape provides a realistic model of the patchy distribution of cyanobacterial mats in which the diatom assemblages reside. Because the results are based on simulations with known differences in the underlying species sorting and dispersal parameters ( Table 4 and Supplementary Data Sheet 4), outcomes can explicitly link spatial biodiversity patterns to specific metacommunity dynamics. Similar to inferences made from empirical observations by Sokol et al. (2013) and Sakaeva et al. (2016), the most plausible simulation outcomes in this study suggest species sorting by environmental filters, alone, does not provide a suitable mechanism to maintain microbial biodiversity in the MDV. Rather, simulations that included mechanisms for dispersal to influence the composition of the recruitment pool at each patch in the landscape provided the best candidate models for the dynamics that govern the diatom metacommunity in Fryxell Basin.
In this study, a metacommunity simulation was considered "plausible" if it was capable of maintaining both patch-scale (alpha) and metacommunity-scale (gamma) diversity (thus, also beta-diversity) similar to what was observed in the MCM LTER record for Fryxell Basin. NM and MSS metacommunity simulations both offered plausible candidate models that maintained diversity over time, but the dispersal characteristics in the most plausible NM simulations were different than those in the most plausible MSS simulations. It is important to note that the composition of diatom assemblages can still change over time in simulations that maintain diversity. Even the most plausible NM simulations generally showed moderate change in composition over time (patch-scale temporal beta-diversity was >0.30) because NM dynamics allow for drift in species composition (Chase, 2007). We did not use temporal trends in composition to evaluate the plausibility of alternative scenarios because we know the assumption of invariant community composition over time would be incorrect (Esposito et al., 2006;Stanish et al., 2011). The best simulation scenarios identified in this study (summarized in Figures 5-8) provide plausible and testable explanations about the types of metacommunity dynamics that might explain how spatial variation in stream hydrology and cyanobacterial mat coverage are linked to emergent patterns in diatom biodiversity in the MDV.

Ecological Interpretation of Plausible NM Simulations
The most plausible NM simulations suggest both regional diversity (gamma-scale) and patch level (alpha-scale) diversity can be maintained in the Fryxell Basin diatom metacommunity (Figure 6) if (1) there is a disparity in immigration rates between different mat types (higher immigration in orange mats) and (2) if the probability of invasion of novel taxa from outside the metacommunity (e.g., from other terminal lake basins located up-valley) is low (Figure 8). This outcome suggests that adding complexity to dispersal dynamics in the diatom FIGURE 7 | The distribution of patch-scale temporal beta-diversity (x-axis) for all simulations (red) and the top 5% most plausible simulations (blue) for each metacommunity type: neutral model (NM, top), moderate species sorting (MSS, middle), and strong species sorting (SSS, bottom). Patch-scale temporal beta-diversity was quantified for each simulation as the median of the Bray-Curtis dissimilarities between the initial and final (t = 100) time steps at each patch, using Hellinger transformed abundances. metacommunity (i.e., different immigration rates in different mat types) can maintain biodiversity, and the functional complexity required for species sorting dynamics is not necessary for the maintenance of biodiversity. However, complexity in dispersal dynamics alone was not sufficient because NM simulations in which the relative magnitude of immigration rates were reversed (i.e., black mats had higher immigration rates than orange mats) were not equally plausible. This outcome suggests that both the spatial arrangement and relative abundance of patch types with different immigrant fluxes (i.e., values of m in the model) are important characteristics that govern biodiversity at both local and regional scales.
This result provides a potential mechanism by which hydrologic controls over mats in Antarctic streams can control diatom biodiversity at both local and regional scales by modifying bulk diatom dispersal and colonization dynamics. Indeed, stream harshness, characterized by increased flow intermittency in addition to diel and seasonal flow variability, has been linked to the characteristics of both the MDV stream cyanobacterial mats (e.g., orange vs. black) (Kohler et al., 2015) and their resident diatom communities (Esposito et al., 2006;Stanish et al., 2011Stanish et al., , 2012. Because the most plausible NM simulations had larger values of m for orange mats than black mats, immigrants had a higher success rate for establishment in orange mats than in black mats in these scenarios. We interpret this property of the most plausible NM scenarios to indicate there is likely some characteristic of orange mat habitats that make them more suitable for colonization by immigrant diatoms dispersing from neighboring patches. Alternatively, simulation outcomes suggest recruitment from one generation to the next in black mats is more likely to be biased toward resident taxa that are already established in their respective mats (i.e., priority effects) (Chase, 2003). Past empirical work has shown that orange and black mats have different physical characteristics that might relate to diatom dispersal and colonization dynamics in this system. Orange mats are typically located in stream channels and have greater standing biomass and coverage (McKnight and Tate, 1997;McKnight et al., 1998), whereas black mats are located along wetted margins and more loosely attached to the substrate, and thus, more likely to be dislodged during high wind and/or high flow events (Kohler et al., 2015). Therefore, during high wind and/or flow events, black mats may be more likely to be sources of emigration and orange mats may be more likely to be immigrant sinks. Further, simulation outcomes suggest relative dominance of black mat communities in the landscape may determine the importance of priority effects in the metacommunity.
Future empirical and modeling studies should test if and how these mat characteristics predict colonist success versus resident taxon success. For example, do orange mats indeed experience higher immigrant fluxes? If so, is it because orange mats are often located in thalweg where they are bathed in a stream of diatoms dispersing via drift? Alternatively, do orange mats just provide a more suitable habitat for a wider variety of diatoms to successfully colonize, irrespective of whether they are wind-dispersed or drifting in the water column? Would repeated observations reveal that priority effects are indeed stronger in black mat diatom assemblages? Regardless, the NM simulation outcomes suggest that by controlling the distribution of different mat types, physical processes (e.g., hydrology) can effectively control the dispersal dynamics of the diatom metacommunity by way of controlling the areal extent and spatial arrangement of patch types with different immigrant fluxes. If this mechanism provides a major control over the structure of the diatom metacommunity, then it also suggests local and regional diatom biodiversity will be sensitive to hydrologically driven shifts in the presence of and coverage of different mat types predicted to occur as the MDV climate changes (Gooseff et al., 2017).

Ecological Interpretation of Most Plausible MSS Simulations
Importantly, NM scenarios did not provide the only plausible metacommunity dynamics for this study system. Simulations with MSS dynamics were nearly as successful as NM scenarios at maintaining diatom diversity at both the patch (alpha) and metacommunity (gamma) scales (Figure 6). However, in contrast with the most plausible NM simulations, the most plausible MSS simulations (1) had a small bias toward higher immigration rates in black mats, (2) had higher rates of invasion by novel taxa from outside the metacommunity, and (3) relied on a broader dispersal kernel. The most plausible MSS simulations demonstrate that interactions between a moderately strong environmental filter and dispersal dynamics in a complex landscape can create similar emergent biodiversity patterns to NM simulations and provide a plausible alternative hypothesis about the metacommunity dynamics that maintain diatom biodiversity in the MDV.
A major difference between NM and MSS simulations was how mat identity (orange or black, Figure 4A) was linked to immigration and recruitment dynamics at each patch. In both NM and MSS simulations, the magnitude of m was defined by mat identity. However, in MSS simulations, mat type also influenced the environmental filter [λ(E), Eq. 3] because the patch specific values for E were a function of both mat identity and chlorophyll-a standing biomass (x-axis in Figure 4C). While λ(E) had no effect on the magnitude of the immigrant flux, it did control which taxa were more likely to be successfully recruited at a given patch based on the local environment (E) and the habitat preferences (µ) and niche breadths (σ) of the resident and immigrant taxa in the recruitment pool.
With the addition of an environmental filter affecting species sorting, MSS simulations changed the nature of dispersal dynamics required to maintain biodiversity at both local and regional levels relative to NM simulations (Figure 8). In contrast with the NM simulations, scenarios with moderate environmental filtering do not support the predictions that black mats are more insular (e.g., more likely to show priority effects) or that orange mats are more likely to be sink habitats for dispersing propagules. Further, simulation outcomes demonstrated that metacommunities with moderate environmental filtering required flatter dispersal kernels (e.g., curves for small values of w in Figure 4B) to maintain biodiversity. Thus, the most plausible MSS scenarios require that propagules in the immigrant pool are well dispersed among the patches in the landscape.
MSS metacommunities could include dynamics that resemble the mass effects (ME) metacommunity paradigm. Under ME dynamics, environmental filters determine the composition of source communities in productive habitats, but emigration from source communities into nearby sink communities can overwhelm the influence of environmental filtering at those less productive habitats (Leibold et al., 2004). Thus, dispersal can moderate the spatial extent of the influence of environmental filters under ME dynamics. While both MSS and ME dynamics involve an interaction between dispersal and species sorting to determine community assembly outcomes, the ME paradigm represents a specific case of MSS dynamics, and not all MSS metacommunity scenarios would necessarily conform to the ME paradigm.

Testable Predictions Derived From Simulation Outcomes
The most plausible NM and MSS simulations identified in this study predict quite different dispersal dynamics for diatoms in the MDV. The NM simulations predict diatom assemblages in orange mats act as sinks in the metacommunity, and priority effects are stronger in black mats. These predictions could be tested by comparing field observations of diatom assemblages in the different mat types with observations from potential sources of immigrants. Two likely vectors for diatom dispersal in this system are aeolian transport Diaz et al., 2018) and drift in the water column in streams (Cullis et al., 2014;Kohler et al., 2018). If NM dynamics govern the structure of the diatom metacommunity, and wind/saltation is the primary driver of diatom dispersal, then these modeling outcomes suggest the composition of the pool of dispersing diatoms in both aeolian and drift samples will more closely represent the extant community in nearby, downwind or downstream orange mats than black mats.
Alternatively, if local environmental filters are influencing species sorting, the most plausible MSS simulation outcomes suggest that some specific dispersal dynamics are required to maintain the observed levels of alpha, beta, and gamma diversity. Specifically, Figure 8 shows that as filter strength increases (from NM to MSS to SSS), the dispersal kernel slope (w) must decrease (i.e., mean dispersal distances increase). Immigration rates (m) can still be low, but propagules must be capable of dispersing broadly for species sorting dynamics to maintain the observed biodiversity patterns. However, species sorting alone cannot maintain local diversity when environmental filters are too strong (as the SSS scenarios demonstrate). MSS model dynamics predict that the diatom assemblages in either aeolian samplers or drift nets distributed across Fryxell Basin would have very low levels of beta diversity compared to samples collected for orange and black mats. However, the flat dispersal kernel and larger values for ν identified in the most plausible MSS simulations suggest aeolian transport is a more likely vector than drift. Because of the influence of environmental filtering in MSS simulations, the MSS hypothesis predicts orange and black mat beta diversity should correlate somewhat with environmental variables, whereas beta diversity from samples of the dispersing pool should not correlate with environmental variables (or that correlation should be much weaker than currently observed in the extant communities). Lastly, MSS simulations suggest orange and black mat communities are equally insular, or orange mats are slightly more insular than black mats. Thus, we would not expect dispersing diatoms in either upstream drift or aeolian samples to be more similar to nearby orange mats than nearby black mats.

Moving Beyond the Model Assumptions
It is important to note the limitations on how these simulation outcomes can be interpreted. While the most plausible metacommunity models were able to maintain alpha, beta, and gamma diversity, we did not assess how simulated trends in community composition compared to the long-term record maintained by the MCM LTER. Past empirical studies have provided quite strong evidence that variation in stream hydrologic regimes serves as an environmental filter (Poff and Ward, 1990) to control diatom community composition in Fryxell Basin (Stanish et al., 2011(Stanish et al., , 2012. Thus, a next step is to refine our modeling approach to test predictions about community composition. However, identifying metacommunity dynamics that can maintain biodiversity was an important first step. Importantly, our simulation outcomes indicate that species sorting alone does not maintain biodiversity in this ecosystem that has been shown to be relatively stable. Our results may enhance previous findings (Stanish et al., 2011(Stanish et al., , 2012 by pointing to specific mechanisms by which hydrology might exert control over diatom community composition in MDV streams. We implemented filtering in MCSim using chlorophyll-a concentrations as a proxy for local stream harshness, based on the empirical relationships described in Stanish et al. (2011Stanish et al. ( , 2012, Kohler et al. (2015), and used empirically observed diatom occurrence data in the OMI analysis ( Figure 4C) to determine the habitat preferences for each taxon in the environmental filtering process (Eq. 3). However, a large proportion, but not all of the taxa in the regional species pool are aerophilic (Esposito et al., 2008;Spaulding et al., 2010;Stanish et al., 2011). In the current study, we did not link functional traits to dispersal-taxa were neutral with respect to dispersal ability-though the simulation outcomes already suggest dispersal dynamics are important. A crucial next step in modeling metacommunity dynamics in this system is to use diatom functional traits to provide variation in dispersal abilities in the next generation of metacommunity models. It is likely that simulations that incorporate dispersal traits will better predict composition than the models used in the current study, but we expect that such simulations will involve dispersal dynamics similar to those identified in the present study.
Lastly, the simulation scenarios used in this study did not include temporal variability in the landscape. We attempted to indirectly include the influence of spatial variation in hydrologic regimes in the landscape by using spatial variation in chlorophyll-a standing biomass as a proxy, but there is obviously an important temporal component to this that interacts with the recruitment process that was not captured in the metacommunity dynamics that were modeled in these simulation scenarios. We suspect that adding realistic temporal variation to the distribution of mat types in addition to diatom dispersal traits will be key to more accurately modeling diatom metacommunity dynamics.

CONCLUSION
Overall, the results presented in this study point to the importance of dispersal dynamics for maintaining both patchscale and metacommunity-scale measures of biodiversity in the diatom metacommunity in Fryxell Basin in the MDV. Both remote imagery and the MCM LTER record provided crucial information for characterizing the patchy mosaic of cyanobacterial mats in which the diatom assemblages reside. Results from this study suggest that shifts in the distribution of different mat types in Fryxell Basin, which can potentially be tracked using remote imagery (Salvatore, 2015;Salvatore et al., 2020), can be consequential for both local and metacommunity scale biodiversity. Importantly, this study highlights the necessity of an accurate understanding of how dispersal affects the scaling of community assembly processes, which will be important for predicting how diatom biodiversity will change in future climate scenarios (Shindell and Schmidt, 2004;Gooseff et al., 2017). Such an approach may also prove useful for understanding how dispersal affects other microbial metacommunities (Martiny et al., 2006;Heino et al., 2010;Soininen, 2012), both in the MDV, and beyond.

DATA AVAILABILITY STATEMENT
Code to run the simulations used in this study, input data used to initialize the simulations, and simulation results are available as an archived data package hosted by the Environmental Data Initiative (EDI) at: https://portal. edirepository.org/nis/mapbrowse?scope=edi&identifier=597 (doi: 10.6073/pasta/928879746d04e847c1b5acfa31f1cbfa). The code and input data are also available on GitHub at: https: //github.com/sokole/diatom_metacommunity_simulations.

AUTHOR CONTRIBUTIONS
ES designed the study, created the MCSim modeling platform, conducted the analyses, and wrote the manuscript. ES, JB, MS, and LS secured funding to support the analysis of remote imagery used in the study. MS analyzed the remote imagery and created the visible and NDVI maps used in this study. DM and JB are PIs with the MCM LTER program, which provided the algal mat and diatom observational data used in the study. DM, TK, and LS designed and executed the collection of the MCM LTER data used in this study. All authors contributed to interpretation of data and writing the manuscript.

FUNDING
This work was funded by United States National Science Foundation OPP (Antarctic Organisms and Ecosystems) awards 1744785, 1744849, and 1745053 and NSF OPP award 1637708 supporting the MCM LTER Project.