Isaac Overcast
email: isaac.overcast@maine.edu

Víctor Noguerales

Emmanouil Meramveliotakis

Carmelo Andújar

Paula Arribas

Thomas J Creedy

Brent C Emerson

Alfried P Vogler

Anna Papadopoulou

Hélène Morlon

Inferring the ecological and evolutionary determinants of community genetic diversity

Keywords: Modeling community biodiversity genetics Biodiversity genomics, community ecology, metabarcoding, population genetics, machine learning, eco-evolutionary modeling

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Biodiversity in ecological communities is structured hierarchically across spatial and temporal scales. Many open questions remain as to how this structure accumulates. For example, what are the relative contributions of dispersal versus in situ speciation? Or, how important are stochastic drift versus deterministic processes? Until recently, these questions have been investigated by isolated disciplines (e.g. macroecology, comparative phylogeography, macroevolution) using tools and data that tend to focus on only one axis of community-scale data (e.g. phylogenies, relative abundances, and/or trait information). Yet we know that there are feedbacks among processes that respond on short, medium, and long time scales: local changes of abundance, accumulation of population genetic variation, and speciation processes, respectively [START_REF] Johnson | An emerging synthesis between community ecology and evolutionary biology[END_REF]. For example, changes in population size can affect population genetic variation, which can in turn promote or impede speciation. Only recently have disciplines arisen to simultaneously evaluate community-scale data from multiple levels of biological organization. Community phylogenetics, community genetics [START_REF] Barker | Linking plant genes to insect communities: Identifying the genetic bases of plant traits and community composition[END_REF], and comparative phylogeography are mature disciplines which demonstrate the power of merging timescales. Theory and models of community-scale abundances and phylogenies are well developed, yet only recently have models of community-scale genetic diversity been developed [START_REF] Overcast | An integrated model of population genetics and community ecology[END_REF][START_REF] Overcast | A unified model of species abundance, genetic diversity, and functional diversity reveals the mechanisms structuring ecological communities[END_REF]. Genetic diversity of species within ecological communities can provide a unifying view of the ecological and evolutionary processes that structure and maintain biodiversity, so continued development of such models is critical.

Using traditional methods, obtaining genetic samples from multiple individuals from all species within a community requires a massive amount of effort in field sampling, laboratory work and sequencing. Whole organism community DNA metabarcoding [START_REF] Andújar | Why the COI barcode should be the community DNA metabarcode for the metazoa[END_REF][START_REF] Creedy | Coming of age for COI metabarcoding of whole organism community DNA: Towards bioinformatic harmonisation[END_REF][START_REF] Yu | Biodiversity soup: metabarcoding of arthropods for rapid biodiversity assessment and biomonitoring[END_REF], eDNA metabarcoding [START_REF] Cordier | Ecosystems monitoring powered by environmental genomics: a review of current strategies with an implementation roadmap[END_REF][START_REF] Deiner | Environmental DNA metabarcoding: Transforming how we survey animal and plant communities[END_REF][START_REF] Taberlet | Environmental DNA[END_REF], and mitochondrial mitogenomics [START_REF] Crampton-Platt | Mitochondrial metagenomics: letting the genes out of the bottle[END_REF][START_REF] Gómez-Rodríguez | Validating the power of mitochondrial metagenomics for community ecology and phylogenetics of complex assemblages[END_REF], 2017) have significantly changed this equation, allowing community-scale sequence data to be obtained at scale, through space and time [START_REF] Davies | The founding charter of the Genomic Observatories Network[END_REF][START_REF] Marquina | Establishing arthropod community composition using metabarcoding: surprising inconsistencies between soil samples and preservative ethanol and homogenate from Malaise trap catches[END_REF].

Global scale biodiversity inventorying efforts like the Global Malaise Trap Program and BioScan [START_REF] Hobern | BIOSCAN-revealing eukaryote diversity, dynamics, and interactions[END_REF] along with repositories of spatially tagged sequence data like GeoME [START_REF] Deck | The Genomic Observatories Metadatabase (GeOMe): A new repository for field and sampling event metadata associated with genetic samples[END_REF] and massive, community curated reference barcode databases like BOLD [START_REF] Ratnasingham | bold: The Barcode of Life Data System[END_REF] provide unprecedented access to community-scale genetic data allowing to investigate patterns of genetic diversity at the global scale [START_REF] Miraldo | An Anthropocene map of genetic diversity[END_REF][START_REF] Pelletier | Geographical range size and latitude predict population genetic structure in a global survey[END_REF][START_REF] Theodoridis | Evolutionary history and past climate change shape the distribution of genetic diversity in terrestrial mammals[END_REF]. DNA barcoding as initially practiced did not consider intraspecific genetic variation as a quantity of interest, the original intention being species identification either by matching to a database or clustering by a sequence similarity threshold to discover operational taxonomic units (OTUs; [START_REF] Blaxter | Defining operational taxonomic units using DNA barcode data[END_REF]. While estimating abundances from metabarcoding data is a controversial topic [START_REF] Deagle | Counting with DNA in metabarcoding studies: How should we convert sequence reads to dietary data[END_REF][START_REF] Elbrecht | Can DNA-Based Ecosystem Assessments Quantify Species Abundance? Testing Primer Bias and Biomass--Sequence Relationships with an Innovative Metabarcoding Protocol[END_REF][START_REF] Piñol | The choice of universal primers and the characteristics of the species mixture determine when DNA metabarcoding can be quantitative[END_REF][START_REF] Zhou | Ultra-deep sequencing enables high-fidelity recovery of biodiversity for bulk arthropod samples without PCR amplification[END_REF], recent work has demonstrated metabarcoding datasets can be used for estimating intraspecific genetic diversity at the scale of the whole community [START_REF] Adams | Beyond Biodiversity: Can Environmental DNA (eDNA) Cut It as a Population Genetics Tool?[END_REF][START_REF] Andújar | Validated removal of nuclear pseudogenes and sequencing artefacts from mitochondrial metabarcode data[END_REF][START_REF] Arribas | The limited spatial scale of dispersal in soil arthropods revealed with whole-community haplotype-level metabarcoding[END_REF][START_REF] Elbrecht | Estimating intraspecific genetic diversity from community DNA metabarcoding data[END_REF][START_REF] Sigsgaard | Population-level inferences from environmental DNA-Current status and future perspectives[END_REF][START_REF] Turon | From metabarcoding to metaphylogeography: separating the wheat from the chaff[END_REF].

Model based methods for inferring processes from patterns in genetic data of multiple codistributed species have been in development for decades [START_REF] Arbogast | Comparative phylogeography as an integrative approach to historical biogeography[END_REF]), yet numerous limitations and opportunities remain to fully exploit the signal in the genetic data. For example, comparative phylogeographic methods leverage the power of aggregated population genetic inferences of demographic history from multiple co-distributed species to answer fundamental questions about processes underlying community diversification, assembly, and macroecology [START_REF] Edwards | The Evolution of Comparative Phylogeography: Putting the Geography (and More) into Comparative Population Genomics[END_REF][START_REF] Hickerson | Phylogeography's past, present, and future: 10 years after Avise[END_REF][START_REF] Papadopoulou | Toward a paradigm shift in comparative phylogeography driven by trait-based hypotheses[END_REF].

