Biogeography at the limits of life: Do extremophilic microbial communities show biogeographical regionalization?

Aim: Biogeographical regions are the fundamental geographical units for grouping Earth’s biodiversity. Biogeographical regionalization has been demonstrated for many higher taxa, such as terrestrial plants and vertebrates, but not in microbial communities. Therefore, we sought to test empirically whether microbial communities, or taxa, show patterns consistent with biogeographical regionalization.


| I NTR OD U CTI ON
The classification of Earth's biota into biogeographical regions separated by dispersal barriers has captivated ecologists for centuries (Sclater, 1858;Wallace, 1876). The concept of biogeographical regionalization has yielded insight into the origins of biodiversity and areas of endemism (Lamoreux et al., 2006), informed us of species' conservation status (Buckley & Jetz, 2007) and revealed historical connectivity between communities (Cowen, Paris, & Srinivasan, 2006). However, early attempts to define these regions have been superseded by more quantitative methods, improving the robustness and reproducibility of region delineations (Dapporto, Ciolli, Dennis, Fox, & Shreeve, 2015;Kreft & Jetz, 2010;Vilhena & Antonelli, 2015). Coupled with these new methods, the ever-increasing availability of species distribution data has renewed interest in the concept of biogeographical regionalization. As a result, a far greater range of taxa have been studied than ever before in order to define biogeographical regions (Holt et al., 2013). Yet, our knowledge about how Earth's biota may be divided into biogeographical regions is still overwhelmingly based on multicellular (and usually large) eukaryotes. Many inconspicuous, but functionally critical organisms, such as microorganisms, remain poorly studied.
Consequently, it is unknown whether microbial communities may be grouped into biogeographical regions, similar to those observed for higher taxa.
Microorganisms are arguably the most functionally diverse and important organisms on Earth (Dinsdale et al., 2008;Fierer et al., 2012), driving every biogeochemical cycle (Falkowski, Fenchel, & Delong, 2008;Zak, Holmes, White, Peacock, & Tilman, 2003). Originally, microorganisms were assumed to have cosmopolitan distributions, with their small size and high population densities making them effective passive dispersers (Baas Becking, 1934;Finlay, 2002). From this assumption, it follows that biogeographical regionalization may not be possible because dispersal limitation is required for areas of endemism to occur (Ficetola, Mazel, & Thuiller, 2017) and produce regions with distinct communities. In contrast to cosmopolitanism, many recent studies have documented relationships between community turnover (the replacement of species) and geographical distance, indicative of dispersal limitation (e.g., Dumbrell, Nelson, Helgason, Dytham, & Fitter, 2010;Lear, Bellamy, Case, Lee, & Buckley, 2014), hinting that biogeographical regionalization of microbial communities could be possible. However, whilst the composition of microbial communities has been shown to differ over biogeographical regional scales, a formal test of whether microbial communities exhibit biogeographical regionalization is lacking.
In order to test for the presence of biogeographical regionalization in microbial communities, an ideal model community should have relatively low diversity, inhabit isolated environments and show a-priori evidence of dispersal limitation. The halite-associated Archaea fulfil these criteria. These Archaea typically belong to the class Halobacteria (more commonly referred to as haloarchaea) and are a major component of halite endolith communities (Henriet, Fourmentin, Delinc e, & Mahillon, 2014). Their entombment into the brine inclusions of halite crystals is believed to be an escape mechanism from desiccation and the increasingly chaotropic conditions present in evaporating brines (Hallsworth et al., 2007). Within these pockets they are able to survive over geological time scales (Gramain, Díaz, Demergasso, Lowenstein, & McGenity, 2011;McGenity, Gemmell, Grant, & Stan-Lotter, 2000). As with many extremophilic microbial communities, the halite-associated Archaea are typically less diverse than other microbial systems, facilitating more exhaustive sampling of the total diversity and improving detection of the less abundant endemic taxa, which are indicative of biogeographical regions. Furthermore, these Archaea occupy isolated 'habitat islands' that are physicochemically distinct from the surrounding environment. Many haloarchaea are obligately halophilic and lyse in less saline conditions (Oren, 1994), such as seawater, rendering the surrounding environment a physiological dispersal barrier. Finally, halite crystals form in highly similar conditions worldwide (i.e., saturated NaCl), thus ensuring that species filtering by the environment is low compared with many other environments. Any physicochemical differences between halite crystals (e.g., caused by underlying geology or climate) should themselves be spatially autocorrelated, meaning that species filtering by the environment should enhance, rather than obscure, biogeographical clustering. Such systems are therefore ideal for studies of community turnover and biogeography (Santos, Field, & Ricklefs, 2016). Previous studies of halophilic microbial communities have found evidence of community turnover at regional scales (Pagaling et al., 2009;Zhaxybayeva, Stepanauskas, Mohan, & Papke, 2013), suggesting the potential for biogeographical regions to form.
Overall, these properties render the halite-associated Archaea an ideal system in which to test for biogeographical regionalization of microbial communities.
Therefore, we examine the regional turnover (replacement of species over biogeographical regional scales) of halite-associated archaeal communities to test whether communities group together in a manner consistent with biogeographical regionalization. Using high-throughput Illumina HiSeq amplicon sequencing, we characterize the archaeal communities of halite from 17 locations, spanning three geographical regions. We apply robust biogeographical clustering methods to determine the extent to which archaeal communities, and taxa, show spatial patterns consistent with biogeographical regionalization. We propose the following three hypotheses: 1. (a) There will be a significant relationship between community turnover and geographical distance, and (b) the rate of community turnover will be greater at biogeographical regional scales than at within-region scales.
2. Communities will form biogeographical clusters that are statistically well supported and spatially coherent.
3. The presence and abundance of some archaeal taxa can predict the (bio)geographical origin of each community.

