Frontiers of Biogeography Species richness and composition of Caribbean aquatic entomofauna: role of climate, island area, and distance to mainland Frontiers

From a literature review, we constructed a database comprising >1000 freshwater insect species (especially Odonata, Coleoptera, Trichoptera, Ephemeroptera; OCTE) in 26 Geographical Caribbean Units (GCU) and quantified local filtering (climate heterogeneity, annual rainfall, annual temperature), geography (area, distance from the mainland) and emergence age as a proxy for island ontogeny. We investigated the relative role of these variables on the species richness, endemism and composition of the units using island species-area relationship (ISAR), generalised linear modelling (GLM) and nonmetric multidimensional scaling (NMDS). In addition, we analysed the spatial patterns of species richness and composition using Moran’s I index. ISAR generally demonstrated one or two thresholds and continuous or discontinuous responses according to OCTE groups. A small island effect could be detected for Trichoptera and Ephemeroptera richness, whereas Odonata and Coleoptera only demonstrated differences in slope between smaller and larger GCUs. Area, climate heterogeneity, maximal rainfall and distance from mainland were major drivers of species composition in GCUs, whereas local climate variables were of main importance for the endemism rate. Due to the potential complexity of the Caribbean island ontogeny, middle-stage islands had an expected higher freshwater invertebrate richness than younger ones but an unexpected lower richness compared to older islands. Finally, the degree of colonization of islands was linked to the dispersal ability of species, with Odonata and Coleoptera having larger distribution ranges than Trichoptera and Ephemeroptera, which were more restricted by their comparatively narrow ecological niches. The high endemism (>60%) found in the Caribbean freshwaters calls for more conservation effort in managing these highly threatened freshwater environments.


