Community structure of vascular epiphytes: a neutral perspective

Vascular epiphytes form a diverse group of almost 30 000 species, yet theory concerning their community structure is still largely lacking. We therefore employed the simplest models of biodiversity, (near-)neutral models, to generate hypotheses concerning their community structure. With recently developed tools for (near-)neutral models we analyzed species abundance data from many samples in Central and South America which we divided into four metacommunities (Mesoamerica, Central America, Amazonia and Parana), where for each metacommunity we considered two subsets differing in dispersal syndrome: an animal-dispersed guild and a wind-dispersed guild. We considered three models differing in the underlying speciation mode. Across all metacommunities, we found observed patterns to be indistinguishable from patterns generated by neutral or near-neutral processes. Furthermore, we found that subdivision in different dispersal guilds was often supported, with recruitment limitation being stronger for animal-dispersed species than for wind-dispersed species. This is the first time that (near-)neutral theory has been applied to epiphyte communities. Future efforts with additional data sets and more refined models are expected to further improve our understanding of community structure in epiphytes and will have to test the generality of our findings.


Introduction
Ecologists have been interested in the distribution of vastly different abundances of species for decades. In a review of species abundance distributions (SADs), McGill et al. (2007) pointed out that an SAD contains crucial information about the forces that shape ecological communities. However, they argued that there are many theories predicting roughly the same SAD shape; for inferences to be made about underlying mechanisms, we should look beyond the single-sample SAD. McGill et al. (2007) also suggested that partitioning, i.e. taking subsets of, the SAD may provide insights into differences in community structure. Here, we follow up on both suggestions by studying the SADs of multiple community samples of vascular epiphytes, where we further partition the data into two guilds of species: a guild consisting of animal-dispersed Community structure of vascular epiphytes: a neutral perspective Thijs Janzen, Gerhard Zotz and Rampal S. Etienne 2 species, and a guild consisting of wind-dispersed species. These two dispersal modes have been shown repeatedly to lead to different spatial patterns of aggregation (Seidler andPlotkin 2006, Li et al. 2009).
Vascular epiphytes, a diverse group of almost 30 000 species with mostly tropical distributions (Zotz 2013a), are especially suitable for this approach, given high functional equivalence amongst species (Catchpole and Kirkpatrick 2011). Because theory on the mechanisms underlying community structure in vascular epiphytes has only recently started to emerge (Kitching 2006, Burns 2007, Burns and Zotz 2010, Mendieta-Leiva and Zotz 2015, Taylor et al. 2016, our current understanding of the structure of vascular epiphyte communities is mainly based on a few deterministic axes. For example, at a larger geographical scale, epiphytes are known to be more responsive to moisture availability than any other life-form (such as trees, shrubs or herbs), with the highest epiphyte diversity found in wet forests with low seasonality (Gentry andDodson 1987, Kreft et al. 2004). At a local scale, establishment, growth and survival of epiphytic plants are affected by several factors, with their relative importance still being debated. Primary factors are vertical abiotic gradients (Griffiths and Smith 1983, Hietz and Briones 1998, Petter et al. 2016, host tree identity (Laube andZotz 2006, Wagner et al. 2015), and host tree size and age (Flores-Palacios andGarcía-Franco 2006, Taylor et al. 2016). Even taken together these factors still hardly suffice to explain the many aspects of community structure in epiphytes, such as frequently observed extreme local patchiness or very high local or point diversity (Catchpole and Kirkpatrick 2011), which are indicative of a high degree of ecological equivalence among many species. Hence, stochasticity (e.g. neutral processes) may be highly relevant to understand community structure and dynamics in vascular epiphytes, but also for the linkage among communities, i.e. as one determinant of β-diversity.
Against this background, neutral theory seems to be an ideal starting point for studying community structure in vascular epiphytes. The spatially implicit neutral model (Hubbell 2001) for multiple local communities (Etienne , 2009a has a parameter θ that measures the regional diversity due to the biogeographical processes of speciation and extinction. Furthermore, it uses parameters I i for each local community sample i that measure the degree of recruitment limitation (i.e. both dispersal limitation and establishment limitation (Jabot et al. 2008)) of immigrants arriving at these local communities. Fitting the neutral model to the SADs of multiple samples of vascular epiphytes, either partitioned into dispersal guilds or not, will therefore inform us on any differences in recruitment limitation between sites and between dispersal guilds (Janzen et al. 2015). Hence our goal is not so much to test the neutral theory, but rather to use it as a starting point to study the community structure of vascular epiphytes.
Using the neutral theory, we postulate a number of hypotheses regarding vascular-epiphyte communities. Firstly, given the indicators of ecological equivalence amongst species, we expect for many communities to find patterns that are indistinguishable from neutral patterns. Yet, regional differences are likely, and some communities may show deviations from the neutral model, suggesting other processes driving local community assembly as well, such as competition or environmental factors benefitting specific species or species groups (e.g. tank bromeliads (Gilmartin 1983)). Furthermore, local environmental conditions, such as flooding regime, might favor specific hosts, thus cascading differences in associated epiphyte flora. Therefore, using the neutral model will be a first step in disentangling these effects from patterns caused by neutral processes (e.g. stochasticity). Secondly, we expect dispersal limitation to vary across metacommunities due to variation in topography, where we expect the flat Amazonian forests to be less dispersal limited than the mountainous regions in Central America or the Andes (Pérez-Escobar et al. 2017). Furthermore, because species relying on wind dispersal (e.g. orchids, ferns) are more prevalent in Amazonia, we expect dispersal limitation to be lower in Amazonia, where the small seeds of wind-dispersers can travel long distances. Thirdly, we expect to find support for a subdivision into dispersal associated guilds and expect that species relying on wind dispersal are less dispersal limited. Lastly, given evidence of potential hyperdominance in epiphytes (Zotz 2007), we expect that the diversity-dependent speciation model will, at least for some datasets, provide a better explanation of the data than the standard neutral model.
Here we analyze all available vascular-epiphyte SAD data sets from moist lowland forests that have absolute abundances for all species and total samples sizes of at least a few hundred individuals as a prerequisite to obtain reasonably accurate estimates of the neutral model parameters. We first describe these data sets, and briefly review the neutral model and the ingredients that are important for our analyses. Then we report the fit of the neutral model to the data, associated parameter estimates, and draw first conclusions on local and regional differences and guild differences in recruitment limitation and speciation rates.

Data sets
The term 'epiphyte' is inconsistently used in the literature (discussed by Zotz 2016). Apart from plants that germinate and establish on trees without ever making contact with the forest soil (true epiphytes), some studies also include hemiepiphytes (which establish as epiphytes, but later establish contact with the forest soil) or nomadic vines (which germinate on the ground and may or may not become epiphytic later in their ontogeny) as 'epiphytes'. As argued by Zotz (2013b), nomadic vines diverge strongly in their ecology from true epiphytes and should not be merged with them. Consequently, we eliminated nomadic vines from all data sets in our study. Hemiepiphytes, on the other hand, share the vulnerable  juvenile stage with true epiphytes and were therefore considered in our analysis. However, we ran two separate analyses, one in which both epiphytes and hemiepiphytes were included and one with only epiphytes. Data sets were either available in the original publications (Table 1) or supplied by the authors (see Acknowledgements). With one exception, all censuses were plot-based, although both plot numbers and sizes varied substantially (details in Table 1). Based on the location of the sampling sites, we grouped sites close to each other into metacommunities. The resulting four metacommunities are referred to as Mesoamerica, Central America, Amazonia and Paraná. Summary statistics and locations of the sites are given in Table  1, Fig. 1, respectively. There are two data sets for Mesoamerica, from a so-called 'Tintal' wetland area and from a mature, seasonally dry forest on Yucatán Peninsula with more than 3000 individual epiphytes (Goode and Allen 2008). Two censuses are available, before and after a hurricane, which are denoted by Mesoamerica pre-hurricane and Mesoamerica post-hurricane. We studied both censuses because we anticipated that community assembly might differ after a hurricane.
There are two data sets for Central America. The first (San Lorenzo) encompasses more than 13 000 individuals in an area of 0.4 ha of lowland rain forest at the San Lorenzo Canopy Crane in Panama (Zotz and Schultz 2008). The second data set (Río Changuinola) stems from a census of around 0.1 ha of undisturbed forest (500 m a.s.l.) near Río Changuinola (Bocas del Toro Province, Panama), yielding almost 9000 individuals (Wester et al. 2011).
Eight epiphyte datasets are from forests in Amazonia. The first data set (Surumoni) represents 1.5 ha of moist lowland forest in the Surumoni region (Venezuela) with about 1000 individuals (Nieder et al. 2000, Schmit-Neuerburg 2002. In this case data from two consecutive censuses are available (1996 and 2000). The second data set (Tiputini) has over 8000 individuals of epiphytes in 0.1 ha of lowland forest near the Tiputini Biodiversity Station in Ecuador (Kreft et al. 2004). The third data set (Mocagua) represents an epiphyte community on an island in the Amazon River in the southernmost tip of Colombia and contains almost 800 individuals recorded in several small plots in different forest types (Higuera Díaz 2003). The fourth set (Caquetá) represents epiphyte communities from four different forest types (total area 0.5 ha) in the middle Caquetá area (Colombia) with almost 4000 individuals (Benavides et al. 2005). The fifth data set (Adolpho Ducke) encompasses more than 20 000 individuals of vascular epiphytes representing 82 species from lowland forest (Boelter et al. 2014). A sixth data set comes from Amacayacu National Park, Colombia with close to 700 epiphytic plants representing 54 species (Benavides et al. 2006). Data set 7 represents 2413 individual epiphytes from a terra firme forest near Coari in the Urucu River basin (Irume et al. 2013). The last data set from Amazonia (Quaresma et al. 2017) represents almost 3000 epiphytes from two forests, a 'várzea' forest at the Mamiraua Reserve and an 'igapó' forest in Parna (Jau National Park).
A final data set used in our analysis was published by Geraldino et al. (2010). They report about 700 epiphytes in 37 species from a 30 ha transition area between semi-deciduous forest and Araucaria forest near Campo Mourão, Paraná.
We compiled all data sets in a single data file. Species names were standardised against The Plant List using the R package Taxonstand (Cayuela et al. 2012). All species were compared with an updated version of the global checklist of vascular epiphytes (Zotz 2013a(Zotz , 2016 and accidental epiphytes and nomadic vines were removed from the list before the analysis.

Tree community data
To compare community assembly between trees and epiphytes, we compared tree community data from the sites in the Central American metacommunity with the abovementioned epiphyte community data. Tree abundance data was collected in the same area, although the overall plot size was larger (Condit et al. 2002, 2004, Daguerre andCondit unpubl.), including trees not sampled for epiphyte abundance. Only trees with a diameter at breast height > 10 cm were included, and resulting sample sizes were 515 individuals in Río Changuinola (1 ha plot) and 3363 individuals at the San Lorenzo site (6-ha plot). Associated species richness was 139 and 129 species, respectively.  (7), one in Amacayacu National Park (9), one in Adolpho Ducke Reserve (10), one in Coari (11), two in Parna Miraraua (12/13). Paraná (diamond) with one sample: Campo Mourão (14).

Models
The most widely used spatially implicit neutral model of biodiversity (Hubbell 2001, Volkov et al. 2003, Etienne 2005, Rosindell et al. 2011 assumes that individuals of all species in a local community are equally subject to stochastic birth and death events, while individuals of new species can arrive through immigration from the larger metacommunity that is governed by speciation and extinction processes. This metacommunity is summarized by the fundamental biodiversity parameter, which is a compound of the speciation probability per birth ν and the metacommunity size J M : Immigration into local community i, which includes both dispersal and establishment, is characterized by the fundamental immigration number I i (Etienne and Alonso 2005): where m i is the probability of immigration into community i and J L is the local community size. We use I instead of m because I is relatively insensitive to sample size when estimated from SAD data . While early versions of this model only considered one local community, the model concept has subsequently been extended to include multiple samples , Forster and Warton 2007, Munoz et al. 2007, 2009b. A full joint likelihood of model parameters based on the distribution of species abundances across these samples, keeping track of species identity, was provided in  and made practically applicable in Etienne (2009b) without restricting assumptions that were still present in . Differences in dispersal between groups of species that share the same life-history might impact estimates as well. A full sampling formula estimating dispersal separately for different dispersal-guilds was proposed by Janzen et al. (2015), but this formula was not yet able to accommodate multiple samples at the same time. The joint likelihood for multiple samples and multiple guilds was developed by Haegeman and Etienne (2017), and this is the sampling formula that we use here. Note that it assumes that I may differ between guilds and between sites, θ is fixed across guilds. We apply this likelihood formula (called a sampling formula) to the epiphyte data. For details of this formula we refer to the above-mentioned papers. The standard neutral model applies a point-mutation speciation model in the metacommunity, where novel species are introduced as singletons at a fixed rate. Previous studies have shown that this speciation model might not always be appropriate Haegeman 2010, Rosindell et al. 2010). Therefore, we also applied two alternative speciation models and compared fits across models.
Firstly, we applied the protracted speciation model, which assumes that species first go through an incipient stage before speciation is completed (Rosindell et al. 2010). Thus, speciation in the protracted speciation model takes time, in contrast to the instantaneous process reflected by the point-mutation speciation model. The protracted speciation model uses an additional parameter (φ) that indicates the expected rate of completion, which is the inverse of the expected duration of the incipient stage (Haegeman and Etienne 2017). Hence, when φ = ∞, speciation is instantaneous and the protracted speciation model is identical to the point-mutation speciation model.
Secondly, we applied a density-dependent speciation model, where speciation is not constant per individual (Etienne et al. 2007, Haegeman andEtienne 2009). Whereas in the point mutation speciation model, species that are more abundant in the metacommunity are more likely to speciate, the density-dependent speciation model modifies this effect with an additional parameter α (Haegeman and Etienne 2017), that reflects the degree of density-dependence. For α = 0 the model reduces to the standard neutral model with constant per-individual point mutation. For α > 0, positive density-dependence reduces the probability of speciation for abundant species, where for α = 1 the model reduces to a model with constant per-species speciation. For α < 0, negative density-dependence increases the probability of speciation for abundant species, where for α = −1 the model reduces to the random fission model of speciation (Hubbell 2001, Etienne and Haegeman 2010, Haegeman and Etienne 2017. While originally introduced as density-dependence 'speciation', the model can also be interpreted as reflecting density-dependent 'birth' in the metacommunity (Haegeman and Etienne 2017). The effect is the same: positive densitydependence makes abundant species more likely to occur, and negative density-dependence makes abundant species less likely to occur. This model is not strictly neutral for two reasons: firstly because the speciation rate applied to individuals can be weakly sensitive to species identity resulting in density dependence and secondly because there can be multiple guilds of individuals with different dispersal abilities. Within each guild, this is nevertheless a symmetric model where the species labels have no biological meaning. Because of these reasons, we refer to this model as a near-neutral model.
We also explore extensions of these models to models where within each site, two guilds of species coexist, that only differ in their dispersal ability. Although both guilds exist separately in the metacommunity, they share the same speciation dynamics.
In summary, we explore three different speciation models in the metacommunity (point-mutation, protracted-speciation and density-dependence), whilst keeping local community dynamics identical. The model with density-dependent speciation involves departures from the strict neutrality assumption, and are referred to here as near-neutral models (given that all other aspects of these models closely resemble the standard neutral model). For each model we explore an extension with subdivision into two dispersal guilds.

6
We used the maximum likelihood framework as implemented in the R package SADISA (Haegeman and Etienne 2017). This framework optimizes the likelihood for the parameters given the data. Because maximum likelihood optimization can be sensitive to the starting conditions, we used three different starting conditions, and report here the optima with the highest likelihood across these three starting conditions. The starting conditions used were [θ, I] = [10, 1000], [1000, 10] and [500,500]. Previous applications of the neutral model of biodiversity and biogeography have found that there can be multiple optima of comparable likelihood, one optimum with low θ and high I and one optimum with high θ and low I (Etienne and Alonso 2006, Jabot and Chave 2009. To ensure we recover the global optimum, we chose two sets of initial parameter values that are biased towards either of these two potential optima, and one additional initial parameter set in between these two first sets. For the protracted speciation model two initial values of φ were explored, one large value (1000) and one very small value (0.001). These two values represent either a model where speciation is almost instantaneous (for large values the model reduces to the point mutation model), or where speciation is extremely protracted (for small values). This brings the number of initial parameter sets to 6 for the protracted speciation model. For the density-dependent model two values of α were explored, one negative (-0.5) and one positive (0.5). This also brings the number of initial parameter sets for the density-dependent model to 6. An overview of all starting values can be found in the Supplementary material Appendix 1 Table A1.
We divided the species into two guilds: wind-and animal-dispersed species for each of the metacommunities. In many cases, species were assigned to either group using basic biological information at the family level, e.g. all orchids and ferns were considered wind-dispersed although very rare exceptions are known for these groups (Benzing and Clements 1991), while, e.g. aroids or Piperaceae were invariably considered to be animal-dispersed. In other groups, e.g. bromeliads, dispersal syndromes were deduced from a range of sources, mostly notes in the taxonomic literature. The status of all species (wind-or animal-dispersed) is reported in the data file.

Test of model fit
In line with the 'exact' test of neutrality of Etienne (2007), we tested whether the best fitting model is a good description of the data, because, even though a model may fit the observed data best out of the all the fitted models, the observed data may still deviate substantially from this model. To test for this, we simulated 100 datasets using the best fitting model (using the maximum likelihood parameter estimates), and then again performed maximum likelihood on the simulated datasets. If the observed data were a product of the same process as the simulated data, we expected the maximum likelihood of the best fitting model to the observed data to lie within the distribution of maximum likelihoods obtained for the 100 simulated datasets. To obtain a p-value we counted the number of datasets for which the likelihood was lower than the likelihood of the empirical data. Datasets were simulated using the R package SADISA (Haegeman and Etienne 2017).
As a second test for goodness-of-fit, we studied whether observed community compositions (rather than overall abundances) within sites reflect the best fitting models' distribution. Here we based the p-value on the distribution of Bray-Curtis dissimilarities, as described in Aduse-Poku et al. (2018). To do so, we first calculated the average Bray-Curtis dissimilarity (BC) between a focal site and all other sites in the metacommunity. Then, for each simulated dataset, we also calculated the average BC between the focal (simulated) site and all other sites in the (simulated) metacommunity. This yielded 100 average BC values (one for each simulation replicate). We then compared the average BC in the empirical data with the distribution of BC values across the simulations. The BC p-value is then the fraction of simulated datasets that has a lower average BC than the average BC recovered in the empirical data. The advantage of using BC is that it allows for studying individual sites, whereas the traditional exact test of neutrality based on likelihood distributions can only test the entire metacommunity as a whole.

Error in parameter estimates
To obtain an uncertainty estimate for our parameter estimates, we used the 100 simulated neutral datasets created for the tests of goodness-of-fit, and calculated the 95% percentile of the distribution of inferred parameter values.

Model selection
Within each family of community models we studied three models for each community: 1) a 'free' model where all parameters were inferred independently, 2) a model where we fixed dispersal across guilds, such that species on different sites belonging to the same guild shared the same dispersal rate (we called this 'Fix guilds') and 3) a model where we fixed dispersal across sites, such that species on the same site, regardless of guild, shared the same dispersal rate. Using the maximum likelihood estimates, we computed the Akaike information criterion (AIC) for each model, where the AIC is given by AIC = 2k − 2ln(L) (Akaike 1974), where k is the degrees of freedom and L is the likelihood of the model, given the data. Degrees of freedom for the different models were: 1) for the 'Free' model: 1 + (the number of sites times the number of guilds), i.e. one θ and all the corresponding I values, one for each site-guild combination, 2) for the 'Fix guilds' model: 3, i.e. one θ + two I values, one for each guild and 3) for the 'Fix sites' model: 1 + the number of sites, i.e. one θ and one I value for each site. The protracted-speciation model and the density-dependent model each have one more parameter. AIC values were compared with each other via AIC weights (Wagenmakers and Farrell 2004), which converts the differences between AIC values into the relative probability of a model being the best model out of the tested models.

Influence of undetermined species
For some specimens it is difficult, or even impossible, to determine it to species or even genus level. Across all datasets, we found many samples labeled as 'Genus sp.' or 'unidentified' (4661 samples out of 67 530 samples, Supplementary material Appendix 1 Table A4). Although we are confident that specimens were consistently labeled to their respective unidentified label within each sampling effort, when comparing datasets it becomes impossible to verify whether specimens labeled in multiple datasets as 'Genus sp.' actually belong to the same species, which might inflate the overall species diversity in our communities. Therefore, we repeated our analysis removing all undetermined species. For the Mesoamerican and Paraná datasets undetermined species that were shared with other datasets were not found (these datasets each stem from a single paper, see also Supplementary material Appendix 1 Table A4).

Comparison of dispersal rates
Local per-community differences in dispersal across guilds are reflected by a higher likelihood of the model with guild structure. However, to assess whether there were also differences in dispersal across communities and metacommunities, we also performed an ANOVA to test for differences in dispersal. We started out with a full model with three factors: guild, metacommunity and site. We then removed components in a stepwise fashion to obtain a model only containing significant components. To compare dispersal estimates across datasets (e.g. datasets 'Epiphyte' and 'Epiphyte + Hemiepiphyte'), we performed a similar analysis, including dataset as an additional factor.

Influence of sampling area
To assess whether sampling area influenced our estimates of dispersal, we performed two one-way ANOVAs, with dispersal or θ as response variable and size of the sampling area in hectares as independent variable.

Results
Across the studied metacommunities (Central America, Amazonia, Paraná and Mesoamerica), the density-dependent model provided the best fit for two metacommunities (Central America, Amazonia) (Table 2, Fig. 2). For the remaining metacommunities, the point mutation speciation model provided the best fit. Across all metacommunities, the highest AIC weights were associated with the 'free' models, in which each site and each guild has a separately inferred dispersal rate. Only for the Mesoamerican datasets we found models with dispersal rates fixed across sites and guilds to perform slightly better (Table 2).
Dispersal rate was higher for wind-dispersed species (p < 0.001 for epiphytes, two-way ANOVA, Table 3, Fig. 3, also evidenced by the fact that the model with dispersal fixed across guilds had a lower AIC weight), except for the Amacayacu and Coari dataset. It is unclear why these two sites within the Amazonia metacommunity showed opposite patterns; local species diversity and numbers of individuals are very similar to other sites in the Amazonia metacommunity. Inclusion of hemiepiphytes did not change this pattern (p = 0.006, two-way ANOVA), and we did not find significant differences between the two datasets (two-way ANOVA, p = 0.59, Fig. 3). Inclusion of hemiepiphytes increased the estimate of θ (no statistical testing performed, too few values), which is in line with expectations because the inclusion of hemiepiphytes increases both the number of species and the number of individuals. For the three metacommunities, where the density-dependent speciation model best fit the data, α values were positive.
Goodness-of-fit tests yielded p-values well above the threshold of 0.05, suggesting no evidence to reject the best fitting model as a good explanation for our data (Table 3, Fig. 4). Furthermore, when we performed the more detailed Bray-Curtis dissimilarity test, we still did not find any site that was significantly different from expectations (Table 3, Fig. 5) (but see Surumoni, animal-dispersed guild: p = 0.06). Including hemiepiphytes did not change this pattern. Lastly, using the 1996 data of the Surumoni plot did not change this pattern either (Supplementary material Appendix  1 Table A2, A3, p = 0.48 for epiphytes and p = 0.49 for epiphytes + hemiepiphytes).
Uncertainty in parameter estimates was high for the Mesoamerican datasets, particularly for dispersal estimates after the hurricane. This is likely due to the low number of individuals in the dataset. For the Central American datasets we only recovered high variation in the θ estimate, but not in the I estimate. This was potentially driven by a small number of outliers among the simulated datasets for which very high θ values were estimated, both for the epiphyte dataset and the dataset including hemiepiphytes. For the Amazonia dataset we observed reasonably low variation for estimates of θ, α and I, reflecting the high quality of this dataset. For the Paraná metacommunity uncertainty was high in all parameter estimates. Although the Paraná dataset contains many sampled individuals, this dataset represents only a single plot, resulting in the observed uncertainties.
The size of the sampled area and the estimate of θ were unrelated (ANOVA, p = 0.78), but the dispersal parameter I depended weakly on the sampling area (ANOVA, p = 0.04, intercept = 8.8, slope = 0.175, adjusted R 2 = 0.04). This weak relationship suggests that larger areas tend to be associated with higher rates of dispersal, although the explanatory power is low (R 2 = 0.04).
When we excluded undetermined species in order to correct possible inflation of species diversity due to the inclusion of undetermined species, we found very similar results (Supplementary material Appendix 1 Table  A5) and parameter estimates changed very little, i.e. the inclusion of undetermined species hardly affected our estimates (Supplementary material Appendix 1 Table  A6). Furthermore, when excluding undetermined species we also found for both the Central American and the Amazonia datasets that we cannot reject the neutral model, both based on likelihoods and based on Bray-Curtis dissimilarities (Supplementary material Appendix  1 Table A6). This further emphasizes that the inclusion of undetermined species does not affect our findings, in line with previous findings (Pos et al. 2014).
Results for trees at the Central American sites were very similar (Table 4, 5). For these data, the density-dependent model provided the best fit to the data (identical to the epiphyte data fits), and we also found support for positive density-dependence (Table 4). These findings of positive density-dependence are in line with previous parameter estimates for six communities of tropical trees (Haegeman and Etienne 2017). Also in line with our I estimates for epiphytes and epiphytes + hemiepiphytes, estimates of I for the Río Table 2. Comparison of likelihoods and AIC weights for all models and for all datasets. Indicated are the models (pm = point mutation, pr = protracted speciation, dd = density dependent), the linkage between parameters (free = all parameters independent, 'sites' = same dispersal within sites across guilds, 'guilds' = same dispersal within guilds across sites). Shown are the degrees of freedom (df), the likelihood (LL), the resulting AIC, the ΔAIC and finally the AIC weight. Models with the highest AIC weights are highlighted in bold. When AIC weights differed by < 0.5%, both models are highlighted. Changuinola site were higher than for the San Lorenzo site (Table 5). I estimates for trees are much higher than for epiphytes (Table 3, 5), suggesting higher dispersal limitation for epiphytes compared to trees. Model fit tests for the tree data cannot reject the best fitting model as a good description of the data, neither using the likelihood nor using Bray-Curtis dissimilarity (Table 5).

Discussion
This study was initiated to use neutral theory to study the community structure of vascular epiphytes, and more specifically to 1) identify possible regional differences in diversity and the integration of local communities in a larger metacommunity, and 2) to detect possible differences among Figure 2. Rank abundance plots (black) for the Epiphytes dataset (excluding hemi epiphytes) for each site, separated into animal (gold) and wind (blue) dispersed guilds. Shown are results for the five different metacommunities (Table 1). Blue lines indicate the average abundance across 100 replicate simulations, using the best fitting model. Shaded areas indicate the 95% confidence intervals. epiphytes within a region in respect to recruitment limitation (animal based vs. wind based dispersal). Recruitment limitation may be expected to be a function of topography in these lowland forests with highest I-values in the vast Amazonian forests and lower I-values in Central America. However, our results are inconclusive, e.g. I-values in Amazonia are sometimes higher, sometimes lower than in Central America. This finding may simply be the consequence of the low number of available data sets or of differences in sample area (larger areas are expected to have lower I-values ), which make any conclusion preliminary, but may well reflect some unexpected biological reality.
A highly noteworthy finding is the higher recruitment limitation in epiphytes compared to trees for the Central America metacommunity. This is unexpected because in contrast to trees, the vast majority of epiphytes are anemochorous (Zotz 2016), which should translate into lower recruitment limitation. Yet, based on biogeographical approaches it has been noted that epiphytes can differ systematically from ground-rooted plants in their range, which may be directly related to differences in dispersal ability among life forms. In most studies, epiphytes were found to show a lower degree of endemism and larger average range compared to terrestrial plants (Ibisch et al. 1996, Nieder et al. 1999, Kessler 2000, but other studies report the opposite (Kreft et al. 2004). Most of the cited studies were conducted in montane regions or with entire national floras in countries with a large proportion of montane regions, while the study of Kreft et al. (2004) was conducted in lowland Amazonia. Interestingly, our results are from lowland sites as well, and it is tempting to speculate that relative changes in the ranges of vascular epiphytes versus terrestrial plants are a function of elevation. Consistent with this notion, Kluge et al. (2008) report more narrow elevational ranges of epiphytic ferns in the lowlands compared to higher altitudes in Costa Rica. Conversely, the mean elevational ranges of terrestrial species in that study showed a sharp decline at high elevation, while those of epiphytes remained rather constant. Clearly, further studies in Figure 3. Distribution of inferred dispersal rates for both dispersal guilds, shown for the dataset including only epiphytes, and for the dataset covering epiphytes and hemiepiphytes. Wind-dispersed species have a significantly higher dispersal rate than animal-dispersed species, whereas the inclusion of hemiepiphytes did not lead to significant differences. Figure 4. Log-likelihood of the Maximum likelihood estimate for 100 datasets simulated using the best-fitted model ( Table 3). The blue arrow indicates the log-likelihood value of the empirical data. The p-value reported in Table 4 is the fraction of simulated datasets with a likelihood smaller than the blue arrow. Shown are results for the datasets without hemiepiphytes. other regions are needed before more general conclusions are possible.
Many of the results of this paper are consistent with previous observations. For example, compared to Central America, orchids are known to be less dominant in many Amazonian epiphyte communities, while aroids become more important (Leimbeck andBalslev 2001, Benavides et al. 2005). This shift in dominance automatically leads to changes in the Figure 5. Distribution of Bray-Curtis (BC) dissimilarities for 100 datasets simulated per metacommunity, using the best fitted model ( Table  3). The blue arrow indicates the average BC dissimilarity for that specific site and guild in the empirical data. The BC p-value in Table 4 then reflects the fraction of simulated datasets with a BC dissimilarity smaller than the blue arrow. Shown are results for the dataset without hemiepiphytes.
13 relative diversity of wind-and animal-dispersed taxa in the metacommunity (Table 1). We expected that the dust-like seeds of orchids, the minute spores of ferns, and the winged or plumed diaspores of other wind-dispersed taxa, e.g. many bromeliads, would translate in more effective long-distance dispersal in Amazonia. This is exactly the pattern that we observe, i.e. values of I for wind-dispersed species are higher in Central America than in Amazonia.
Arguably, treating local epiphyte communities as units of investigation still represents a rather coarse scale because any epiphyte community is characterized by an inherent β-diversity due to the vertical gradient in abiotic conditions. This gradient, which has been repeatedly shown to have rather deterministic consequences, results in a strong stratification of species by physiological and morphological traits related to drought or shade (Hietz and Briones 1998, Stuntz and Zotz 2001, Petter et al. 2016, and offers itself to define meaningful ecological scales in the analysis of epiphyte communities (Mendieta-Leiva and Zotz (2015). Unfortunately, quantitative information on this vertical structure is rarely available for entire communities. Assuming a sufficiently large sample, it should be possible to separate a local epiphyte community into several vertical zones or more refined guilds (Johansson 1974, Kernan and Fowler 1995, Zotz and Schultz 2008 and analyze whether we can identify groups of functionally equivalent species associated to specific vertical zones (neutral guilds), sensu Aduse-Poku et al. (2018).
For the Amazonian and Central American metacommunities the best fitting model is a density-dependent model with a positive α value. This indicates that abundant species are less likely to undergo speciation, or have a competitive advantage compared to rare species, which suggests hyperdominance of some species in the community. This is consistent with available data on the population structure of 45 species in an epiphyte community in lowland Panama (Zotz 2007), where more abundant species tend to have a higher proportion of juveniles. Interestingly, we also found positive density-dependence for the tree communities in Central America, which is in line with previous support for positive density-dependence in tropical tree communities (Haegeman and Etienne 2017).
Detecting deviations from the neutral model is a data intensive process. Here, we have presented several datasets of varying size. Across these datasets, we cannot detect deviations from the standard neutral model in three metacommunity datasets (both Mesoamerican datasets, and Parana). This implies that the effect size of non-neutral processes is too small to be statistically significant. This can be either due to a small effect size per se, or a small sample size, which are two sides of the same coin. We note that in the two other metacommunities, where we do detect deviations from neutrality (in the form of non-neutral speciation dynamics), datasets represent large sample sizes. Furthermore, AIC weights for the three datasets indicate some support for alternative models. Combined, this suggests that a lack of information content is driving our inability to detect deviations from the standard neutral model for these three metacommunities.
With this paper we have shown that vascular epiphyte communities in many different localities can be described and analysed with neutral, or near-neutral models. Given the still sketchy database on epiphyte communities and the heterogeneity in the structure and quality of available data we did not expect to find ultimate answers on epiphyte community structure. Rather, our study is one of several current efforts to produce testable ideas (Burns and Zotz 2010, Mendieta-Leiva and Zotz 2015, Taylor and Burns 2015, Spruch et al. 2019. The presented analysis highlights important properties of epiphyte communities, e.g. regional differences in recruitment limitation, differences between wind-and animal-dispersed taxa and positive density-dependence. Testing these still preliminary findings with additional data sets, and refining the models by separating 'communities' into more refined guilds (e.g. acknowledging vertical zonation) will be a task of the future.

Speculations
We found evidence for hyperdominance in epiphytes, i.e. some species have higher abundance than expected from strict neutral theory. The model underlying this conclusion Table 4. Comparison of loglikelihoods and AIC weights for all models applied to the tree dataset. Indicated are the models (pm = point mutation, pr = protracted speciation, dd = density dependent), the linkage between parameters (free = all parameters independent, 'sites' = same dispersal within sites). Shown are the degrees of freedom (df), the loglikelihood (LL), the resulting AIC, the ΔAIC and finally the AIC weight. The model with the highest AIC weight is highlighted in bold. assumes that either per species speciation rate decreases with abundance or that the birth rate in the metacommunity is subject to positive density-dependence; both assumptions lead to mathematically equivalent results. Accepting the second assumption suggests a minor role of parasites and/or seed predators in promoting epiphyte diversity in contrast to the Janzen-Connell hypothesis, which assumes negative densitydependence. This conclusion largely agrees with findings from a number of empirical studies with epiphytes (summarized in Zotz 2016) on the role of herbivores and frugivores. However, we also find hyperdominance for trees, which is in line with other assertions (Ter Steege et al. 2013). If parasites and seed predators are to play a substantial role, then the first assumption must be accepted: species are less likely to speciate when they are highly abundant. This in turn suggests that species become more plastic or generalist as they reach very high abundances, making them less prone to speciate.