| M A TE RI A LS AN D M ET HOD S
We obtained 27 halite samples (in triplicate) from 17 locations in the years between 2006 and 2013 ( Figure 1 and Supporting Information Appendix S1). A photographic record of the samples and further details can be found in Supporting Information Appendix S2. We recorded the grain size, which reflects the time taken for the crystals to form, and the impurity colour, which provides a qualitative measure of the types of impurities and physicochemical environment present within the crystal (Sonnenfeld, 1995). Samples were stored in the dark at room temperature.

| Bioinformatic analyses
Owing to the length of the amplicon, it was not possible to pair-end align the forward and reverse sequences; therefore, all analyses were based on forward sequences only. This approach has been shown to have little effect on ecological conclusions (Werner, Zhou, Caporaso, Knight, & Angenent, 2012), and in our case, the forward sequence spans the V3 region of the 16S rRNA gene, which has been shown to perform well for profiling archaeal communities (Yu, García-Gonz alez, Schanbacher, & Morrison, 2008). Sequences were processed according to guidelines outlined by Dumbrell, Ferguson, and Clark (2016). Briefly, we quality trimmed sequences using Sickle (Joshi & Fass, 2011) at a threshold of Q20, trimming only the 3ʹ end of the sequence and discarding sequences with ambiguous nucleotides. Quality-trimmed sequences were error corrected using the BayesHammer algorithm implemented in SPAdes version 3.10.1, with default parameters (Bankevich et al., 2012;Nikolenko, Korobeynikov, & Alekseyev, 2013).
We removed primer sequences, calculated library sizes for each sample, and discarded sequences < 230 nucleotides in length using Linux shell We used VSEARCH (Rognes, Flouri, Nichols, Quince, & Mah e, 2016) to cluster sequences into operational taxonomic units (OTUs). First, sequences were de-replicated and singleton sequences discarded, as they are more likely to be artefacts (Flynn, Brown, Chain, MacIsaac, & Cristescu, 2015). We then clustered sequences into OTUs at 97 and 99% sequence similarity (referred to as 97% dataset and 99% dataset).
The 97% similarity threshold is the most frequently used, corresponding approximately to intragenus-level similarity (Yarza et al., 2014). The 99% threshold approximates to species-level similarity. We screened OTUs for chimeras against the RDP database (Cole et al., 2009) using VSEARCH, and discarded putative chimeras.