Introduction
Oceanic islands have inspired the development of several theories in ecology and evolution (Valente et al. 2014), and many scientists consider these regions as "laboratories of evolution" (Whittaker et al. 2008). One of the most influential theories is MacArthur and Wilson's (1967) theory of island biogeography (TIB). Based on the premise that island species richness depends on a dynamic balance between species colonization and extinction processes, TIB assumes that larger islands or islands close to the mainland should contain a greater number of species than small islands or more distant islands. However, this dynamic balance between immigration and extinction alone cannot explain the richness of insular communities, and local environmental filtering and biotic interactions are also seen as significant evolutionary drivers of island species richness (Lack 1976, Gros-Désormeaux 2012. In addition, not all species are ecologically equivalent, and their colonization of an island may be successful or unsuccessful depending on their ecological requirements and dispersal ability (Nekola and White 1999). Therefore, Lomolino (2000) proposed a biogeographic conceptual model, i.e., the Immigration-Extinction-Evolution model, which includes the degree of isolation and the area of ecosystems as well as the characteristics of immigration filters (e.g., wind direction and speed, current velocity), geology and climate. Whittaker et al. (2008) further developed the General Dynamic Model of Ocean Island Biogeography (GDM), which broadened and updated MacArthur & Wilson's TIB by integrating the processes of immigration, speciation and extinction in a time frame related to the dynamic, transient nature of volcanic islands (island ontogeny). Island ontogeny is very important for the evolution of the diversity of island habitats. After their complete emergence, islands are transformed by erosion processes, which causes a reduction in topographical complexity and habitat diversity affecting the rate of immigration, speciation and species extinction (Fernández-Palacios de Nascimento et al. 2011, Whittaker et al. 2007. As a result, lower richness should occur on the newly formed as well as on the sinking and eroding late stage islands, whereas higher richness is expected in the middle-stage islands (Kraemer et al. 2022, Whittaker et al. 2017, Whittaker et al. 2008. Likewise, endemism rate should decrease with the increase in extinction resulting from ontogenetic changes (Chen andHe 2009, Whittaker et al. 2007). Finally, below a certain size, island area becomes a poor predictor of species richness, whereas habitat (stochastic events, natural disturbances) may be more influential on species richness (Triantis et al. 2006). The difference in the island area-species richness relationship between small and larger islands is known as the small island effect (SIE, Lomolino andWeiser 2001, Matthews andRigal 2021).
Due to its geological history and geographic location, the Caribbean archipelago is an excellent study system to ask certain biogeographic questions. In this region, islands are the result of various geological events (e.g., sediment uplift, accretion prism and volcanic eruption; Cruse 2014) occurring at distinct geological time scales. They are located close to the American continent (Central, North and South America) and include a variety of landscapes associated with a wide range of island areas (from 13 to >110,000 km 2 ) and elevations (from 20 to 3087 m). This topographical variability, even within a single island, creates varied climates, which are favourable to species richness. In addition, it may create barriers that favour endemism. For example, the Caribbean biota presents exceptionally high biodiversity with an endemism rate of up to 95.4% in some taxonomic groups of some islands, such as amphibians in Cuba (Rivalta-Gonzalez et al. 2014). Biogeographic studies in the Caribbean archipelago are numerous but do not concern all taxa. Most studies so far have focused on vertebrates (Hedges et al. 1992, Tavares and Simmons, 2000, Dávalos 2004, Camargo et al. 2009, Alonso et al. 2011) and terrestrial arthropods (Hollocher and Williamson 1996, Oneal et al. 2010, McHugh et al. 2014, Crews and Esposito 2020. These studies have shown that nonflying vertebrates (mammals, reptiles and amphibians) have close relatives in South America, that freshwater fish and flying vertebrates (birds and bats) have strong ties to North and Central America (Hedges 1996), and that terrestrial arthropods originate mainly from South America (McHugh et al. 2014) and partly from Central and North America (Crews andEsposito 2020, Davies andBermingham 2002). Rodriguez-Silva and Schlupp (2021) demonstrated that after reaching the Caribbean archipelago, new colonizers experienced high speciation because of the new habitat conditions they had to adapt to, which led to significant endemism. The adaptive radiation of species in the Caribbean archipelago is probably best illustrated by Anolis lizards, with more than 150 species (Losos and Ricklefs 2009). In addition, Calisto spp., a Lepidoptera Nymphalidae endemic to the Caribbean archipelago, shows remarkable radiation, with ~52 species distributed across Hispaniola (36 species), Cuba (11), Puerto Rico (1), Anegada (1), Jamaica (1), and the Bahamas (2) (Aguila et al. 2017).
In contrast, little is known about the spatial distribution of aquatic invertebrates in the Caribbean. This is surprising since a majority of Caribbean islands include a variety of large lowland and mountain rivers and many lakes and wetlands (CEPF 2010). Despite their small size, some of these remote islands are recognized for their high species richness. For example, the Virgin and Cayman Islands are exhibiting a particularly high level of biodiversity/endemism (CEPF 2010). A majority of biogeographical research on Caribbean aquatic invertebrates consists of scattered taxonomic studies providing species checklists. The only study dealing with aquatic insect biogeography in the Caribbean islands was by Flint (1978) and was limited to Trichoptera and Odonata. The author showed that Trichoptera species were very highly endemic, with 80% of the species being known from only one island, whereas the Odonata species were much more widespread, with only 15% of the species restricted to a single island. It should be noted, however, that this study was restricted to the Greater Antilles (Cuba, Hispaniola, Jamaica and Puerto Rico) and the mediumsized islands (Guadeloupe, Martinique, Grenada, Dominique and St Vincent). Flint (1978) (Bonada et al. 2012). This variation in species dispersal among taxa may partly explain the large-scale distribution of species and constitute a source of variation in the taxonomic composition and richness among Caribbean island aquatic communities (Bonada et al. 2012, Marquèz andKolosa 2013).
The main aim of this study was to quantify the relative contributions of geographical situation (area, distance from the mainland) and time of emergence of the islands, local environmental conditions (elevation, air temperature, rainfall) and dispersal ability of taxa in determining the overall species richness, endemic species richness and species composition of the aquatic entomofauna of the Caribbean islands. We achieved this based on a database of 1027 aquatic insect species across 26 GCUs assembled from the literature. Based on island biogeographic models and theories discussed above, we used our data to test four predictions: (1) Following TIB, we expected island area and distance to island to be the main driver of species richness with a higher richness in larger islands and a lower richness for remote islands. In addition, due to the SIE, if speciation rate increases with island area, a transition zone (threshold) should occur with a steeper relationship for larger than for smaller islands. (2) Due to ontogenetic changes (see GDM), we expected older and younger islands (higher rate of extinction) to exhibit lower species richness than middle-stage islands. (3) The size of the distribution ranges of species should reflect their varying dispersal abilities, with good dispersers (e.g. Odonata) having larger spatial distribution than weaker dispersers (e.g. Ephemeroptera). (4) Local environmental variables should secondarily influence species richness, endemism and composition of the Caribbean islands.

The Caribbean archipelago
The Caribbean is a vast archipelago that extends between Florida and Venezuela and consists of approximately 7,000 islands and islets of varying sizes covering an area of 234,000 km 2 (Ridvan (2007) (Fig. 1). The Caribbean archipelago encompasses three groups of islands: (1) the Greater Antilles (Cuba, Hispaniola, Puerto Rico and Jamaica), located at the north end and which represent 90% of the total area (GA in Table 1); (2) the Lesser Antilles, located at the centre and southern end, include small and medium-size islands, which are subdivided into middle stage (early emerged) and younger (lately emerged) islands (noted LA a and LA b respectively in Table 1); and (3) the Cayman Islands, the Bahamas and the Turks and Caicos Islands (OT in Table 1). Note that Aruba, Bonaire, Curacao, Trinidad and Tobago are usually included in the Lesser Antilles only due to their size, but here they were included into the third group for the purpose of our study (OT in Table 1). Most Caribbean islands arose from volcanic activities, with the exceptions of Barbados (accretion prism), Bahamas (accumulation of carbonate marine sediment) and Trinidad (fragmentation of South American continent) (Cruse 2014), and emerged at distinct geological times (Graham 2003; Table 1). Relief greatly varies across the Caribbean archipelago, comprising low-altitude islands (e.g., Cayman Islands    with an elevation of 20 m maximum) and islands with peaks beyond 3000 m a.s.l. (e.g., Hispaniola, 3087 m; countrystudies.us; https://www.worldatlas.com). In addition, climatic conditions vary from arid tropical to humid tropical with temperatures ranging between 15°C and 35°C and rainfall ranging between 200 and 9000 mm per year (https://www.meteo.cw/; http:// en.climate-data.org).

Study sites
Since there is a large number of very small islands in the Caribbean archipelago, we combined groups of islands into GCU. For example, we considered Cuba and its 4,194 satellite islands and islets as one GCU. Likewise, the Bahamas, which comprises 700 secondary islands, was considered a single GCU. Finally, twin islands such as St Kitts and Nevis, St Vincent and Grenadines, and Antigua and Barbuda were considered together. In total, we retained 26 GCUs (omitting Anguilla and the Turks and Caicos Islands due to a lack of available data and not shown on Fig. 1).

Species inventories
We gathered species inventories on each island based on research works including doctoral theses and scientific papers, Proceedings and United Nations Convention on Biological Diversity (CBD) reports and checklists using "Caribbean aquatic insects" as keywords. Several online databases (e.g., https://www. odonatacentral.org; https://sites.google.com/site/ distributionaldatabase) or Web portal (e.g., https:// www.gbif.org; https://www.dutchcaribbeanspecies. org, especially for ARU, BON, CUR, SAB, STE and STM) were also considered to complete our database. Part of the gathered information was verified by a leading scientist (F. Meurgey, pers. com.). References that were used to gather data are available as Supplementary Information (Appendix S1). Ten orders of aquatic insects were initially considered: Coleoptera, Diptera, Ephemeroptera, Hemiptera, Hymenoptera, Lepidoptera, Megaloptera, Odonata, Plecoptera and Trichoptera. Due to insufficient information for all groups in all GCUs (Appendix S2), we retained only the four best documented orders, namely, Odonata, Coleoptera, Trichoptera and Ephemeroptera (OCTE), which also represented different dispersal abilities (see Bonada et al. 2012). The validity of the name of each species was checked using the Interagency Taxonomic Information System (ITIS) to avoid synonymous species names, and we checked whether each species lived in aquatic habitats with the Catalogue of Life (https:// www.catalogueoflife.org/) and the Dutch Caribbean Species Register (https://www.dutchcaribbeanspecies. org).

Geographical and local variables
We documented two variables associated with TIB: unit area (AREA in Table 1, Meditz and Hanratty 1987; https://www.countrystudies.us; https://www. worldatlas.com) and the minimum distance from the mainland (DFM in Table 1). This value was obtained from the measurement of the shortest approximate distance of each unit from South America, North America and Central America using Google Earth™. For each GCU, we also documented the maximum elevation (ELEV in Table 1; Meditz and Hanratty 1987; https:// www.countrystudies.us; https://www.worldatlas. com), climate heterogeneity by counting the number (from 1 to 5) of different climates occurring in each GCU (CLIM in Table 1, see Appendix S3 for references), annual minimum and maximum air temperature (TMin and TMax, respectively, in Table 1), annual minimum and maximum rainfall (RMin and RMax, respectively, in Table 1), number of permanent streams (PERM in Table 1; see Appendix S3 for references), and the approximate period of emergence of units (AGE in Table 1; see Appendix S3 for references).

Data analyses
We investigated the island species area-relationships (ISAR) using piecewise regression models (Matthews and Rigal 2021). Seven models were considered: continuous one-and two-threshold, discontinuous one-and two-threshold, left-horizontal one-and two-threshold and linear. Models with the smallest BIC were selected as the best models for depicting the ISAR of overall OCTE species richness and of each OCTE order separately, respectively.
To analyse the spatial patterns in species richness and composition, we constructed a neighbouring graph (Fig. S1) from a spatial matrix that contained as many rows and columns as the number of GCUs, assigning a value of '1' if two units were neighbours among the 26 GCUs and '0' for 'otherwise' in the ith row and the jth column of the matrix. We linked a unit to its nearest neighbour considering the shortest path between units. We further quantified the spatial autocorrelation among GCUs in terms of species richness using Moran's I (1950) statistics (Cliff and Ord 1973). Positive and significant Moran's I values indicate that nearby GCUs have similar spatial patterns, whereas negative values indicate the opposite. All Moran's I values were further tested using Monte Carlo permutations (999 runs). The above computations were performed for the whole OCTE data and for each OCTE order separately, respectively.
We computed both overall OCTE species richness and species richness of each OCTE group separately using the species presence-absence table elaborated from our literature review. In addition, to assess patterns of endemism, we calculated the number of species that were unique to each GCU, i.e., endemic species. We calculated the endemism rate for combined OCTE and for each OCTE order separately, respectively, as the ratio between the number of endemic species and total number of species.
Prior to statistical analyses, we used a variance inflation factor (VIF) to assess multicollinearity in our variables (Zuur et al. 2010). We separated the variable AGE (Table 1) into three components, i.e. the start of emergence (not used in analysis, see below), the end of emergence and the duration of emergence (difference between the two previous values). We used generalized linear modelling (GLM) to assess the effects of area, distance from the mainland, time of emergence, and local variables and their interactions on freshwater insect species richness. We used a negative binomial error for overall OCT and each OCT order separately, and Poisson error for Ephemeroptera, which gave a more normal distribution of residuals. Models were selected using a stepwise approach based on Akaike's information criterion (AIC). Likelihood ratio tests allowed us to compare the deviance of each selected model and a null model (no effect). To assess the effects of area, distance from mainland, time of emergence, and local variables on overall OCTE endemism rate and on the endemism rate of each OCTE order, we used the logit model with a quasi-binomial error distribution to account for potential overdispersion (see McCullagh and Nelder 1989). For model selection here, we used the quasi-likelihood form of Akaike's information criterion (QAIC). The statistical significance of model parameters was estimated by an analysis of deviance (validated by "F" tests, McCullagh and Nelder 1989) between a null model (no effect) and a model that included the selected explanatory variables. To address the strength of variables in models, we computed effect sizes as the exponential of GLM coefficients, which represents incidence rate ratios (IRR) for negative binomial regression and odds ratios (OR) for quasi-poisson or poisson regression. A value close to 1 represented a small effect size for a given variable. The magnitude of a value above 1 demonstrated a proportional size effect of the variable.
Finally, for assessing the similarity of GCUs in terms of species composition, we used nonmetric multidimensional scaling with the Jaccard distance index (NMDS, Minchin 1987). We further considered two types of grouping: the three groups of GCUs as described in the material and methods (see above) and a grouping provided by cutting the tree generated with Jaccard distance according to the Ward's maximization criterion (Murtagh and Legendre 2014). We further computed an Analysis of Similarities (ANOSIM), which provided an ANOSIM R statistic for each grouping and an associated simulated probability (permutation procedure) to assess the strength of groupings. Geographical variables (area, distance from mainland) and environmental variables (climate heterogeneity, temperature and rainfall), as well as the end and duration of emergence were fitted to the NMDS ordination to allow visualisation of their contribution to the NMDS ordination. The goodness of fit was assessed by a squared correlation coefficient and the significance of the variance explained was obtained from the permutation of environmental variables.

Species richness
Our database of OCTE Caribbean aquatic entomofauna included 1027 species comprising 192 Odonata, 360 Coleoptera, 405 Trichoptera, and 70 Ephemeroptera. Cuba was the richest GCU (390 species, Fig. 2), followed by Hispaniola (290), whereas some islands, such as St Eustatius or Saba, hosted only a few species (6 and 4 species, respectively). Odonata were present in all GCUs, whereas Coleoptera, Trichoptera and Ephemeroptera were absent from 2, 10, and 12 GCUs, respectively (Fig. 3). ISAR analyses demonstrated at least one threshold for overall OCTE species richness and for each OCTE insect orders (Fig. 4, Table 2). Continuous one-threshold models fitted overall OCTE and Coleoptera species richness, respectively, whereas discontinuous two-threshold models fitted Ephemeroptera species richness, respectively (Fig. 4, Table 2). Continuous one-threshold model fitted Odonata species richness if accounting for all GCUs but if omitting the extremely high value of Odonata richness in TOB, the model turned to a linear relationship with island area. Discontinuous two-threshold models fitted Trichoptera species richness but when omitting the extremely low value of Trichoptera richness in BAH, the best fitted model for this order was linear (Fig. 4, Table 2).
Looking at spatial patterns, the overall richness had a significantly positive Moran's I value (MI=0.473, simulated P=0.020). Coleoptera and Odonata had significant positive Moran's I values (MI=0.425,P=0.016,MI=0.463,P=0.008,respectively), whereas Ephemeroptera and Trichoptera demonstrated no  Table 2. Results of the best fit piecewise regression model (according to BIC (Bayesian Information Criterion) for the overall and for each taxonomic group of aquatic invertebrates studied (CO = Continuous One-Threshold, DT=Discontinuous Two-Thresholds, LI=Linear, LT=Left-Horizontal Two-Threshold). A second model is given (in parenthesis) for Overall OCTE, and for Odonata and Trichoptera after omitting Tobago (TOB) and BAH (Bahamas) values respectively. The log likelihood (LL), quality (AIC (Akaike Information Criterion), AICc (Akaike Information Criterion with a correction for small sample sizes), BIC), r-squared (R 2 ) and adjusted r-squared (R 2 adj) are given for each model as well as the log values of area threshold (Th1, Th2). in Table 3) was a main driver of the overall and Odonata and Coleoptera species richness. The end of island emergence (AGE_END in Table 3) had a significant positive effect on the overall richness and separate groups with the exception of Coleoptera and Ephemeroptera Among local variables, climate heterogeneity (CLIM in Table 3) was a main driver of the overall and separate insect group richness, with the exception of Coleoptera. In addition, a minor interaction occurred between DFM and CLIM for the overall richness suggesting that the effect of CLIM was less important if DFM increased and vice versa ( Table 3). Other local variables had a less pronounced effect. Elevation (ELEV in Table 2) influenced the overall richness and separate insect group species richness with the exception of Odonata and Trichoptera. Minimum rainfall (RMin in Table 3) significantly influenced the overall and Odonata. Minimum air temperature (TMin in Table 3) had a very marginal effect on Ephemeroptera species richness.
According to GLM results, geographical variables had no or minor effects (DFM for Odonata and Coleoptera) on the endemism rate ( Table 5). Climate heterogeneity was prominent in terms of explaining the endemism rate overall and that for separate insect groups, with exception of Trichoptera. The other local variables had no or minor effects (ELEV and RMin in Table 5) on the endemism rate.

Species composition
NMDS performed on the species presence/absence provided a fairly good representation of units into two dimensions (Stress= 0.113). Grouping units into five groups provided a better separation (ANOSIM r=0.857; P=0.001) than using the four groups presented in Table 1 (r=0.272; P=0.001). To interpret the grouping, geographical and local variables were fitted to the NMDS ordination ( Fig. 5A and Table 5). Island area (AREA in Table 5), climate heterogeneity (CLIM) and maximum precipitation (RMax) were the main variables associated to NMDS1. They were positively correlated with elevation (ELEV) and negatively Table 3. Effect sizes (IRR (Incident Rate Ratio)) and confidence intervals of variables (in brackets) and their interactions explaining species richness overall and for each order of aquatic invertebrates studied obtained from GLM. Akaike's information criterion (AIC), the likelihood ratio (statistics with *) and the associated probability (P) are given for each GLM model except Ephemeroptera. For this latter group a Poisson model gave a better fit and we considered analysis of deviance (statistics with **). Taxonomic groups are arranged by decreasing dispersal ability (DFM, distance from the mainland; AREA, covered by the GCU; AGE_END, end of emergence process; ELEV, elevation; CLIM, climate heterogeneity; TMin, minimum air temperature; RMin and RMax minimum and maximum rainfall, respectively).

Variable
Overall  Table 4. Number of endemic species and species richness in each GCU. Values are given for overall and for each OCTE group. The number before "/" represents the number of endemic species, and the number after "/" represents the total richness (see acronyms of units in Table 1 Table 5. Effect size (OR (Odd Ratios) and associated confidence intervals (in brackets) of geographical and local variables explaining overall endemism rate and that for each order of aquatic invertebrates studied obtained from GLM. Taxonomic groups are arranged by decreasing dispersal ability (see Table 1 for the acronyms of variables). The quasi-likelihood form of Akaike's information criterion (QAIC), the value of the analysis of deviance (F) and associated probability (P) are given for each GLM model. End or duration of emergence (AGE_END and DURA in Table 5) were not related to the NMDS ordination. The four variables having the highest fit (R 2 values in Table 5) included by decreasing order AREA, CLIM, RMax and DFM. As a result, two gradients could be put to the fore from the NMDS ordination (Fig. 5A) distance from mainland. A second gradient separated the smallest islands (group 2; <80 km 2 ) from mid-size (group 1; <300 km 2 ; group 4; <700 km 2 ) and largest units (>26,000 km 2 ). Group 5 also included large-size islands (~1200 km 2 ).

Discussion
Altogether, our results partially or completely supported our expectations. First, the overall species richness of the Caribbean aquatic entomofauna taxa (OCTE) studied depended on the area of GCUs and, with few exceptions (i.e. Trichoptera if one GCU was removed), demonstrated one or two thresholds in the context of the Small Island effect. Second, ontogenetic stage of islands only partially impacted aquatic insect richness in the expected way. Third, Coleoptera and Odonata species, which have good dispersal abilities, were more widely distributed than Trichoptera and Ephemeroptera, which have poorer dispersal abilities. Fourth, OCTE species richness and composition in the GCUs were weakly influenced by geographic distance from the mainland, whereas local environmental conditions were good predictors of species richness, endemism and composition of OCTE taxa

Species area relationships
Overall, while unit area was a major driver of species richness, our results demonstrated varying responses according to the taxonomic order considered, which  Table 1) along the first and second NMDS dimensions (noted NMDS1 and NMDS2, respectively). Blue vectors correspond to the position of geographical and local variables (see acronyms in Table 1) having a statistically significant r-squared (see Table 6). The five clusters are identified by a number (see text for further details). Mapping of the first (B) and second (C) NMDS unit scores. The size of squares is proportional to the absolute value. Open and closed symbols represent values below and above zero, respectively. Table 6. Relative contributions of the geographical and local variables to NMDS ordination arranged by decreasing squared correlation coefficient (R 2 ), and the simulated probability (P), i.e., statistical significance of R 2 for each variable (ns: P >0.05; see Table 1 for the acronyms of variables). Note that only statistically significant variables (P<0.05) are displayed in Fig. 5 confirm other authors' work showing a taxonomic and context dependency of the ISAR (e.g. Gao et al. 2019). Trichoptera and Ephemeroptera demonstrated no increase in richness with increasing area suggesting the existence of a small island effect (SIE, Lomolino and Weiser 2001) for island below ~760 km 2 surface area. When considering the second "best model", overall OCTE species richness also shows such a SIE. SIEs have been observed several times in the literature and are generally associated with stochastic events that structure small island communities (Gao et al. 2019) contrasting with communities of medium-sized and large islands mainly determined by ecological processes and speciation, respectively (Lomolino and Weiser 2001). In contrast, overall species richness across OCTE taxa (if considering the first "best" model), as well as of Odonata and Coleoptera demonstrated an initial slow rate of species richness increase (below ~79,000 km 2 ) with increasing area followed by a steeper slope concerning essentially Hispaniola (HIS) and (Cuba). These differences in ISAR models among taxonomic groups may reflect ecological differences among the orders. Differences in relief (elevation) regardless of area can specifically explain the SIEs observed for Trichoptera since elevation is known for driving Trichoptera species distribution (Bass 2003). For example, though being the third largest island, BAH being the second lowest GCU in elevation (see Table 1) had only one Trichoptera species. In addition, the Cayman Islands having the lowest elevation do not contain any Trichoptera species. For Ephemeroptera, differences in climate among GCUs regardless of area can specifically explain the observed SIEs, since most Ephemeroptera species present in the Caribbean were found on GCU with high climate heterogeneity (see Table 1).

Ontogenetic changes
The emergence of Caribbean islands occurred at different geological periods and their initial areas depended on the magma production rate and the magnitude of eruptions (Jackson 2013). This geological complexity makes it difficult to detect the effect of ontogenetic changes on species richness as expected in our study. For example, our results showed that the older volcanic islands had a greater specific richness than the younger islands. This is in contradiction with island biogeographical models, which generally predict progressive increase in species richness from younger volcanic islands (early colonisation) toward middle-age islands and a decrease in species richness in older islands due to a progressive decrease in area through time (erosion leading to reduction in habitat diversity) (Valente et al. 2014). Given their ongoing volcanic activities, the small recent volcanic islands are considered by some authors at the beginning of their emergence (Jackson 2013). As a result, their young age could explain their lower specific richness. Alternatively, the positive relationship between age and overall species richness on larger islands may be due to a high speciation rate common in insular tropical environments (Whittaker et al. 2008). This pattern should compensate and even exceed the rate of extinction on the larger islands (Chen andHe 2009, Valente et al. 2014). Given this complexity, only phylogenetic studies could shed light on the effects of ontogenetic changes by estimating colonization and speciation times.

Dispersal abilities of the Caribbean aquatic entomofauna
The values of Moran's I obtained by our analyses were consistent with our predictions on the distribution patterns of our study groups associated to their varying dispersal abilities. Odonata and Coleoptera are considered good dispersers, whereas Ephemeroptera are weak dispersers and Trichoptera are moderate dispersers (Bonada et al. 2012, Sarremejane et al. 2020, and our study showed that Odonata and Coleoptera had higher Moran's I values than Trichoptera or Ephemeroptera. Species richness of all OCTE taxa combined as well as each OC order displayed a significant positive spatial autocorrelation in the Caribbean archipelago, suggesting that close GCUs generally have similar freshwater invertebrate species richness. This resemblance in species richness among close GCUs suggests easier dispersal, which facilitates the movement of individuals from one GCU to another (Bass 2003) since good dispersers are generally good colonizers (e.g., McLachlan 1985). Geographically close units may also share a range of similar ecological niches that favour the establishment of species after arrival (Nekola and White 1999). The fact that the statistical test for Trichoptera gave a positive Moran's I value, whereas Ephemeroptera did not, may be due to their wider range of ecological requirements (Bonada et al. 2004). Indeed, Trichoptera are known to colonize springs, small streams, large rivers, lakes and wetlands (Hering et al. 2009), whereas many Ephemeroptera require more specific habitats (Sartori and Brittain 2015). In our study, among the two taxa with lower dispersal abilities, the taxon with the highest ecological requirements (Ephemeroptera) had fewer species than the richer group (Trichoptera). Thus, we can assume that there is a combination of endogenous (dispersal ability of taxa) and exogenous (environmental conditions, topographical variables) factors that regulate the species composition of the Caribbean OCTE groups, in agreement with the findings of Bonada et al. (2012).

Species richness and endemism
Our results partially confirm the theory of classical biogeography in that island area is a major driver of species richness of OCTE insect communities, whereas the smallest distance to the mainland has only a marginal effect on Trichoptera and Ephemeroptera richness. As a first explanation, the geographical conformation of the Caribbean archipelago adjoining three main lands may cause the OCTE species to have diverse origins. For instance, Crews and Esposito (2020) found that of the 10 genera of terrestrial arthropods Frontiers of Biogeography 2022, 14.3, e54479 © the authors, CC-BY 4.0 license 13 studied, five originated from South America, one from North America, one from Central America, and three from two or all three of these main lands. As another example, in a biogeographical study of two species of Caribbean Heliconiidae (Lepidoptera), D. iulia was shown to originate from South America, whereas Heliconius charithonia originated from Central America (Davies and Bermingham 2002). Only a phylogeographic study of the taxonomic composition of Caribbean aquatic entomofauna and of the adjacent mainland would allow us to clarify this issue, and this was out of the scope of this research. Among local variables, climate heterogeneity was a fairly good predictor for the overall OCTE species richness and for separated OTE orders and the major predictor of the overall and separated endemism of the Caribbean OCE orders. This result that was somewhat expected because several studies have already shown that tropical climatic conditions are one of the major elements responsible for diversification and speciation on islands (Hua and Wiens 2013) as shown for reptiles (Gifford and Larson 2008), amphibians (Rodríguez et al. 2010, Alonso et al. 2011, and freshwater fishes (Doadrio et al. 2009).
The absence of influence of climate heterogeneity on the Coleoptera richness may result from their stronger resistance to environmental constraints due to their physiological and morphological adaptations, such as sclerotization (Hammond 1979, Lawrence andBritton 1994), and they are known to display a large variety of ecologies and life habits (Miller and Bergsten 2016). We found a positive minor relationship between elevation and Trichoptera species richness, which suggests a decrease in richness similar to that observed elsewhere (Miserendino and Brand 2007, Blinn and Ruiter 2009, Maltchik et al. 2009). Indeed, as indicated above, Trichoptera are known to colonize high elevation habitats (Haidekker andHering 2008, Cogo et al. 2020) with higher dissolved oxygen content and lower temperature (Serna et al. 2015) that may favour species diversification (e.g. Bass 2003, Statzner andDolédec 2011). The almost total lack of effects of temperature on the richness of OCTE insects may be explained by the common low temperature range in tropical regions. Only Ephemeroptera species richness was positively but marginally influenced by the minimum air temperature. Ephemeroptera species are sensitive to high temperatures and droughts (Sartori and Brittain 2015). This result may explain their narrower distribution and their absence from small Caribbean islands with more severe climatic conditions. For example, Callibaetis floridanus (Baetidae), known to live in lentic and even temporary habitats (Merritt et al. 2008), was the only Ephemeroptera species present on the Virgin Islands, which have no or few permanent water courses. Barber-James (2007) also recorded the absence of Ephemeroptera on remote islands such as the Tristan da Cunha Archipelago and Gough Island.

Species composition
The species composition of Caribbean islands is the product of multiple factors, including their periodic disturbance by tropical storms, island age and geological origin, which drive the physical and chemical properties of freshwater bodies (Bass 2003). NMDS analysis revealed a spatially organized OCTE species composition among Caribbean GCUs according to geographical and local filtering conditions, but not island age as expected from our fourth hypothesis. The initial grouping based on area and emergence age hardly separated groups of GCUs. In contrast, the species composition of the largest GCUs, with high climatic diversity and elevation, was clearly distinct from that of the smallest GCUs. In addition, the species composition of islands closer to South America differed significantly from those farther (see Fig. 4). We found that the species composition may greatly vary among smaller islands, as already observed by Bass (2003), which underlines the uniqueness of these environments. This difference among groups of geographic units is expected since differences in environmental conditions favour the establishment and diversification of species on larger and more heterogeneous units (Schluter 2000). Oceanic barriers may also limit the movement of individuals according to their dispersal ability (see above), as it has been shown that island isolation limits gene flow among populations and favours differentiation, (Benavides et al. 2007).Whereas geographical variables (area and distance from mainland) are usually seen as major drivers of island fauna (MacArthur and Wilson 1967), we found that local filtering (climate heterogeneity and maximum rainfall) explained as much of the species composition. These results are in the line of the findings of other authors (Lack 1976, Gros-Désormeaux 2012. In addition, the amount of freshwater environment is obviously linked to rainfall (Bass 2003).

Limitation of the study
Our study provides an analysis of the four best documented insect groups (Appendix S2). Only the information for 26 GCUs could be gathered. From the species accumulation curves computed for the OCTE groups, Odonata and Coleoptera seem rather well documented in our database, whereas the three other groups seem far from reaching a "plateau", suggesting the need for including more sites (Fig. S2). However, we accumulated species by GCU, which represented various sampling units. As a result, the plateau may be reach if considering all sampling units. Unfortunately, it is very difficult from checklists, reports or even papers to get the exact locality of samplings, which would have helped to provide more accurate species accumulation curves. We collected information on 6 additional insect orders, using the same documents as for OCTE orders (Diptera, Hemiptera, Hymenoptera, Lepidoptera, Megaloptera, Plecoptera). Only one Hymenoptera family occurred on one unit (DOM). Lepidoptera identified at the family or genus level appeared on 8 units, and one Megaloptera species on one unit (CUB) was found in our literature review. The data that we gathered for Diptera concerned only species known as vectors of disease, which may bias the distribution among islands. Diptera has long been considered as indicating poor water quality whereas this group has a large range of sensitivity to impairment (Serra et al. 2017) and could be given more focus in the area. The two Plecoptera species found in our literature review concerned only Trinidad and Tobago, and should be given more focus as this group is known to be highly sensitive to water quality impairment of (Lenat 1988). In contrast Hemiptera species were well documented and inhabited nearly all GCUs. We decided not to include them in this study to keep a balance between good (Odonata, Coleoptera) and weak dispersers (Trichoptera and Ephemeroptera). In addition, the GCUs were variously documented (see Appendix S2). Cuba was the best documented followed by Guadeloupe, Puerto-Rico, Jamaica, Martinique, Hispaniola and Dominique which cumulate more than 50% of the documents covered. Notably, the island area hardly affected our results since these GCUs represent a range between small (Dominique) and large (Cuba) GCUs. Another limitation of our study is that we did not account for the effect of human populations on species richness patterns in our GCUs. Gao et al. (2019) explained that because larger islands are generally more populated, human-mediated dispersal of species among large islands are facilitated resulting in a slope lowering of SAR. However, this seemed not to affect our results because the slope of SAR for larger islands was steep. In addition, Hispaniola (HIS) comprises two countries (Haiti and Dominican Republic). Haiti has been highly deforested, which could have affected the integrity of rivers whereas the Dominican Republic has better preserved its forested areas. In spite of this, in our data set, Hispaniola remains the second richest island in OCTE.
Building a database like ours was time consuming and we advocate aquatic entomologists working in the Caribbean to implement a database of freshwater invertebrates such as those developed elsewhere (e.g. Schmidt-Kloiber et al. 2019). There is an urgent need to collect this information since freshwater ecosystems are among the most endangered in the world (Dudgeon et al. 2015) and due to its high endemism, the Caribbean rivers require the highest concern.
Conclusions Finally, the combination of geographical and local filtering variables in our study, with the exception of Coleoptera that was mainly influenced by island area and Ephemeroptera that was mainly influenced by local variables, partially supports the biogeographic conceptual model of Lomolino (2000). This result is also in line with a recent study by Vilmy et al. (2021), who showed a generally more balanced situation between dispersal and niche processes for aquatic macroinvertebrates in comparison to other taxa. However, in our study, the local variables were restricted to climate, but other variables, such as the chemical properties of water and the nature of the substrate, require consideration. Our study remains preliminary since our database needs to be further developed by completing our database with other insect groups and conducting phylogeographic research. Such research would allow confirming endemic species and specifying the level and nature of speciation to better understand the origin of Caribbean biota, which is still under debate. Our results confirm that the observation that climatic conditions are one of the major elements responsible for diversification and speciation on islands holds for freshwater invertebrates. The high overall endemism rate found in the Caribbean freshwater biota urges for more conservation effort of this area under the threat of rising anthropogenic disturbances.

Acknowledgements
Chevelie Cineas received funds from the French Embassy in Haiti in the context of a co-tutelle PhD thesis between the University of Lyon 1 (France) and the State University of Haiti (Haiti). This work was also performed under the auspice of the EUR H2O'Lyon (ANR-17-EURE-0018) of Université de Lyon (UdL) within the program 'Investissements d'Avenir' operated by the French National Research Agency (ANR). This article was previously edited for proper English language, grammar, punctuation, spelling and overall style by two highly qualified native English-speaking editors at American Journal Experts (certification verification key: 747A-2FFD-117D-C8B5-6BC3).

Author Contributions
CC and SD conceived the ideas; CC elaborated the database; CC and SD analysed the data; and CC led the writing with assistance from SD.

Data Accessibility
The database supporting the findings presented in this study are accessible via the Dryad digital repositery at https://doi.org/10.5061/dryad.vhhmgqntx.

Supplementary Material
The following materials are available as part of the online article at https://escholarship.org/uc/fb Figure S1: Location of the 26 selected GCUs with a blue line joining them by considering the shortest path Figure S2. Species accumulation curves for each insect order under study. Sites stands for GCUs. [[Q2: Q2]] Appendix S1. List of references and websites used to elaborate the database of aquatic Caribbean insects Appendix S2. Number of documents considered for each freshwater invertebrate group in each GCU Appendix S3. List of documents and websites used to elaborate the environmental data set Appendix S4. List of endemic species of the four insect orders under study and their GCU