However, comparative phylogeographic methods (e.g. [START_REF] Huang | MTML-msBayes: approximate Bayesian comparative phylogeographic inference from multiple taxa and multiple loci with rate heterogeneity[END_REF][START_REF] Oaks | Full Bayesian Comparative Phylogeography from Genomic Data[END_REF] typically do not model ecological processes, making it difficult to connect results with specific predictions of species interactions and coexistence. These methods have tended to focus on general models of shared demographic histories [START_REF] Burbrink | Asynchronous demographic responses to Pleistocene climate change in Eastern Nearctic vertebrates[END_REF][START_REF] Satler | Do ecological communities disperse across biogeographic barriers as a unit?[END_REF][START_REF] Stone | Reconstructing community assembly in time and space reveals enemy escape in a Western Palearctic insect community[END_REF], rather than on models that are explicitly parameterized from ecological community assembly theory (but see [START_REF] Bunnefeld | Whole-genome data reveal the complex history of a diverse ecological community[END_REF], and they also do not account for processes above the species level (e.g. macroevolution). Bottom-up models of community assembly have begun to move beyond focusing on predictions of abundance distributions, to considering intraspecific genetic polymorphism within a dynamic non-equilibrium assembly framework [START_REF] Laroche | A neutral theory for interpreting correlations between species and genetic diversity in communities[END_REF][START_REF] Vellend | Drawing ecological inferences from coincident patterns of population-and community-level biodiversity[END_REF], allowing to investigate the correlation between species diversity and genetic diversity in ecological communities [START_REF] Lamy | The contribution of speciesgenetic diversity correlations to the understanding of community assembly rules[END_REF][START_REF] Papadopoulou | Testing the Species-Genetic Diversity Correlation in the Aegean Archipelago: Toward a Haplotype-Based Macroecology?[END_REF][START_REF] Vellend | Species diversity and genetic diversity: parallel processes and correlated patterns[END_REF]. Despite these important efforts, and with some key exceptions [START_REF] Aguilée | Clade diversification dynamics and the biotic and abiotic controls of speciation and extinction rates[END_REF][START_REF] Bunnefeld | Whole-genome data reveal the complex history of a diverse ecological community[END_REF][START_REF] Gascuel | How Ecology and Landscape Dynamics Shape Phylogenetic Trees[END_REF][START_REF] Overcast | An integrated model of population genetics and community ecology[END_REF][START_REF] Overcast | A unified model of species abundance, genetic diversity, and functional diversity reveals the mechanisms structuring ecological communities[END_REF], opportunities remain for developing a dynamic eco-evolutionary model of community assembly which makes predictions of community scale genetic data, as well as other relevant predictions such as abundances, trait values, and phylogenetic relatedness, at the community scale or larger.

Here, we present the Island Biodiversity Genomics (iBioGen) toolkit, a software package for modeling and testing hypotheses using multi-dimensional data including community-scale genetic diversity such as can be obtained from contemporary metabarcoding experiments. The iBioGen package generates predictions of multi-dimensional data to understand how the strength of various ecological and evolutionary processes affect community composition. Model simulations can be paired with a machine learning inference framework to test hypotheses about community assembly mechanisms and estimate key ecological and evolutionary parameters from empirical metabarcoding datasets. Methods within the iBioGen package are taxon and scale independent, equally applicable to any island community system, from vertebrate island assemblages to human microbiome datasets. In its current state of development and depending on the purpose of the analysis and the the data in hands, the iBioGen package can be used, for example, to differentiate neutral assembly processes from competition-driven and/or abiotic niche filtering, and to estimate the equilibrium state of a community as well as the strength of filtering effects, and colonization and speciation rates. XXXX. We illustrate the application of the iBioGen package on a recently published dataset of soil microarthropod metabarcode data from 44 communities sampled in the montane forests of Cyprus [START_REF] Noguerales | Community metabarcoding reveals the relative role of environmental filtering and spatial processes in metacommunity dynamics of soil microarthropods across a mosaic of montane forests[END_REF]. Using correlative statistical approaches, [START_REF] Noguerales | Community metabarcoding reveals the relative role of environmental filtering and spatial processes in metacommunity dynamics of soil microarthropods across a mosaic of montane forests[END_REF] found that mid-altitude communities were dominated by neutral processes while higher elevation communities were dominated by environmental filtering. This dataset can thus conveniently be used to validate our model-based approach and to illustrate its added value, such as providing measures of confidence around best supported models, as well as estimates of community size and colonization rates. We find that neutral communities are larger and characterized by higher colonization rates than communities structured by environmental filtering.

Materials and Methods

Model overview

Our model takes its roots in the theory of island biogeography (ref MacArthur & Wilson), in the spirit of the continent-island model, where the composition of an evolving metacommunity (sensu [START_REF] Hubbell | The Unified Neutral Theory of Biodiversity and Biogeography (MPB-32)[END_REF] influences the ecological assembly of a small, isolated, spatially-implicit local community which is connected to the metacommunity by immigration (Fig. 1). Under this framework, the metacommunity is a pool of species that can each experience speciation and go extinct. In our model, the metacommunity evolves according to a lineage-based birth/death process with small shifts in speciation rate at branching events (following [START_REF] Maliet | A model with many small shifts for estimating species-specific diversification rates[END_REF] and is characterized by the species phylogeny. Additionally, lineage abundance in the metacommunity is a function of the population growth rate, a lineage-dependent trait which evolves along the branches of the phylogeny by Brownian motion (BM). Local community assembly follows an individual-based birth/death/colonization/speciation process forward in time, in the tradition of the unified neutral theory of biodiversity (UNTB; [START_REF] Hubbell | The Unified Neutral Theory of Biodiversity and Biogeography (MPB-32)[END_REF] and its descendents (e.g. [START_REF] Overcast | A unified model of species abundance, genetic diversity, and functional diversity reveals the mechanisms structuring ecological communities[END_REF][START_REF] Rosindell | Protracted speciation revitalizes the neutral theory of biodiversity[END_REF][START_REF] Rosindell | A unified model of island biogeography sheds light on the zone of radiation[END_REF]. Following [START_REF] Overcast | A unified model of species abundance, genetic diversity, and functional diversity reveals the mechanisms structuring ecological communities[END_REF], genetic variation per species within the local community is modeled by a backward-time coalescent process, which accounts for historical demography (e.g. population size change histories) as well as recurrent migration from the metacommunity. Importantly, the timescales of the macroevolutionary and ecological processes within the metacommunity and local community are unified, thus the metacommunity continues to evolve for the duration of the local assembly process.

The metacommunity

We model diversification within the metacommunity using a modified birth-death process with small shifts in speciation rate, following the ClaDS model of [START_REF] Maliet | A model with many small shifts for estimating species-specific diversification rates[END_REF]. Our implementation of this model comprises three key parameters: the ancestral speciation rate (λ0), that is, the rate of speciation at the root of the tree; the trend in change of speciation rate at branching events (ClaDS α); and the stochastic variation in rate change about the trend (ClaDS σ; for details see [START_REF] Maliet | A model with many small shifts for estimating species-specific diversification rates[END_REF]. [START_REF] Maliet | A model with many small shifts for estimating species-specific diversification rates[END_REF] also defined a compound parameter (ClaDS m) expressed as α*exp(σ 2 /2), to quantify the mean relative change in speciation rate at branching events, irrespective of the parental rate. Beyond branch-specific speciation rates (λi), we extend the ClaDS model to further characterize for each lineage i its abundance at time t (Ai(t)), a population growth rate (ri(t)) and a morphological trait (xi(t); such as body size), the trait evolution dynamics of which are detailed below. Departing from the ClaDS model, in which extinction is an independent parameter (µi) which may be fixed to 0, homogenous and constant across the tree, or lineage specific with a constant turnover, in our model extinction is an emergent property of lineage-specific abundances and growth rates, whereby if the abundance of lineage i falls below some threshold value of population viability (by default 1), lineage i goes extinct.

We now consider the trait evolution model of the lineage-specific growth rates (ri) and morphological traits (xi), and the temporal dynamics of abundances (Ai). We begin with an ancestral population growth rate (r0) and trait value (x0) at the root of the tree, and global rates at which these traits evolve (σr and σx, respectively), all of which are configurable model parameters.

Growth rates and morphological traits per lineage evolve concurrently and independently along the branches of the phylogeny following a BM process as the macroevolutionary process unfolds.

In other words, the growth rate for lineage i after a given time interval (Δt) beginning with rate ri at time t with rate of change σr is given by ri(t + Δt) ~ N(ri(t), σr 2 Δt) with a similar formulation for morphological trait evolution. In our model, the trait values (xi) have no effect on lineage diversification within the metacommunity (but see [START_REF] Aristide | Understanding the effect of competition during evolutionary radiations: an integrated model of phenotypic and species diversification[END_REF][START_REF] Fitzjohn | Quantitative traits and diversification[END_REF] for how this might be implemented), though they will play an important role in the local community assembly process (see below). Lineage-specific abundances also fluctuate throughout the diversification process, and begin with an ancestral abundance (A0) assigned at the root. Lineage-specific abundances (Ai) grow or decline along the branches as a function of growth rate (ri) and time, following a continuous model of logistic growth, with abundance of lineage i after time interval Δt given by

𝐾 1 + ((𝐾-𝐴 ! (𝑡))/𝐴 ! (𝑡))𝑒𝑥𝑝(-𝑟 ! (𝑡)𝛥𝑡)
where K is the per lineage carrying capacity, which is fixed and constant, and Ai(t) and ri(t) are the abundance and population growth rate of lineage i at time t, respectively [START_REF] Murray | MathematicalBiology I. An Introduction[END_REF]. The global limit on lineage abundance (K) captures the fact that, while growth rate may vary among lineages, resource availability constrains the maximum sustainable population size. As mentioned above, lineages go extinct if their abundance crosses a lower-bound threshold for population viability. Finally, at each speciation event, in the spirit of allopatric speciation, ancestral abundance is split by random fission [START_REF] Etienne | The neutral theory of biodiversity with random fission speciation[END_REF][START_REF] Hubbell | The Unified Neutral Theory of Biodiversity and Biogeography (MPB-32)[END_REF] between the two descendent lineages, while ensuring that each offspring lineage retains a population size above the chosen threshold value of population viability.

The metacommunity evolves until a specified stopping condition is met, which conditions include either that: 1) a given amount of time has passed (as specified in millions of years); or 2) a given number of extant tips have been realized (but see [START_REF] Maliet | Fast and Accurate Estimation of Species-Specific Diversification Rates Using Data Augmentation[END_REF] for an unbiased approach to conditioning the phylogeny on the number of extant tips). At this point (t0 in Fig. 1), the local community is initialized (see below), and from this moment the local community and the metacommunity evolve in lockstep. This is an improvement over previous UNTB-derived models (e.g. [START_REF] Overcast | A unified model of species abundance, genetic diversity, and functional diversity reveals the mechanisms structuring ecological communities[END_REF] where the metacommunity is simulated first and remains static for the duration of the local community simulation. To insure consistency across scales (for example, a speciation event occurring in the local community contributes a new species in the metacommunity, and should thus be recorded in the metacommunity phylogeny), we establish a mechanism for maintaining synchronization. First, we homogenize timescales: we scale time in the metacommunity (measured in millions of years) to generations (the timescale of the local community, see below) by dividing by the average generation time of the local community taxa (in years per generation), with generation time being a model parameter which is typically fixed given a priori knowledge of the empirical system. Next, we define 'epochs', which are the successive times between speciation events (happening either in the metacommunity or in the local community). Within each successive epoch, the local community assembles. When a speciation event occurs, the metacommunity phylogeny is updated (explained above in the case of a speciation event in the metacommunity, and below in the case of a speciation event in the local community). This process continues until the local community stopping criterion is satisfied (see below), such that the number of epochs is determined by the time required for the local assembly process to reach this criterion and is not directly parameterized. At this point the simulation is complete, and simulation parameters and state variables are preserved.

The local community

We model the local community as an individual-based, spatially implicit system in the manner of the ecological neutral theory [START_REF] Hubbell | The Unified Neutral Theory of Biodiversity and Biogeography (MPB-32)[END_REF]. The basic propositions of the UNTB include that the local community is composed of a fixed and constant number of individuals (J), that it is a dispersal-limited sample of the metacommunity, and that all individuals of all species are ecologically equivalent (an assumption which will be relaxed in the following section). Individualbased models of community assembly inspired by the UNTB have tended to treat the metacommunity as a static pool of species, following some well-studied abundance distribution (e.g. log-series), which is very large and fixed with respect to the timescale of the local community (e.g. [START_REF] Overcast | An integrated model of population genetics and community ecology[END_REF][START_REF] Rosindell | A unified model of species immigration, extinction and abundance on islands[END_REF][START_REF] Rosindell | A unified model of island biogeography sheds light on the zone of radiation[END_REF], but see [START_REF] Davies | Neutral biodiversity theory can explain the imbalance of phylogenetic trees but not the tempo of their diversification[END_REF][START_REF] Jabot | Inferring the parameters of the neutral theory of biodiversity using phylogenetic information and implications for tropical forests[END_REF][START_REF] Manceau | Phylogenies support out-of-equilibrium models of biodiversity[END_REF]. Our metacommunity model makes dynamic predictions of phylogenetically informed abundances, motivated by the axiom that population growth rate, and therefore species abundance in the metacommunity, should correlate with life history traits that are phylogenetically conserved (generation time, fecundity, metabolic rate, etc.).

The local community of size J is initialized (at t0 in Fig. 1) by sampling J individuals following a multinomial distribution parameterized by the metacommunity species abundances at t0. The local community process then proceeds in time units of "generations". For a community of size J, a generation is composed of J discrete timesteps, with each timestep representing the death and replacement of one individual in the local community. Upon the death of an individual, the vacated ecological space is immediately occupied with probability 1 -m by the offspring of an individual randomly sampled from the local community, and with probability m it is replaced by a migrant from the metacommunity. At a given timestep, if migration from the metacommunity occurs, the migrating lineage is sampled following a multinomial distribution parameterized by the metacommunity species abundances. In other words, metacommunity species are sampled for migration with probability proportional to their metacommunity abundance. If individuals of the same species as the migrant are already present in the local community, this is considered as recurrent migration, and the count of post-colonization migrants for that species is incremented.

These migrant counts will be averaged over time to calculate a continuous rate of migration per species for the genetic simulations (see below). This process of individual birth and death, adopted by [START_REF] Hubbell | The Unified Neutral Theory of Biodiversity and Biogeography (MPB-32)[END_REF] and widely studied in UNTB-derived models (see Rosindell et al., 2011 and references therein), constitutes a model of overlapping generations akin to the discrete-time Moran process from population genetics theory [START_REF] Moran | Random processes in genetics[END_REF].

We consider three distinct community assembly models which represent different processes relevant for ecological modeling: neutrality, competitive exclusion, and environmental filtering.

The differences between these three models are realized by differentially weighting individual death probabilities per species with respect to their trait values. In the simplest case (the 'neutral' model), all species are ecologically equivalent, so that individual probability of death is 1/J for all individuals, regardless of species identity or trait value. The other two assembly models capture the importance of biotic or abiotic interaction (which we collectively term 'ecological interactions'), reflecting models of competitive exclusion [START_REF] Macarthur | The Limiting Similarity, Convergence, and Divergence of Coexisting Species[END_REF] and ecological filtering [START_REF] Bazzaz | Habitat Selection in Plants[END_REF], respectively. Death probabilities are all the more affected by individuals' traits that the strength of ecological interactions (es; or 'ecological strength') is high. Together we refer to these as 'non-neutral' models because within each, the neutral assumption of ecological equivalence is violated. Both non-neutral models are adapted from coevolutionary models relating probability of persistence with trait mediated interactions among lineages [START_REF] Nuismer | Coevolution and the architecture of mutualistic networks[END_REF][START_REF] Nuismer | Predicting rates of interspecific interaction from phylogenetic trees[END_REF], which have been previously adapted to study ecological interactions at macroevolutionary [START_REF] Aristide | Understanding the effect of competition during evolutionary radiations: an integrated model of phenotypic and species diversification[END_REF][START_REF] Drury | Estimating the Effect of Competition on Trait Evolution Using Maximum Likelihood Inference[END_REF][START_REF] Drury | An Assessment of Phylogenetic Tools for Analyzing the Interplay Between Interspecific Interactions and Phenotypic Evolution[END_REF] and community ecological scales (Overcast et al., 2020;[START_REF] Ruffley | Identifying models of trait-mediated community assembly using random forests and approximate Bayesian computation[END_REF]. In the model of environmental filtering (hereafter 'filtering'), we express the probability of death of an individual from species i as:

𝜹 ! = 1 -𝑒𝑥𝑝[-𝑒 " (𝑧 ! -𝑧 # ) 2 ]
where zi is the trait value of species i, zE is the optimal trait value within the environment, and es is the strength of ecological interactions. This model favors individuals with trait values closer to the environmental optimum, and so should reduce the variance in local trait values, and increase both age and abundance of increasingly 'fit' species. In the model of competitive exclusion (hereafter 'competition'), the probability of death of an individual from species i is:

𝜹 ! = 𝑒𝑥𝑝[-𝑒 " (𝑧 ! -𝑧) 2 ]
with similar terms as above, excepting 𝑧 which indicates the local community trait mean. This model favors individuals with trait values that are farther from the local trait mean, in essence, individuals that dwell in unoccupied niche space, and should increase the variance in local trait values (see for example [START_REF] Govaert | Measuring the contribution of evolution to community trait structure in freshwater zooplankton[END_REF]. The functional form of both equations are structured such that the effect of different es values are comparable between the two non-neutral models.

Increasing values of es correspond to increasing departures from neutrality, with trait differences among species having a strong impact on death probability. In practice, as for most of the other model parameters, eS will rarely be known a priori, but rather will be estimated from the data (see [START_REF] Drury | An Assessment of Phylogenetic Tools for Analyzing the Interplay Between Interspecific Interactions and Phenotypic Evolution[END_REF][START_REF] Ruffley | Identifying models of trait-mediated community assembly using random forests and approximate Bayesian computation[END_REF]. At each time step, an individual is sampled for death following a multinomial distribution parameterized by the death probabilities characterizing each species, weighted by their local abundance.

New species may arise within the local community by speciation, which we model in two ways, either through divergence in isolation or through in situ cladogenesis. In the first case, we consider a variant of the protracted speciation model, whereby migrant lineages diverge from their metacommunity sister lineages given sufficient passage of time (following Rosindell & Phillimore, 2011 and the notation therein). Upon migration of a new species unique within the local community, we initially consider this lineage to be a "population" of the metacommunity species, as insufficient time has elapsed for it to form a proper endemic species. If the local population of this lineage persists for a given amount of time (τ generations; i.e. the speciation duration), then we recognize it as a new species, and reassign all individuals of this lineage to a new species identity. Recognizing that repeated migration from the metacommunity may dampen the speciation process by preserving genetic connectivity, for each post-colonization migration event the τ value of the focal species is incremented by a small fixed amount (G; the gene flow effect).

Therefore, given a sufficiently high rate of post-colonization migration, gene flow will effectively arrest the speciation process. Upon speciation completion, a new branch is grafted onto the branch of the ancestral species in the metacommunity phylogeny, with branch length equal to the realized speciation duration (including time accumulated through the gene flow effect), and speciation rate determined by the rate of the immediately ancestral branch (λi), as modulated by the ClaDS model parameters. The presence of this new branch in the metacommunity phylogeny incorporates it within the global diversification process, and it may now speciate in situ by random fission, with individuals randomly assigned to each of the two new species.

The process of individual death and replacement by birth or migration, and speciation continues in this manner until one of two stopping conditions are satisfied (tn in Fig. 1), either 1) a fixed amount of time has elapsed (as measured in generations); or 2) a threshold equilibrium state has been exceeded. As stated above, in the local community process, generations are measured in terms of J timesteps, effectively controlling the timescale for variation in local community size among simulations. In other words, communities i and j of very different sizes (Ji << Jj) will experience a similar degree of change in local community composition over the course of the same number of generations. The second stopping criterion threshold is the "equilibrium state" of the community (Λ; [START_REF] Overcast | An integrated model of population genetics and community ecology[END_REF], following Rosindell & Harmon 2013), which we define as the complement of the proportion of the signal of initial community state retained in the present state (see [START_REF] Overcast | An integrated model of population genetics and community ecology[END_REF] section 2.5 for further treatment). Equilibrium state (Λ) can range from 0 to 1 (inclusive), with small values characterizing young communities, and values approaching 1 characterizing older communities. In practice, because community assembly is a stochastic process, measuring community age as equilibrium state, rather than as time in generations, tends to better reflect the true underlying community history [START_REF] Overcast | A unified model of species abundance, genetic diversity, and functional diversity reveals the mechanisms structuring ecological communities[END_REF].

Upon completion of the forward-time community assembly process we undertake a backwardtime coalescent simulation per species in the local community to generate community-scale patterns of genetic variation [START_REF] Kelleher | Efficient Coalescent Simulation and Genealogical Analysis for Large Sample Sizes[END_REF]. Abundance histories and colonization times per species are tracked during the forward-time process, then rescaled to effective population size (following the α parameter of [START_REF] Overcast | A unified model of species abundance, genetic diversity, and functional diversity reveals the mechanisms structuring ecological communities[END_REF] and divergence time to parameterize each coalescent simulation. To reduce computational complexity, we calculate the harmonic mean of effective population size trajectories per species, a transformation for populations with fluctuating sizes which is justified if reasonable conditions are met [START_REF] Karlin | Rates of Approach to Homozygosity for Finite Stochastic Models with Variable Population Size[END_REF]). Additionally, forward-time post-colonization migration events are recorded and the frequency of these is used to calculate a continuous migration probability, to inform the model of divergence with a rate of ongoing gene flow. Finally, given these parameters, one coalescent tree is simulated per species and mutations are applied to it following an infinite-sites model [START_REF] Kelleher | Efficient Coalescent Simulation and Genealogical Analysis for Large Sample Sizes[END_REF]. For simulation of the genetic data, the mutation rate and simulated sequence length are global parameters, and assumed to be fixed and homogenous for all species. They will typically be fixed to match characteristics of the metabarcoding data, with an underlying assumption that the simulated gene(s) is/are neutral.

Statistical inference, raw data, and summary statistics

We use a machine learning approach to perform classification (i.e. model selection) and parameter estimation (scikit-learn; [START_REF] Pedregosa | Scikit-learn: Machine Learning in Python[END_REF]. Of the available models within scikitlearn, we tested both random forest (RF; [START_REF] Breiman | Random Forests[END_REF] and gradient boosting (GB; [START_REF] Friedman | Greedy function approximation: a gradient boosting machine[END_REF] methods, and used a random search to identify the best machine learning algorithm hyperparameters [START_REF] Bergstra | Random search for hyper-parameter optimization[END_REF]. We generated 10,000 simulations under each of the neutral, environmental filtering, and competition models, and we tested the machine-learning approach on both raw data output and on summary statistics, focusing on the local community.

The community-scale raw data output of any given simulation includes the abundance, trait value, and nucleotide diversity (π; [START_REF] Nei | Mathematical model for studying genetic variation in terms of restriction endonucleases[END_REF] for each species within the local community, as well as a phylogeny pruned to contain only local species. The local phylogeny is encoded using a novel tree representation format, Compact Bijective Ladderized Vectorization (CBLV; [START_REF] Voznica | Deep learning from phylogenies to uncover the transmission dynamics of epidemics[END_REF], which has been shown to have useful statistical properties in machine learning applications. In order to design summary statistics, we used the framework of Hill numbers on the abundances, trait values, and genetic diversities; Hill numbers provide a unified way to quantify the evenness and magnitude of different axes of biodiversity data [START_REF] Chao | Unifying Species Diversity, Phylogenetic Diversity, Functional Diversity, and Related Similarity and Differentiation Measures Through Hill Numbers[END_REF]). Here we focus on genetic diversity Hill numbers of order one ( 1 GD), two ( 2 GD), and three ( 3 GD), whereby Hill numbers of increasing order place increasing weight on diversity values of higher magnitude.

By design, the simulated data can be partitioned in numerous ways, to allow for direct comparability with the particular biodiversity data available within the empirical system under investigation. For example, if only genetic data and trait values are available, then the local phylogeny and abundances of the simulated data can be masked out, to allow performing inference using the exact data in hand. In many cases metabarcoding datasets may be limited to only having community-scale genetic diversities, and our model and inferential framework can be used in such cases, as well.

Model validation simulation experiments

We evaluate the precision and accuracy of classification and parameter estimation using simulation experiments. We specified the form of the simulated data (i.e. number of OTUs, number of individuals per OTU, length and mutation rate of barcode sequence, etc.) to match that of our empirical metabarcoding example, to better understand uncertainty in inference given the constraints of the real data. Parameters of the simulation experiments were chosen based on previous estimates for these from literature and/or to adhere as closely as possible to those most relevant for application to the empirical metabarcoding data. We set loguniform priors on ClaDS σ (0.1-0.4), ClaDS α (0.6-0.9), the size of the local community (40,000-150,000), colonization rate (0.0001-0.005), ecological strength (0.1-10), and equilibrium state stopping time for the local community assembly process (0.5-1.0), fixing all other parameters to biologically reasonable values (see Table S1 for the full list of parameter values used). Additionally, we applied a transformation to the simulated sequence data in order to more accurately represent the empirical metabarcoding data that we use as a test case in this manuscript (see the following section). The empirical test data retained only amplicon sequencing variants (ASVs; i.e. unique haplotypes; see [START_REF] Callahan | Exact sequence variants should replace operational taxonomic units in marker-gene data analysis[END_REF] clustered into operational taxonomic units (OTUs; putative species). To better represent this, we deduplicated simulated sequences (i.e. we removed all but one copy of identical haplotypes) within each species prior to calculation of π for all species in the local community. This sequence deduplication step is optional and can be omitted when generating simulations for empirical datasets for which reliable haplotype abundances are available, potentially improving performance.

In our experimental approach, we evaluated the precision and accuracy of our machine learningbased inference procedure under several possible combinations of data availability, to better understand how each biodiversity data axis contributes to both assembly model classification and parameter estimation. In one set of experiments we used abundances, trait values, and genetic diversities per species in the local community, as well as the CBLV representation of the local community phylogeny to investigate machine learning inference under the best-case data availability scenario. Assuming that each data axis retains some unique signal of community assembly history, this represents an upper bound on power and accuracy. In a second set of experiments we used only genetic data from the local community, to better understand machine learning performance given the data available for our empirical example. In each of these sets of experiments we also evaluated machine learning performance using community-scale genetic diversity either as raw π values or as summarized by the first three genetic diversity Hill numbers (as described above). When utilizing raw π values as the data (for both cross-validation experiments and empirical analysis), we first sorted the vector of π values within each simulation, in order of decreasing magnitude. This ordering scheme is a common method for reducing combinatorial sample space and improving computational efficiency during simulation based inference where the distribution of π values is of primary interest, and not the specific π value of any given OTU (e.g. [START_REF] Overcast | Strategies for improving approximate Bayesian computation tests for synchronous diversification[END_REF].

Model classification performance was visually inspected using confusion matrices [START_REF] Susmaga | Confusion Matrix Visualization[END_REF]) and was quantified using precision, recall, and F1 measures [START_REF] Lipton | Optimal Thresholding of Classifiers to Maximize F1 Measure[END_REF]. Precision (complement of the false positive rate), recall (complement of the false negative rate), and F1 (the harmonic mean of precision and recall) all take values between 0 and 1, with values approaching 1 indicating increasingly improved model performance. Parameter estimation performance was evaluated independently for simulations from each of the assembly models, given that assembly models may generate drastically different distributions of data given identical input parameters.

Parameter estimation accuracy was visualized by plotting true values against predicted values, and was quantified by computing the coefficient of determination using a cross-validation approach. Classification and regression performance were both evaluated by averaging scores using 5-fold cross-validation, and feature importance was evaluated using a permutation approach, to better control for bias in feature relevance [START_REF] Altmann | Permutation importance: a corrected feature importance measure[END_REF].

Application to Cyprus soil microarthropod metabarcode data

As an empirical system, we focus on the soil microarthropod assemblages of the montane forests of Troodos within the oceanic island of Cyprus. These assemblages form a highly isolated montane metacommunity distributed across a mosaic of different forest types over a topographically complex landscape (Fig. 4a). We used the community-scale metabarcoding data of [START_REF] Noguerales | Community metabarcoding reveals the relative role of environmental filtering and spatial processes in metacommunity dynamics of soil microarthropods across a mosaic of montane forests[END_REF] who extracted microarthropods sampled from soil at 44 sites and targeted a 418 bp fragment of the cytochrome c oxidase subunit I (COI) barcode region. Raw reads were subject to a stringent filtering strategy, resulting in 907 ASVs, corresponding to 386 OTUs after 3% sequence similarity clustering. The average nucleotide diversity within OTUs was 0.003 (std = 0.006). Correlation analysis with spatial and environmental variables found a joint effect of neutral processes (dispersal) and environmental filtering in driving metacommunity structure, which together with the extensive dataset, constitute an interesting system to apply our model. Specifically, the analyses of [START_REF] Noguerales | Community metabarcoding reveals the relative role of environmental filtering and spatial processes in metacommunity dynamics of soil microarthropods across a mosaic of montane forests[END_REF] showed a primary role of forest type in structuring regional soil microarthropod metacommunities, with largely distinct community composition among the following four forest types: (i) the Calabrian pine forest-Pinus brutia (Pb) the dominant forest type, forming continuous and extensive forests across Troodos; (ii) the endemic golden oak-Quercus alnifolia (Qa), with a broad but highly fragmented distribution across Troodos; (iii) the narrow endemic Cyprus cedar-Cedrus brevifolia (Cb), with a highly restricted distribution in Western Troodos; and finally (iv) the black pine-Pinus nigra pallasiana (Pn) forest merged together with the stinking juniper-Juniperus foetidissima forest, as they are both narrowly distributed at the top of the highest peak in Central Troodos and have similar community composition [START_REF] Noguerales | Community metabarcoding reveals the relative role of environmental filtering and spatial processes in metacommunity dynamics of soil microarthropods across a mosaic of montane forests[END_REF]. Because these communities were found to structure more by forest type than by geographic space, we analyze each of the four forest types separately, leveraging this dataset as an illustrative test case to investigate how elevation, habitat distribution, and connectivity impact community composition. For the empirical analysis we generated 20,000 simulations under each of the three community assembly models, sampling parameters uniformly from prior ranges (Table S1). We used principal component analysis to perform a prior predictive check of the appropriateness of the sampled parameters by randomly selecting 1000 simulations and projecting both simulated and empirical nucleotide diversity values into PC space. Visual overlap in PC space indicates that simulations can accurately capture important properties of the empirical data. We then applied our machine learning classification procedure to identify the most probable community assembly model per forest type. Following this, we retained simulations for the most probable model and performed parameter estimation using gradient boosting and quantile regression [START_REF] Waldmann | Quantile regression: A short story on how and why[END_REF], to characterize prediction intervals (PIs).

Results

Model validation simulation experiments

Assembly model classification showed higher accuracy given the full complement of data, though classification remained reasonably accurate when constrained to using genetic data alone (Table 1). In all cases, gradient boosting methods outperformed random forest, so only results using this method are reported. Classification based on raw data values was slightly improved, with respect to that based on summary statistics (Table 1; Fig. 2), indicating that the machine learning method was better able to extract information from unprocessed data [START_REF] Schrider | Supervised Machine Learning for Population Genetics: A New Paradigm[END_REF]. Within classification experiments there was variation in precision, recall, and F1 among community assembly models, with competition more accurately predicted, neutrality less so, and filtering obtaining intermediate accuracy. Reduced accuracy of neutral models can be attributed to the chosen range of values adopted for the ecological strength parameter in simulations, as nonneutral assembly models with very low ecological strength may be indistinguishable from purely neutral models. Evaluation of feature importance for classification based on summary statistics showed decreasing importance with increasing order of genetic diversity Hill numbers (e.g. 1 GD was most important and 3 GD least; Fig. S1a). Recalling that when using raw π values for machine learning the 'features' correspond to the vector of π values per OTU ordered from greatest to least magnitude, evaluation of feature importance identified primarily OTUs with intermediate values of π as the most important for distinguishing community assembly models (Fig. S1b).

Parameter estimation accuracy varied substantially across tested parameters, assembly models, and data combinations (Table 1; Figs. 3,S2). Again, random forest underperformed with respect to gradient boosting methods in all cases (e.g. Fig. S3), so only the latter results are reported.

Colonization rate was the best estimated parameter across all treatments (global mean R 2 =0.7; std=0.13). Equilibrium state was generally better estimated under competition models (mean R 2 =0.42, max. R 2 =0.58), and ecological strength better estimated under filtering models (mean R 2 =0.15; max. R 2 =0.19). The most accurate estimates of equilibrium state were obtained under neutral models when using summary statistics on the full complement of data axes (R 2 =0.87).

Metacommunity phylogeny parameters could not be estimated with confidence given only genetic data from the local community (neutral models R 2 ~0), though non-neutral models showed small positive R 2 values for ClaDS α (trend of shift in speciation rate) and ClaDS m (relative change in speciation rate) indicating promising, but weak signal in this data (R 2 ~0.08). Utilizing all available local data dramatically improved estimation of ClaDS α (R 2 neutral=0.78; competition=0.67; filtering=0.74) and ClaDS m (R 2 neutral=0.8; competition=0.69; filtering=0.75). On the other hand, ClaDS σ (stochasticity of rate shift) was extremely difficult to estimate from local community data, even when using all available data axes, including the local phylogeny (global mean R 2 =0; std=0.02; max=0.04). In most cases parameter estimation cross-validation procedures using raw data values outperformed those using Hill numbers as summary statistics, whether considering the full complement of data or when limited to local genetic diversity data alone (Table 1). Feature importance for parameter estimation regression followed broadly similar patterns as with classification (see above).

Application to Cyprus soil macroinvertebrate metabarcode data

Machine learning classification of the Cyprus soil macroinvertebrate communities identified neutral processes as driving communities from broadly distributed, low-mid elevation forests (Pb, Pinus brutia and Qa, Quercus alnifolia) and environmental filtering as driving communities from high elevation, geographically restricted forests (Cb, Cedrus brevifolia and Pn, Pinus nigra). Pinus brutia had higher support for ecological neutrality (0.84 prediction probability) with respect to Quercus alnifolia (0.75 prediction probability), while the remaining fraction of non-neutral prediction for both forests was roughly split between competition and filtering models (Pb: 0.09/0.07; Qa: 0.11/0.15 for competition and filtering, respectively). Cedrus brevifolia and Pinus nigra both received greater than 0.95 prediction probability supporting ecological filtering as the generating model (Cb: 0.98; Pn: 0.97).

Parameter estimation was performed independently for each forest type, restricting the simulations used for inference to only those which corresponded to the most probable assembly model per forest (e.g. Pb and Qa used neutral simulations). Communities classified as neutral were predicted to have been sampled from larger local communities (Table S2; J: Pb=117561; 110492-139455 50% PI) with higher colonization rates (Pb=0.0036; 0.0026-0.0045 50% PI) than communities structured by environmental filtering (J: Cb=81995; 70046-109969 50% PI; colonization rate: Cb=0.0019; 0.0012-0.0027 50% PI). All communities showed broad overlap in predicted equilibrium state and ecological strength, as expected from the limited ability to infer these parameters, though neutral community PIs (e.g. Pb=0.69-0.95 50% PI) were slightly higher than the non-neutral communities (e.g. Cb=0.61-0.93 50% PI), and the narrow endemic Cb forest generated slightly higher predicted ecological strength (2.41; 0.54-5.49 50% PI) than the high elevation Pn forest (1.01; 0.56-2.95 50% PI). Prediction intervals of ClaDS α and ClaDS m were also broadly overlapping for all communities (Table S2), as expected given that there is little information about these parameters in genetic data alone (Table 1; Fig. 4).

Discussion

In this manuscript we developed a new eco-evolutionary model that can be used in combination with empirical data to test alternative community assembly hypotheses and estimate key biological parameters and which makes several advances over current state-of-the-art (see [START_REF] Overcast | Towards a genetic theory of island biogeography: Inferring processes from multidimensional community-scale data[END_REF]. We make the model and associated inference tool available as a userfriendly package available in both R and python. We have shown that neutral and non-neutral community assembly processes can be distinguished using sequence data sampled from the local community, with inference accuracy increasing with increasing amounts of available community-scale data (e.g. abundances, trait values, phylogenies). We have also shown that key ecological and evolutionary parameters may be inferred with such data, that there is variability in uncertainty of inference given the parameters of interest, and again that the precision and accuracy of inference are a function of the available empirical data. We note explicitly that the iBioGen package is targeted toward investigation of metabarcoding data (see [START_REF] Creedy | Coming of age for COI metabarcoding of whole organism community DNA: Towards bioinformatic harmonisation[END_REF][START_REF] Arribas | Toward global integration of biodiversity big data: a harmonized metabarcode data generation module for terrestrial arthropods[END_REF], and we have shown that this individual axis of data contains valuable information on the processes shaping empirical communities.

As ecological and evolutionary processes are necessarily intertwined, community ecologists and biodiversity modelers have increasingly turned to incorporating alternative axes of data, such as phylogenies and genetic variation, into investigations of community assembly dynamics.

Community phylogenetic methods are prominent among these [START_REF] Webb | Phylogenies and Community Ecology[END_REF], having demonstrated power to distinguish niche-based from neutral processes [START_REF] Kraft | Functional traits and niche-based tree community assembly in an Amazonian forest[END_REF], though these tend to take the observed phylogeny as data and to adopt correlational, rather than predictive, model-based approaches (but see, for example [START_REF] Ruffley | Identifying models of trait-mediated community assembly using random forests and approximate Bayesian computation[END_REF]. Considering speciation processes within models derived from UNTB can make predictions of phylogenies that closely resemble the structure of empirical phylogenies (e.g. [START_REF] Jabot | Inferring the parameters of the neutral theory of biodiversity using phylogenetic information and implications for tropical forests[END_REF], improving the estimation of key biodiversity parameters. However, these models are not appropriate for capturing diversification rate variation among lineages, or the temporal dynamics of speciation [START_REF] Davies | Neutral biodiversity theory can explain the imbalance of phylogenetic trees but not the tempo of their diversification[END_REF]. Our linked model of metacommunity diversification and local community assembly predicts local-scale biodiversity patterns as a function of variable metacommunity diversification dynamics, and reciprocally may be used to estimate parameters of the macroevolutionary process from local-scale data. The speciation process, and thus phylogenies, may also be modeled in a bottom-up fashion as a function of both demographic (population size changes and migration) and mutational processes, considering the accumulation of genetic differentiation among lineages using an individual-based [START_REF] Manceau | Phylogenies support out-of-equilibrium models of biodiversity[END_REF] or populationbased [START_REF] Marin | The genomic view of diversification[END_REF] approach. Linked patterns of community-scale genetic variation and phylogenetic relationships may be predicted by these models, but have yet to be fully explored.

Here we explicitly undertake this investigation, while also fitting our model to empirical metabarcoding data, using community-scale genetic variation to infer ecological and evolutionary processes.

Our model broadly follows that of [START_REF] Overcast | A unified model of species abundance, genetic diversity, and functional diversity reveals the mechanisms structuring ecological communities[END_REF], though with a few important differences.

First, in our metacommunity, abundances are phylogenetically informed as they are an outcome of the macroevolutionary process, whereas in [START_REF] Overcast | A unified model of species abundance, genetic diversity, and functional diversity reveals the mechanisms structuring ecological communities[END_REF] the metacommunity phylogeny is generated with a birth/death process and abundances are randomly assigned at the tips following a log-series distribution [START_REF] Fisher | The relation between the number of species and the number of individuals in a random sample of an animal population[END_REF]. Second, our local community model can optionally follow a Wright-Fisher (WF;[START_REF] Fisher | XXI.-On the Dominance Ratio[END_REF][START_REF] Wright | Evolution in Mendelian Populations[END_REF] process with nonoverlapping generations, whereas the local community in [START_REF] Overcast | A unified model of species abundance, genetic diversity, and functional diversity reveals the mechanisms structuring ecological communities[END_REF] (and studied in our model here) utilizes a Moran model with overlapping generations [START_REF] Moran | Random processes in genetics[END_REF]. In the limit of infinite population size, these two processes are identical given an appropriate rescaling of time, but the WF process our model may adopt greatly improves computational efficiency, allowing simulation of substantially larger local communities in a reasonable amount of time (results not shown). However, characterizing the behavior of the WF implementation is an avenue for future research. Finally, in our model we allow for genetic variation to be calculated based only on unique haplotypes within a given species, facilitating its adoption by researchers implementing community metabarcoding protocols which are able to resolve individual haplotypes, but which typically cannot (at present) accurately estimate haplotype abundances.

Estimating metacommunity parameters from local community data

Given community-scale data, including the CBLV representation of the pruned local phylogeny, we have shown that sufficient information remains within it to adequately infer parameters of the macroevolutionary process (Fig. S2). Here, the question of estimating macroevolutionary parameters from local-scale data resembles the question of estimating diversification rates in the presence of incompletely sampled phylogenies, a problem which is well studied [START_REF] Morlon | Reconciling molecular phylogenies with the fossil record[END_REF][START_REF] Stadler | How can we improve accuracy of macroevolutionary rate estimates?[END_REF]. Indeed, [START_REF] Voznica | Deep learning from phylogenies to uncover the transmission dynamics of epidemics[END_REF] showed that parameters of the macroevolutionary process were well estimated by a neural network approach using their CBLV tree representation, while informing the model with the empirical sampling fraction. In our case, the sampling fraction is implicit within the simulations used as training data, given that simulated metacommunity trees were generated with identical numbers of tips, and the sampling fraction amounts to local species richness. It is interesting to note that the estimates of macroevolutionary parameters (ClaDS α and ClaDS m) are more accurate under a neutral community assembly model, with respect to the two non-neutral models. This may be explained by the fact that the sample of lineages in the local community is unbiased under a neutral model, whereas under nonneutral models there is a bias in this sampling induced by the phylogenetically informed trait values and the local community non-neutral dynamics. A question of perhaps greater interest is the nature of the signal of the macroevolutionary process (albeit faint) in the local-scale genetic data, specifically for non-neutral assembly models (Table 1; Fig. 2). This may be explained by a knock-on effect of diversification dynamics impacting the phylogenetic distribution of trait values, which promote or impede population persistence and abundance in the local community, and ultimately are reflected in patterns of genetic variation. This explanation is evidenced by the absence of information about metacommunity diversification parameters in the local-scale genetic diversity of neutral assembly models.

Processes structuring Cyprus soil macroinvertebrate communities

Machine learning analysis based on our model simulations revealed how habitat-specific communities can be governed by contrasting assembly processes, largely in concordance with previous studies based on correlation analyses [START_REF] Noguerales | Community metabarcoding reveals the relative role of environmental filtering and spatial processes in metacommunity dynamics of soil microarthropods across a mosaic of montane forests[END_REF]. While communities from broadly distributed, low-mid elevation habitats were primarily determined by neutral processes (Pb and Qa), our model estimated a major effect of environmental filtering in high elevation and geographically restricted habitats (Pn and Cb, respectively). A similar pattern has been observed in bird communities distributed along a habitat and elevational gradient [START_REF] García-Navas | Spatial heterogeneity in temporal dynamics of Alpine bird communities along an elevational gradient[END_REF], where stochastic factors were found to be more influential at low-elevation communities.

In our case, this difference could be explained by the widespread distribution and mild climate conditions exhibited by Quercus alnifolia and Pinus brutia forests, both distributed at low-mid elevations and exhibiting higher environmental variation among plots. Communities within these habitats may be composed of more generalist/abundant species, likely with broader altitudinal ranges and wider environmental tolerances, and thus larger geographic distributions determined primarily by stochastic processes at the regional scale, resulting in a more uniform distribution of π values within the community (Fig. 4b). The marked habitat fragmentation of the Quercus alnifolia forests could also contribute to generating stochastic patterns [START_REF] Noguerales | Community metabarcoding reveals the relative role of environmental filtering and spatial processes in metacommunity dynamics of soil microarthropods across a mosaic of montane forests[END_REF], as documented also in herbaceous metacommunities distributed across fragmented forests [START_REF] Jamoneau | Fragmentation alters beta-diversity patterns of habitat specialists within forest metacommunities[END_REF]. Conversely, the narrow environmental conditions of high-elevation or geographically restricted habitats (Pn and Cb), which are characterized by reduced environmental variation among plots (see Fig. S1 in [START_REF] Noguerales | Community metabarcoding reveals the relative role of environmental filtering and spatial processes in metacommunity dynamics of soil microarthropods across a mosaic of montane forests[END_REF]) may impose a strong abiotic filter on species, resulting in a non-neutral signature in community composition [START_REF] Andrew | Beta-diversity gradients of butterflies along productivity axes[END_REF].

Those communities mostly consist of specialist species with small geographic ranges, whose population dynamics may be constrained by climatic and historical processes (Fig. 4; [START_REF] Tonzo | Glacial connectivity and current population fragmentation in sky islands explain the contemporary distribution of genomic variation in two narrow-endemic montane grasshoppers from a biodiversity hotspot[END_REF][START_REF] Polato | Genetic diversity and gene flow decline with elevation in montane mayflies[END_REF]. Indeed, signatures of community genetic variation within highaltitude communities, as exhibited by a skewed distribution of π values (few species with high diversity and most species with little), are consistent with the restricted immigration and ecological selection which can be expected under a model of environmental filtering. Evidence of environmental filtering in high-altitude environments has been found in a wide array of taxa including freshwater fish communities [START_REF] Kirk | Environmental filters of freshwater fish community assembly along elevation and latitudinal gradients[END_REF], bees [START_REF] Hoiss | Altitude acts as an environmental filter on phylogenetic composition, traits and diversity in bee communities[END_REF], and even gut microbiota [START_REF] Li | Environmental filtering increases with elevation for the assembly of gut microbiota in wild pikas[END_REF], and our results contribute to this growing body of literature.

These results emphasize the idiosyncratic dynamics of community assembly in environmentally complex contexts (e.g. [START_REF] Gomez | The biotic interactions hypothesis partially explains bird species turnover along a lowland Neotropical precipitation gradient[END_REF][START_REF] Gweon | Contrasting community assembly processes structure lotic bacteria metacommunities along the river continuum[END_REF], while simultaneously exemplifying the inferential capability of our model to discern between community assembly forces. The aforementioned inferences are in broad concordance with previous findings based on correlation analyses for most of the forest types, except for the Pinus brutia forest in which previous analysis revealed a primary effect of deterministic processes [START_REF] Noguerales | Community metabarcoding reveals the relative role of environmental filtering and spatial processes in metacommunity dynamics of soil microarthropods across a mosaic of montane forests[END_REF] This discrepancy between previous correlational analysis and our current model based inference could potentially be a result of the difference in spatial scale of the approaches. While the present study predicts the processes structuring community composition within habitat types at the regional scale, [START_REF] Noguerales | Community metabarcoding reveals the relative role of environmental filtering and spatial processes in metacommunity dynamics of soil microarthropods across a mosaic of montane forests[END_REF] studied the factors determining beta diversity among sites within the region. Future work examining processes within and among habitats at the regional scale, and/or incorporating additional axes of biodiversity data (e.g. phylogenies), may provide further insight to unify these findings [START_REF] Liu | Mountain metacommunities: climate and spatial connectivity shape ant diversity in a complex landscape[END_REF].

Future directions

We can envision many different ways to develop the set of models included in the iBioGen package. At the scale of the metacommunity, numerous extant modeling approaches can be envisioned to increase the biological realism of the simulations, including constraining covariation among model relevant traits (i.e. body size and population growth rate; [START_REF] Clavel | Mv morph : An r package for fitting multivariate evolutionary models to morphometric data[END_REF]. An allometric theory may be introduced by conditioning colonization rate and probability of local extinction on trait values (again imagining the focal trait as body size; [START_REF] Jacquet | Extensions of Island Biogeography Theory predict the scaling of functional trait composition with habitat area and isolation[END_REF]. Future work may also incorporate more biologically relevant non-neutral models which represent competition as pairwise among lineages [START_REF] Luiselli | Detecting the ecological footprint of selection[END_REF]. From a strictly operational perspective, both the Hill number summarization and the computationally efficient sorting procedure result in breaking the connection between the phylogeny and the other axes of local data (abundances, traits, and genetic diversities). Carefully reconsidering the composition of the simulations to retain the phylogenetic structure of the local data may result in a gain in precision and accuracy for machine learning inference. Additionally, while our simulation experiments fixed some parameter values to benefit computational tractability, fixing parameters to reasonable (though arbitrary) values may bias or artificially improve the inference, in which case further computational effort should be devoted to testing model performance using vastly more simulations and allowing flexibility in all parameters to better assess the accuracy of our approach.

Extensions of the iBioGen model may also be suitable for investigating past population dynamics of whole communities or assemblages, such as community-wide (or some fraction thereof) simultaneous bottlenecks or expansions, as might be expected under a model of historical climate change and expansion/contraction of suitable habitat (following e.g. [START_REF] Hickerson | msBayes: pipeline for testing comparative phylogeographic histories using hierarchical approximate Bayesian computation[END_REF][START_REF] Overcast | Strategies for improving approximate Bayesian computation tests for synchronous diversification[END_REF]. Such tests of shared (synchronous) demographic responses, indeed all inference of community-scale processes, will only gain in power as metabarcoding datasets expand beyond contemporary single locus experiments (e.g. [START_REF] Noguerales | Community metabarcoding reveals the relative role of environmental filtering and spatial processes in metacommunity dynamics of soil microarthropods across a mosaic of montane forests[END_REF], to multilocus DNA metabarcoding [START_REF] Sard | Comparison of fish detections, community diversity, and relative abundance using environmental DNA metabarcoding and traditional gears[END_REF], to whole-genome metagenomics [START_REF] Pérez-Cobas | Metagenomic approaches in microbial ecology: an update on whole-genome and marker gene sequencing analyses[END_REF], and preparing to leverage this increasing depth of sequencing within the model will be a critical avenue for future research.

Conclusions

Here we have presented the iBioGen package, a toolkit for studying community-scale genetic data (such as might be obtained from a typical metabarcoding experiment) to understand the ecological and evolutionary processes structuring biological communities. The iBioGen package can accommodate one, few or all axes of biodiversity data including abundances, trait values, genetic diversities, and phylogenetic relationships, allowing researchers great flexibility to study their systems given their (potentially limited) data in hand. We have demonstrated that, while inference is most powerful given the full complement of input data, reasonable inference can still be obtained with limited data. This is particularly true if the focus of study is inferring neutrality or non-neutrality of the community assembly process, as the single axis of genetic data can distinguish these with sufficiently high accuracy. Thus, the iBioGen package provides a powerful, model-based approach for inferring processes from community metabarcoding patterns, which may complement and greatly extend insights obtained from contemporary correlative approaches.

Figure 1 -Conceptual diagram of iBioGen community assembly simulation process

A diagram representing the dynamics of the iBioGen community assembly model including the joint metacommunity (sensu Hubbell, 2001, a pool of evolving species evolving through speciation and extinction events) and local community processes. The simulation model proceeds in two stages. In the first stage, the metacommunity phylogeny is evolved following a lineage-based speciation/extinction process allowing for shifts in speciation rate at branching events. Abundances additionally evolve along the branches as a function of species specific population growth rates. The abundance of each lineage is indicated by the colored circles at each node, with abundant species in red and rare species in blue. Negative growth rates can result in abundance falling below a viable population threshold leading to extinction (not depicted). Trait values (e.g. body size) also evolve along the branches following Brownian motion, and are indicated in the colored squares adjacent to each tip. The metacommunity phylogeny is evolved until a fixed time has elapsed or a fixed number of lineages are extant (here indicated as T0). In the second stage, the local community is initialized and the metacommunity and local community evolve in lockstep. The local community is of constant size (J individuals) and is spatially implicit (the grid depiction in the figure is for clarity in visual illustration). Between speciation events (branching events in the metacommunity phylogeny), or "epochs', the local community assembles following the birth and death of individuals and colonization from the metacommunity. Speciation can occur in the metacommunity, as well as in the local community by divergence in isolation or in situ speciation (not represented in the figure). The simulation proceeds for a fixed amount of time (n generations), and terminates at time Tn. Three different, independent local assembly models are considered: neutral, environmental filtering, and competition. In the neutral case (top), all individuals of all species are ecologically equivalent and trait values and the environmental context do not contribute to probability of death. In a model of environmental filtering (middle), the environmental context (indicated by the blue background) imposes a constraint on individual death probabilities which are a function of the individual's trait, with trait values closer to the environmental optimum conferring fitness benefits. In the competition model, the environmental context is inconsequential and instead individual death probabilities are a function of the distance of their trait value from all other trait values in the local community. Individuals that are more different have more ecological space and therefore have an increased probability of surviving. . This cross-validation procedure was performed using only raw nucleotide diversity per species in the local community as the test data. Please note that R 2 values in the figure were calculated for the exact cross-validation fold depicted, so these will differ slightly from the average R 2 over all cross-validation folds as shown in Table 1.

Figure 2 -

 2 Figure 2 -Machine learning cross-validation confusion matrix for model classification Confusion matrices indicating the accuracy of our machine learning-based simulation approach. Numerical values within each bin indicate the probability of class label predictions for each true class label, and squares are filled to indicate increasing probability of correct prediction, from low (white) to high (black). Correct model classification falls along the diagonal, with off-diagonal labels indicating incorrect prediction. In a) the cross-validation classification procedure was performed using summary statistics of simulated community genetic data (Hill numbers), whereas in b) the raw nucleotide diversity values per species were used.

Figure 3 -

 3 Figure 3 -Machine learning cross-validation parameter estimation for simulations under a model of environmental filtering Machine learning cross-validation simulations demonstrate the bias and accuracy in estimating model parameters. A machine learning framework is trained on simulated data, and then used to predict assembly model parameters for 1000 simulations generated with known parameter values sampled from a prior distribution. Cross-validation results are depicted indicating the true parameter value on the x-axis and the predicted value on the y-axis. Parameters evaluated in this figure include ecological strength, community equilibrium state, colonization rate, and trend (ClaDS α) and stochasticity (ClaDS σ) parameters of the macroevolutionary process, as well as the relative change in speciation rate (ClaDS m). This cross-validation procedure was performed using only raw nucleotide diversity per species in the local community as the test data. Please note that R 2 values in the figure were calculated for the exact cross-validation fold depicted, so these will differ slightly from the average R 2 over all cross-validation folds as shown in Table1.

Figure 4 -

 4 Figure 4 -Cyprus microarthropod metabarcoding data and inference results a) The geographic distribution of sampling sites within the Troodos mountain range. Each point is a sampling site, and sites are colored according to the primary forest type. b) The observed rank-ordered distributions of nucleotide diversity (π) combining all sites within each forest type. Nucleotide diversity per OTU within each forest type is used for machine learning assembly model classification and parameter estimation. c) Prior predictive check to assess whether model parameter values are appropriate for generating simulations similar to the observed data. Simulated (in gray) and observed data (points colored as in sub-panel b) are projected into principal component space and the first two principal components are depicted. Overlap between the simulated and observed data indicates that the priors are appropriate. d) Predicted community assembly model class probability for each forest type using our inference approach. The fraction of the bar indicates model classification probability for each of the three community assembly models: competition (black), environmental filtering (dark gray), and neutral (light gray). Forest type codes are as in sub-panel b.

Acknowledgements

This work is a product of the iBioGen project, which has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No. 810729. V.N. was supported by a postdoctoral contract under the iBioGen project and a "Juan de la Cierva-Formación" postdoctoral fellowship (grant FJC2018-035611-I) funded by MCIN/AEI/10.13039/501100011033.

Data availability

The iBioGen package is available as both R and python implementations. Source code and jupyter notebooks sufficient to reproduce all simulations and analyses are publicly available on the iBioGen github page: https://github.com/iBioGen/iBioGen. Empirical data from Noguerales et al. (2021) are available from the NCBI Sequence Read Archive under BioProject number PRJNA770337.

Author contributions

IO and HM designed the model. IO coded the model and performed all simulations and analyses. IO drafted the manuscript with support from AP, VN, EM, and HM. All authors contributed substantially to manuscript revisions.

Tables and Figures with captions

Table 1 -Quantifying machine learning performance using iBioGen model simulations Machine learning classification and regression performance were evaluated by taking the average of scores across a 5-fold cross-validation procedure. Classification performance was quantified with average precision, recall and F1. All these measures are bounded between 0 and 1, with 1 indicating perfect performance. Regression performance was quantified with the coefficient of determination (R 2 ), which has an upper bound of 1, in the case of perfect performance. The computational definition of R 2 we used is unbounded on the lower end, with negative values indicating performance worse than random. Results shown are for gradient boosting methods.