| Statistical analyses
We rarefied OTU tables to the smallest library size in each dataset (97% dataset, 27,554 sequences; 99% dataset, 26,578). We checked whether sample age was influencing the OTU richness or community composition using a negative binomial generalised linear model (GLM) and permutation-based multivariate analysis of variance (PERMA-NOVA), respectively.
In order to address our first hypothesis, we quantified community turnover using the b sim index, which purely quantifies community turnover, the process relevant to biogeographical regionalization (Baselga, 2010), and not nestedness, whereby communities are subsets of each other. Geographical distances between sampled communities were calculated as geodesic distances (Hijmans, 2016). We then tested for correlation between community turnover and geographical distance using Mantel tests, with Spearman's correlation coefficient and 10,000 permutations. We fitted piece-wise regressions to determine breakpoints in the relationship, showing the geographical distance at which the slope of the relationship changes (Castro-Insua, G omez-Rodríguez, & Baselga, 2016).
To investigate our second hypothesis, we adopted a clustering approach as described by Kreft and Jetz (2010). Briefly, this approach involves clustering communities based on the b sim turnover matrix, creating a dendrogram. This dendrogram can be split into k clusters representing bioregions. The quality and biological interpretability of the resulting clusters are then checked via statistical metrics and mapping. Biogeographical regionalization may be inferred when clustering solutions are both statistically robust and spatially coherent.
To cluster communities, we used three different clustering algorithms to ensure our conclusions were robust. The unweighted pairgroup method using arithmetic averages (UPGMA) defines the distance between clusters as the average distance between all the communities within each cluster. Kreft and Jetz (2010) found that UPGMA best preserved information present in the original distance matrix. Dapporto, Ciolli et al. (2015) also compared clustering algorithms on datasets of varying completeness. They found that for less well-sampled datasets, the Ward method clustered communities most accurately, whereas for intensely sampled datasets, PAM produced the most accurate clusters. To cluster the communities, we used the methodology described by Dapporto, Ramazzotti et al. (2015). This approach overcomes the biases introduced by having zero similarity or tied values in the dissimilarity matrix (Bloomfield, Knerr, & Encinas-Viso, 2017) (Kaufman & Rousseeuw, 1990). Our second metric, 'explained dissimilarity' (Holt et al., 2013), is a ratio of sums of mean dissimilarity within regions to total dissimilarity over the entire dissimilarity matrix. Explained dissimilarity tends towards one as k tends towards the number of communities. We follow the approach of Holt et al. (2013), who indicated that a threshold of .9 provides sufficient support to infer regionalization.
However, we also examined the cluster solution that produced the greatest incremental increase in explained dissimilarity, which we refer to as the 'knee solution', because this has been proposed to be a more suitable indicator of optimal cluster number (Kreft & Jetz, 2013). After identifying statistically supported clustering solutions, we determined the spatial coherence of clusters by mapping them. To check whether the measured physicochemical parameters (grain size or impurity) explained any clustering patterns observed, we used PERMANOVA.
We included location as the first variable in the model to account for confounding spatial effects. Statistical significance of physicochemical variables was then assessed based on the 'marginal' effects (e.g., after controlling for spatial location), with 999 random permutations. We conducted non-metric multidimensional scaling (NMDS) analysis as a means of visualizing these results.
To test our third hypothesis, we investigated whether the relative abundance of halite-associated archaeal genera could predict the biogeographical origin of a given community using the machine learning method of random forest classification (RFC). Random forests provide an effective method for classification in ecology (Cutler et al., 2007) and are built from an ensemble of classification trees, in which observations of the dependent variable form the leaves and independent variables form the branches. Each tree is trained on a subset of observations and independent variables, and the overall classifier is built by combining predictions from these trees to obtain a more robust classification. We summed the abundances of all OTUs identified to the genus level and converted these abundances to relative abundances. OTUs not identified to genus were excluded from this analysis.
We classified communities (see Supporting Information Appendix S1) based on their biogeographical region (classes: Palearctic, Saharo-Arabian, Madagascan as defined by Holt et al., 2013), geographical region (classes: eastern Europe, western Europe, Mediterranean or west African) and nearest ocean (classes: Atlantic or Indian). We initialized 10,000 trees, and each tree was trained on six archaeal genera.
We normalized the sample size from each class to the size of the smallest class to minimize the effects of class size imbalance (e.g., more observations of European communities than African communities).
Additionally, for the biogeographical and geographical classifiers, we dropped excessively small classes (Saharo-Arabian; n 5 4 and west African; n 5 3), to reduce the imbalance between classes further. We evaluated the overall accuracy of each classifier using the out-of-bag error rate, which quantifies the classifier's ability to classify a given community correctly when it is excluded from the training set. We determined which archaeal genera were the best predictors of biogeographical origin by quantifying variable importance, using the mean decrease in accuracy (MDA) and mean decrease in Gini index (MDGI).
The MDA shows the change in accuracy of the classifier with and without a given variable. Important variables will result in a large decrease in accuracy when they are excluded from the classifier, resulting in large MDA values. The MDGI shows the purity of the groups created when the classifier splits the dataset using a given predictor. A good predictor will create homogeneous groups, in which all data points belong to the same class, resulting in a large decrease in MDGI. We also examined partial dependence plots (Hastie, Tibshirani, & Friedman, 2009). In the context of our study, these plots show how the probability of a community being classified into a given biogeographical region changes in relationship to the relative abundance of a given archaeal genus.
All analyses were conducted in R (R Developement Core Team, 2016), using the vegan (Oksanen et al., 2015), recluster (Dapporto, Ramazzotti et al., 2015) and randomForest (Liaw & Wiener, 2002) packages. However, piece-wise regressions between geographical distance and community turnover suggested that this correlation was largely driven by high turnover at small spatial scales (Figure 2). For both (97 and 99%) datasets, a steep positive relationship was found at smaller spatial scales, with breakpoints estimated at 420.5 km (standard error 5 46.9 km) and 334.6 km (standard error 5 23.7 km), respectively.

| RE S U LT S
After these breakpoints, community turnover was independent of geographical distance (Figure 2). Davies tests confirmed that the prebreakpoint slope was significantly greater than the post-breakpoint

| Do microbial communities cluster into biogeographical regions?
We determined whether archaeal communities group into biogeographical regions by applying three different clustering algorithms (UPGMA, Ward and PAM). To assess the degree of biogeographical clustering within these communities, we first determined the appropriate number of clusters (k) into which our communities should be grouped by examining the cluster quality (using mean silhouette width and explained dissimilarity) for values of k from 2 to 16. For the 97% dataset, statistical support for cluster solutions was poor, as the mean silhouette width never exceeded .25 for any value of k (Figure 3a). In contrast, for the 99% dataset, all three clustering algorithms exceeded .25 for values of k > 12, showing that reasonable statistical support was gained when communities were grouped into > 12 regions. All three clustering algorithms yielded similar results when assessed by the explained dissimilarity metric (Figure 3b). Explained dissimilarity values > .9 were considered to provide good support for a given cluster solution. To satisfy this threshold, communities were grouped into 8-10 (97% dataset) or 9-12 (99% dataset) clusters, depending on the cluster algorithm used. For both 97 and 99% datasets, the Ward algorithm required the fewest clusters to reach this threshold, and UPGMA the most. We also identified the number of clusters (k) that resulted in the greatest increase in explained dissimilarity (knee solutions). For the 97% dataset, this occurred when communities were clustered into three (PAM and Ward) or four (UPGMA) clusters, whereas for the 99% dataset, the greatest increase in explained dissimilarity was found when communities grouped into three (UPGMA and Ward) or four (PAM) clusters.
We examined the spatial coherence of cluster solutions for the minimal number of clusters (k) required to exceed the explained dissimilarity threshold of .9, as well as solutions that yielded the greatest increase in explained dissimilarity. For both 97 and 99% datasets and all three clustering algorithms, mapping revealed poor spatial coherence showed that archaeal communities clustered only weakly by impurity, but not by grain size (Figure 5).

| Can certain haloarchaeal genera be used as indicators of a community's biogeographical origin?
We tested whether the abundance of certain archaeal genera could predict any of three classifiers (biogeographical region, geographical region and nearest ocean) of a community, using random forest classifi- Mediterranean (class error 5 7.7%) communities than east African communities (class error 5 10.5%).
To determine which archaeal genera were the best predictors of a community's oceanic, biogeographical or geographical origin, we quantified the importance of each variable (genus) to each RFC (Supporting Information Appendix S4). Haloquadratum was the best genus for classifying geographical region, followed by Halapricum and Halobaculum.
Partial dependence plots revealed that, as the relative abundance of

| DI SCUS SION
We studied halite-associated Archaea to determine whether archaeal communities can be clustered into biogeographical regions comparable to those observed for most higher organisms. Our results show that, despite community turnover correlating with geographical distance over small spatial scales (< 500 km), communities do not cluster into spatially coherent biogeographical regions. We found little statistical support for clustering communities into few (two or three) biogeographical regions, which would be the number of regions expected for higher organisms, such as terrestrial vertebrates (Holt et al., 2013) or plants (Takhtajan, 1986). Furthermore, when we clustered communities into a greater number of regions, the spatial configuration of these regions was not consistent with biogeographical regionalization. Lastly, we demonstrated that although communities may not show the expected biogeographical patterns, some individual genera do, as their abundances were found to be good predictors of the biogeographical origin of the community.
Numerous studies have demonstrated that microbial communities differ over continental to regional scales (Lauber, Hamady, Knight, & Fierer, 2009;Papke, Ramsing, Bateson, & Ward, 2003;Whitaker, Grogan, & Taylor, 2003), including studies on halophilic microbes (Pagaling et al., 2009;Zhaxybayeva et al., 2013). However, to our knowledge, no studies have tested quantitatively whether such differences are consistent with the concept of biogeographical regionalization, thus it remains unknown whether the processes that shape microbial communities are capable of forming biogeographical patterns over the spatial scales relevant to other organisms. Glassman et al. (2015) examined fungal spore banks of soils across North America, showing that community turnover was significantly related to geographical distance and, using ordination techniques, that fungal communities appeared to group in a regional manner. Consistent with our study, they found that the highest rate of community turnover occurred over subregional scales, as evidenced by their Mantel correlogram, which shows change from positive to negative correlation over spatial scales of c. 500 km. Initially, this might indicate that microbial biogeographical regions are smaller than those defined for higher taxa and more comparable to subregions. However, in our study, this idea is poorly supported by the fact that even for larger values of k (indicating more and smaller regions), the spatial coherence of these clusters was poor.
A global study of soil fungi (Tedersoo et al., 2014) revealed communities that did not cluster in a spatially coherent manner, which is in contrast to the findings of Glassman et al. (2015) and in agreement with our results. For instance, fungal communities of Europe clustered with those of North America, and those of Oceania clustered with South America. Furthermore, a study of the bacterial communities on Tamarix spp. leaf surfaces showed that communities clustered in a manner at odds with their spatial configuration (Finkel et al., 2012). Specifically, communities from around the Dead Sea (Middle East) clustered more closely with those from the Sonoran Desert (North America) than Mediterranean communities. Combined with our results, these studies provide further evidence that biogeographical regionalization may be unlikely in microbial communities.
One possible reason for no evidence of biogeographical regionalization in these communities is that some halophilic Archaea may be dif-  et al., 2015), which may help them to spread between habitat islands. Furthermore, wind-or human-mediated dispersal of halite crystals may disperse entombed haloarchaea. Wind is known to play a role in dispersing free-living microbes over continental distances (Favet et al., 2013;Kellogg & Griffin, 2006) and is likely to disperse small halite crystals, along with endolithic microbes, over such distances. Human transport of salt as a commercial product and as a de-icing agent on roads may also aid the dispersal of halite endolithic communities. However, such dispersal would select for those Archaea capable of survival in halite crystals, filtering out some taxa, as evidenced by the disparity between brine and halite crystal archaeal communities described previously (Henriet et al., 2014). Finally, dispersal via seawater could be possible for some haloarchaeal taxa, because viable cells have been isolated from seawater and coastal sediments (Purdy et al., 2004;Rodriguez-Valera, Ruiz-Berraquero, & Ramos-Cormenzana, 1979).
Ancient halite deposits can become exposed in deep water horizons, where they may dissolve, creating stratified deep-sea brines, which are a potential source of extremely halophilic Archaea (Antunes, Ngugi, & Stingl, 2011). However, although short-term (c. 24 hr) or partial survival at seawater salinity has been found in a number of haloarchaea (Torreblanca et al., 1986), the majority of genera detected in this study, particularly the most abundantly detected genera, are known exclusively from hypersaline habitats, and their cells lyse at seawater salinity. Therefore, seawater is an unlikely medium for their dispersal.
Furthermore, the deposition of cells from ancient halite into modern hypersaline environments would be most likely to occur over regional extents (e.g., owing to oceanic currents), thus increasing the compositional similarity of sites within a region. Finally, even with connectivity between ancient and modern halite, there is no guarantee that those cells will become established and multiply (Jones, Ramoneda, Rivett, & Bell, 2017). Therefore, the influence of ancient haloarchaea on the clustering patterns observed here should be minimal. Even so, the degree to which other potential dispersal vectors contribute to connectivity between sites is unknown and warrants further research, as connectivity between contemporary halite deposits may be a better measure of isolation for these communities than geographical distance alone.
An alternative explanation as to why biogeographical clustering was not observed in these archaeal communities is that environmental filtering, because of physicochemical differences between the halite crystals, could obscure biogeographical clustering. Within hypersaline systems, salinity (concentration of sodium chloride; NaCl) has been shown to be the predominant physicochemical variable causing environmental filtering of microbial communities (Baati et al., 2008;Benlloch et al., 2002;Casamayor et al., 2002;Herlemann et al., 2011).
However, the role of physicochemical differences in structuring microbial communities between hypersaline habitats is less well known, Overall, we found little evidence to support the existence of biogeographical regions in communities of extremely halophilic Archaea.
We demonstrated that, despite finding evidence of a distance-decay relationship in these communities, clustering them into regions did not produce spatially coherent regions. We suggest that the cause of this may be long-distance dispersal of some haloarchaeal taxa, as we identified three particularly abundant and widespread species that were universally detected across all samples. However, certain individual taxa are able to indicate a given community's biogeographical origins accurately, suggesting highly differential dispersal abilities in haloarchaea. Taken together, our results suggest that geographical distance alone may be a poor indicator of isolation in microbial communities and that more work is needed to examine the role of connectivity in microbial biogeography.

ACKNOWLEDGMENTS
We are grateful to Dr Andrew Crombie (University of East Anglia) for providing the Portuguese samples and to Dr Pascal Danthu analysed the data. All authors contributed to the writing of the manuscript.

DATA ACCESSIBILITY
Sequence data from this study are deposited in the European Nucleotide Archive under the project accession number: PRJEB19885. All other data needed to reproduce this study can be found in Supporting Information Appendix S1 (Supplementary file 1).