From Island Biogeography to Conservation: A Multi-Taxon and Multi-Taxonomic Rank Approach in the Tuscan Archipelago

: Investigating the drivers that support species richness (S) in insular contexts can give insights for the conservation of insular biodiversity. Our aim was to decouple the effect of drivers (island area, distance from mainland and habitat diversity) accounted in three hypotheses or a combination of them in explaining S in seven islands of the Tuscan Archipelago: Area (species–area relationship, SAR), area and distance from mainland (equilibrium hypothesis, EQH) and habitat (habitat diversity hypothesis, HDH). We used published and original datasets to assess S (except aliens) for 42 taxa (14 animal and 28 plant taxa) in each island, and we used S as the dependent variable and the drivers as covariates in regression models. In 31 taxa, the data supported one of the tested hypotheses or a combination of them, and the most commonly supported hypotheses were SAR (12 taxa) and EQH (10 taxa). The effect of the area was also evident in SAR + HDH (ﬁve taxa) and EQH + HDH (one taxon), making it the prevailing driver in explaining S. Since distances are relatively short, and three out of four islands are land-bridge islands, the effect of distance was signiﬁcant for 12 taxa. The effects of habitat diversity were evident for just nine taxa. The multi-taxon approach allowed us to understand the differential effect of drivers among taxa in inﬂuencing S in a single archipelago. Moreover, the multi-taxonomic rank approach highlighted how the information contained within higher taxonomic ranks (e.g., Division) can be substantially different from that derived from lower ranks (e.g., Family). These insights are of particular importance from a conservation perspective of the archipelago’s biodiversity, and this approach can be transferred to mainland fragmented systems.


Introduction
Biodiversity on islands has always had a strong appeal to researchers due to its peculiarities generated through the ecological processes of extinction, colonization and speciation [1][2][3]. Indeed, stochastic or deterministic (i.e., natural selection) events, typical of small isolated islands (bottleneck and genetic drift, inbreeding, demographic and environmental stochasticity), can lead to the formation of unique biological entities [4].
Island biogeography is a consolidated field of research that aims to understand the patterns and processes that determine island biodiversity. Nowadays, clarifying the mechanisms that regulate these diversification processes is fundamental in disentangling the effects of anthropogenic threats on selective factors influencing island populations [5,6].
One of the basic objectives of island biogeography is to understand the relationship that exists between species richness and island size, known as the Species-Area Relationship (SAR [7]). SAR is invariably observed in marine island systems and mainland fragmented landscapes (i.e., "islands" of residual habitats within a disturbed anthropogenic "ocean"), making it the most accepted biogeographical model, defined as "one of community ecology's few genuine laws" [8]. Although the SAR can be described via different functions (see [9]), all agree on the fact that an increase in island area results in the increment of species richness. Large islands can support large populations with low probability of extinction (the Equilibrium Hypothesis, EQH [10][11][12][13]), or they have a higher probability in intercepting dispersal propagules (the Passive Sampling Hypothesis or Target Hypothesis [7,14,15]). Furthermore, large islands generally possess more diverse habitats in comparison to small islands, thus increasing the possibility of hosting a higher number of species (Habitat Diversity Hypothesis, HDH) [16]. The HDH predicts that habitat diversity increases with the increment of island area, and species richness increases accordingly with habitat diversity, while area per se has a minor effect on species richness [16]. Finally, the lower species number on small islands could be the result of stochastic perturbations which disproportionately affect small islands with respect to large islands. Without any doubt, small islands are more vulnerable to disturbance and species richness increases as disturbance frequency decreases (the Disturbance Hypothesis [17]). Despite the general validity of SAR, island species richness can also be affected by other geometric characteristics (sensu [18]), such as island distance from the mainland and the relative position of islands with respect to other islands in an archipelago context (EQH [12,13]). All the mechanisms introduced by the previously illustrated hypotheses to explain the SAR are not in opposition to each other but, rather, can act synergically, emphasizing the positive relation between area and species richness. Disentangling the importance of each mechanism that acts locally as a driver of island species richness can be, on the one hand, crucial to understanding how the pools of insular species are dynamically maintained and, on the other hand, useful for countering the present or potential threatening factors that imperil island biodiversity worldwide.
Although species richness of several groups of organisms can be explained by the biogeography hypotheses, the number of studies that have simultaneously considered a multi-taxon approach are scarce and provide limited interpretations. In particular, what is missing (probably due to the lack of data in some groups) is a multi-level biogeographic analysis that can integrate and interpolate information, deriving from different systematic categories (e.g., Family, Superfamily) and ecological groups of both plants and animals. In this context, the Tuscan archipelago (one of the Mediterranean basin's Biodiversity Hotspots [19]), thanks to the numerous fauna and flora studies, is a natural laboratory in which to overcome the limited taxonomic horizon of previous studies.
In this study, we assessed the importance of the factors argued by EQH, that is, the effect of island area per se, and the combined effect of island area along with its distance from the mainland, on island species richness. In addition, we tested the effects of HDH, exploring if the size of a specific island habitat affected each taxon richness. Additionally, we tested the single and joint contribution of each factor argued by both the EQH and the HDH.
Finally, comparing high taxonomic ranks (Division) of plants, with their nested ranks (Families), we investigated if the drivers, and their magnitude, accounted by the biogeographical hypotheses tested are consistent.

Study Area
We conducted the study on the Tuscan Archipelago, protected under the Tuscan Archipelago National Park and recognized as Man and the Biosphere reserve by UNESCO. The Archipelago is located in the northern part of the Tyrrhenian sea at about halfway between Corsica and the Italian mainland and consists of seven main islands (Elba, Giglio, Capraia, Montecristo, Pianosa, Giannutri and Gorgona) and a series of islets and skerries ( Figure 1). The archipelago has a great lithological variety with significant differences from island to island. All the islands, except Pianosa and Giannutri, have very steep relieves, although they reach modest heights, with the only exception of Monte Capanne (1019 m, in western Elba). Gorgona is made up of calcschists, serpentinites and metabasites; Capraia is of entirely volcanic origin, consisting of andesite flows with tuffs, breccias and basaltic rocks. Pianosa consists almost entirely of a Pliocene biocalcarenite discordant on marls and clays of the Miocene age. Montecristo is a granite mass, Giglio is mostly constituted of monzogranite, Giannutri is carbonatic, while Elba, with its extremely complex geology and structure, consists of a large granodiorite pluto to the West and a well-known sequence of tectonic units with sedimentary successions, as well as other products of quartz-monzonite plutonism to the East [20]. The climate of the Tuscan Archipelago is determined by the position of the islands with respect to Corsica and the Italian mainland, by island size, elevation and by the exposure to winds and sea currents. Among the islands of the archipelago, Elba presents an interesting climatic diversification linked both to its greater size and to the varied geomorphology [21].
This archipelago, rich in endemism, Sardinian-Corsican elements and taxa originating from Italian Peninsula, presents a rather complex biogeographic history linked to its origin during the convergence between Europe and Africa [22], the Pleistocene sea regression and driven by the changes that occurred in more recent times (Miocene to present), including those mediated by humans and their activities [23][24][25][26].
All these characteristics, combined with the information deriving from various studies that have focused on the biodiversity and biogeography of different groups of plants and animals (indicated in Table 1), make the Tuscan archipelago a suitable model for multitaxon analysis, useful for understanding the factors that determine island biogeography in the Mediterranean basin. The sources used to determine species richness on each of the seven islands of the Tuscan Archipelago for all 42 taxa (14 animal and 28 plant) are listed in Table 1. The datasets for Orchidaceae and Apoidea Anthophila, were built ex novo by the authors using published and unpublished original data. For each taxon, only native species were taken into account in the calculation of the species richness. Overall, the time span covered by datasets ranges from 1950 to 2020.
For each island, the geometric characteristics (island area, A, and distance from mainland, D) used to test the SAR and EQH were obtained as a GIS (QGIS desktop 3.14) output, starting from the digital cartography available from the Technical Cartography Digital Portal of Tuscany Region (https://www.regione.toscana.it/-/cartografia-tecnicaregionale-e-scarico-dati-geografici; accessed on 15 December 2020). As distance from the mainland, we assessed both the distance between each island and Tuscan or Corsican coastline ( Table 2). To test the HDH we calculated the Shannon's Index of the third level of land-use types of the CORINE Biotopes (obtained from the Technical Cartography Digital Portal of Tuscany Region and from [43]).

Statistical Analyses
One of the most important assumptions of the equilibrium hypothesis is that the population of each species is positively related to the island area, and it does not depend on habitat diversity. This implies that the probability of extinction of island populations is inversely related to the size of the island. With the aim to identify which biogeographical hypothesis is supported by our data, we first tested SAR using several functions [44,45] in order to detect possible deviations from the best-known Power function (S = c A z [7]). For this purpose, we used the sars package [44] in R 4.0.3 [46]. We tested the two parameters' functions since the number of sampling units is limited (7 islands). We also excluded the linear function, since having a bigger island (Elba) comparing to the others would have produced a biased SAR estimation. The functions used to assess SAR were: Power, Exponential (i.e., Logarithmic [44]), Kobayashi logarithmic, Monod, Negative exponential (the last two show an asymptote). Based on the AICc criterion, in 3 out 42 taxa the Negative exponential (S = d [1−exp (−z A)]) resulted the best function in explaining species richness on the archipelago. For the remaining species, the Logarithmic (S = c + z logA) resulted as the most performing function in explaining the SAR (e.g., [47]), although in many cases it was not very different from the others. In order to easily compare the outputs among taxa we decided to use the same function (i.e., the Logarithmic) in our modelling approach.
Following MacArthur and Wilson [12,13], the equilibrium hypothesis envisions island species richness as a balance between colonization events from the mainland and extinctions of established island populations [12,13]. Thus, in addition to the area effect (accounted in the SAR), species richness on an island is potentially affected by the degree of isolation of the island itself, making the colonization events inversely related to the distance between the island and the mainland. However, isolation may affect the extinction rate as well, making unclear what slope ought to be expected [48]. Colonizing propagules belonging to a species already present on the island can reduce the probability of extinction of their insular population (rescue effect). To test the EQH, we added to the two SAR functions the linear effect of the distance from the mainland. We decided to test the linear effect of distance since the distances between the islands of the archipelago and the mainland do not show extreme values, varying within the same order of magnitude. In this way, greater variability is preserved in comparison to the logarithmic transformation. As distance from mainland, we used the minimum distance between the considered island and Tuscan or Corsican coastline.
The effect of habitat diversity on island species richness is not always easy to test. Often, this is due to collinearity between island area and habitat diversity which compromises the statistical analysis [49]. Therefore, the Habitat Unit Model [50] should be applied to verify the role of the habitat diversity hypothesis in explaining the island species richness. If this mechanism is confirmed, the relative areas of the different habitats should be a better predictor of species richness than the total area of the island. However, in order to apply the Habitat Unit Model, it is necessary to know the association of each species with a habitat among those on the island. A further complication arises from the fact that, while specialist species can be affected by both area and habitat diversity, the more generalist species should be mainly affected by the island area per se. Since we did not know the association between species and habitats, we assumed that island habitat diversity potentially affects the species richness of all the taxa analyzed. Thus, we tested the HDH by accounting for the Shannon's Index calculated on the fractional cover of the different land uses at the third level of the CORINE Biotopes (SI CB ), as a proxy of habitat diversity. The models were performed using the nlsLM (minpack.lm package [51]) and glm (stat package [46]) functions in R for those species whose SAR can be explained by the Negative exponential and the Logarithmic function, respectively. For each taxon, as a dependent variable we used the number of species observed in every island. For the glm models we used a logarithmic link function.
For taxa relying on the Negative exponential function, we built one model for each hypothesis (SAR, EQH and HDH) and the best model was selected according the AICc criterion [52]. Thus, we identified the hypothesis that is supported by the richness of Anthophila, Aquatic beetles and Mammalia.
For the taxa relaying on the Logarithmic function we built a full (multivariate) model, adding to the SAR's component model (i.e., island area), the distance component pertaining to the EQH, as the linear term of the distance from the mainland, and the linear component of the HDH, as the linear term of the SI CB . We excluded a strong collinearity in the multivariate models, since covariate correlations were moderate (Pearson's r: logA-D = −0.34; logA-SI CB = +0.46; D-SI CB = −0.37). To the full model, with the drivers explaining the three hypotheses (area, distance, SI CB ), we applied a multi-model inference approach, in order to identify the drivers that affect the species richness of the analyzed taxa. This allowed us to identify which hypothesis, or a combination of them, is supported by our data. Through an automated model selection, we first identified the best model(s) explaining the species richness observed in each taxon. To this aim we use the function dredge in the MuMIn package [53] in R. We selected, as best models, those with a difference AICc value less than 4 (delta AICc < 4) with respect to the best performing model [52]. According to the drivers selected by the multi-inference approach, we identified the supported hypothesis, or a combination of them, that is supported by the taxa species richness. If the null model was included in the best model set, we considered that no hypothesis was supported by our data. Subsequently, in case of more than one best model, to estimate the mean coefficient of drivers we applied a model averaging, using the function model.avg in the MuMIn package.

Results
The SAR in the Tuscan Archipelago did not show strong deviation from the Logarithmic function, except for Anthophila, Aquatic beetles and Mammalia, which were modelled according with the Negative exponential function. For these three taxa, according to AICc, the best model in explaining species richness on islands supported SAR hypothesis ( Figure 2a). The adjusted R 2 for these models varied from 0.912 and 0.972 (Table 3). The coefficients estimated for SAR (Negative exponential function) in Anthophila, Aquatic beetles and Mammalia were all statistically significant. The coefficient of the asymptote (i.e., d according to the conventional notation) varied from 22.1 (Mammalia) to 142.3 (Anthophila), while the area exponent (i.e., z) varied from 0.022 (Aquatic beetles) to 0.039 (Mammalia) ( Table 4). Table 3. Identification of the supported hypothesis by species richness of taxa modeled by the Negative exponential function according to AICc criterion. AICc: Akaike Information Criterion with a correction for finite sample; R 2 adj: Goodness-of-fit of the model according to the adjusted R 2 .   For the remaining 39 taxa whose richness was modelled using the Logarithmic function, according with the multi-model inference analysis, the species richness supported SAR in 12 taxa (Figure 2b), EQH in 10 (with two taxa responding to distance only; Figure 3), HDH in one (Figure 4), SAR + HDH in seven ( Figure 5), EQH + HDH in one (Figure 6), and no hypothesis in eight taxa (Table 5). In the 31 taxa whose richness supported one or a combination of hypotheses, after the model averaging procedure, the intercept (i.e., c) was statistically significant in 28 and varied from 0.925 (Rosaceae) to 5.872 (Magnoliophyta). In those taxa supporting SAR, the slope of the logarithm of the area (i.e., z) always resulted significant and varied from 0.183 (Caryophyllaceae) to 0.488 (Cerambycoidea), while for taxa supporting EQH, resulted as significant in six out of eight cases, varying between 0.101 (Isopoda) and 0.336 (Rhopalocera). In EQH taxa the slope of distance resulted significant in 4 out of 10 cases and varied from −0.0410 (Chrysididae) to −0.0009 (Fabaceae). The single taxon only affected by habitat diversity (i.e., Apiaceae) showed a slope of the SI CB of 1.492 (p = 0.001). In SAR + HDH taxa the slope of the logarithm of the area resulted as significant in two out of seven cases, varying between 0.052 (Orobanchaceae) and 0.429 (Curculionoidea), while the SI CB resulted as significant in one out of seven cases, with values ranging from 0.080 (Poaceae) to 1.756 (Orobanchaceae). Asteraceae was the only taxon supporting the EQH + HDH, in which the slope of the logarithm of the area amounted to 0.135, that of distance to −0.0043 and that of the SI CB to 0.307 (see Table 6).

Discussion
Our results showed that species richness in 34 out of 42 animal and plant taxa on the Tuscan Archipelago was influenced by island area (SAR), island area and distance from mainland (EQH), habitat diversity (HDH), island area combined with habitat diversity (SAR + HDH) and island area and distance from mainland combined with habitat diversity (EQH + HDH). Among the 34 taxa that respond to one of the tested hypotheses, the most supported by our data was the SAR (15 taxa; 44%), followed by EQH (10 taxa; 29%), SAR + HDH (7 taxa; 21%), HDH (1 taxon; 3%) and EQH + HDH (1 taxon; 3%). If we consider the prevalence of SAR and EQH with respect to HDH, our results differ from those of Sfenthourakis and Panitsa [54], who highlighted the crucial role of the habitat diversity in the islands of the Aegean Sea.
In Anthophila, Aquatic beetles and Mammalia, the Negative exponential function explained the SAR with a great explicative power (R 2 adjusted > 0.912). This is particularly noticeable, since the curve function tends to the asymptote, thus describing the maximum number of potentially available species over a defined threshold and corresponding to the pool of species observable on the archipelago. This curve threshold is reached with different velocity among taxa, with Mammalia having the highest z-value (z = 0.039). We must be cautious in interpreting the Anthophila result, since the asymptote could be an artifact due to an under-sampling of the taxon in the bigger island (Elba) (L.F. personal comment). Conversely, the plateau reached by Mammalia can be considered realistic, given the robustness of the data used (multiple studies are available and probably only few elusive species are missing). Mammals' SAR is attributable to biological and ecological traits of the taxon (e.g., K-strategy, wide home ranges), which affect their population size and, consequently, the survival rate in insular contexts [55,56]. Furthermore, our analyses are in contrast with Hortal et al. [57] and Amori et al. [35], in which habitat diversity resulted as the main driver in determining species richness. Even for Aquatic beetles, the plateau reached by the curve can be considered realistic, since data available for this ecological group are exhaustive. Our result contrasts with what indicated by Ribera et al. [58] and Wohlfart and Vamosi [59], where the main driver in species richness of Aquatic beetles was habitat diversity. This discrepancy is probably attributable to the different ecological needs of the species considered, as the group analyzed in our study is not monophyletic but rather an artificial aggregation.
We found a strong difference in the magnitude of area effect on species richness between Division (Magnoliophyta) or Class (Monocotyledonae and Dicotyledonae) and families (Orchidaceae, Cyperaceae, etc.).
In accordance with Triantis and co-authors [9], we did not find any significant difference in the z between invertebrate and plant taxa (F 1,25 = 0.106, p = 0.748). The intercept (c) reflects the biotic richness of the considered taxon; therefore, its value varies across taxa and strictly depends to the considered taxonomic rank [63]. For these reasons, we did not perform any comparative analysis on this parameter.
The strong insularity observed in Cerambycoidea is probably attributable to their ecological requirements in term of larval hostplants. Since this group includes both generalist and specialist species, developing in herbaceous plants or in dead and decaying wood, their richness may be conditioned by the presence of specialist species, whose habitats more easily occur in large islands [64,65].
In contrast with what we observed in the Tuscan Archipelago, the richness of Tenebrionidae in the Aegean Archipelago seems to be affected by habitat diversity [66]. In our context, the low explanatory power of distance from mainland may be due to the high number of endemic species and sub-species (about 26%). This implies that the Tenebrionidae fauna did to not rely on relatively recent immigration from the mainland, as they are composed of apterous and brachypterous species, but instead is characterized by a large endemic component, speciating after an ancient colonization [E.R. unpublished data]. However, their strong insularity, may be due to the capability of large island in intercepting passive-dispersal propagules from mainland (target effect, [67]) and human-mediated species translocation on larger island [68,69].
Orchidaceae showed a particularly strong insularity in full accordance with what was reported in Lussu et al. [70], where area resulted a better predictor of species richness in comparison to elevation in Mediterranean islands, although Gentili et al. [71] highlighted the highest levels of speciation and diversity in Mediterranean mountains. These plants are commonly present at low population densities in a variety of habitats and this feature makes Orchidaceae particularly susceptible to island size and explains their insularity found worldwide [72].
In contrast with Patiño and co-authors, we observed that the degree of insularity is higher in Pteridophyta than in Magnoliophyta [73]. This discrepancy can be linked to the strict habitat requirements of Pteridophyta. Indeed, under the Mediterranean climate, the available habitats for this group (moist and shady biotopes) are limited and thus more likely found on large islands.
The fact that Chilopoda richness supported SAR may be due to the biology of the group; Chilopoda are in fact predators relying on the abundance of prey species [74] and consequently their richness may be mediated by the insularity of their prey and also for their general low population density [75].
The marked degree of insularity of Cyperaceae may be explained by their high-density occurrence in localized moist habitats [76]. This aspect makes them particularly prone to disappear with the reduction of island size.
Ranunculaceae, in the archipelago, are mostly linked to evergreen oak forest and chestnut woodlands [77]. These habitats are generally uncommon on small Mediterranean islands, making the group particularly sensitive to island area.
In the archipelago, Rosaceae and Lamiaceae are relatively common in mesophilous forests or grassland and ecotone habitats. Their strong insularity may probably be due to the intrinsic rarity or selectivity for localized habitat of some species (e.g., Fragaria spp., Geum urbanum, Rosa rubiginosa, Teucrium spp.) included in the two families [78].
Monocotyledonae, Plantaginaceae and Caryophyllaceae are difficult to justify since they are heterogeneous groups (with different reproductive system, dispersal ability, habitat requirement) in which there does not exist any comparative data in literature. In plants, SAR has been preferentially investigated using the main functional groups like herb (annual and perennial), shrub and wood species [41,79] and consequently comparative data at the family level are uncommon [80].
Magnoliophyta and Dicotyledonae are large groups of species for which we found a typical insularity effect, with a significant contribution of distance (EQH). We used these two wide taxonomic groups as reference for the other taxa whose richness was measured at lower taxonomic rank. Therefore, we are confident that the effect of distance in Magnoliophyta and Dicotyledonae would be averaged, between the most and less distance-sensitive taxa. Thus, we defined three classes of distance effects (d) according to its magnitude found across the taxa: strong (d < −0.01), moderate (−0.01 ≤ d < −0.001, including that of Magnoliophyta and Dicotyledonae) and low (d ≥ −0.001).
Polygonaceae includes species linked to human-disturbed areas as well as to uncommon grassland habitats, other than species characterized by different dispersal capabilities. For these reasons it was not surprising to find them typically affected both from area and distance, although area showed an importance about double of that of distance [81].
For Rhopalocera, accordingly with Dennis et al. [82], Dapporto and Cini [23], besides a prevailing effect of area (i.e., taxon showing a strong insularity), we found a moderate but significant effect of island distance from mainland (14.5% of that of the area) on species richness. Surprisingly, we did not find a significant effect of the habitat diversity on the species richness of the taxon, since it is probably masked by insularity, that is the effect played by the geometric characteristics of the islands.
In Land snails the species richness on the Tuscan Archipelago was affected by area and distance from the mainland (EQH), whose effects was in a ratio of 2.6:1. Land snails, despite being characterized by limited dispersal capacity, contain few endemic species (about 5%) [83]. This scarcity of endemics testifies that this taxon may have relied on more recent phenomena of immigration, also considering that three (Elba, Pianosa and Giannutri) out of seven archipelago islands are land-bridge. For this reason, the effect of distance on land snail species richness is still significant, although not the primary driver (see [84]). Conversely, in the Aegean archipelago, habitat diversity resulted as the most relevant driver, given its higher proportion of land-bridge islands [83].
Fabaceae are a large group including several early-successional specialist species, that benefit from habitat disturbances [85]. This condition may be the reason for which the group was sensitive to island area (large islands are more likely to have habitats in different successional stages or disturbed habitats). Many plant species rely on animal dispersal (zoochory), including birds, and this could be the reason for which the group was moderately affected by distance.
Isopoda supported the EQH, with a moderate area effect and typical isolation degree. This result is in agreement with the findings of Sfenthourakis [86] and Gentile and Argano [87] concerning the significant effect of area. However, these authors found an additional effect of the habitat diversity, while the significant effect of distance we found may be explained by the passive dispersal often adopted by species of the taxon (e.g., floating material transported by currents) [88].
Pompilidae [31] and Chrysididae [89,90] are parasitoid Hymenoptera, both with moderate dispersal capability, whose ecology is strictly dependent on their host presence. These groups are apparently neither influenced by habitat diversity nor area, but strongly by distance. The result confirmed what previously found by Dapporto et al. [31].
Apiaceae is a group characterized by species with very different ecology, colonizing a great variety of habitats (grassland to woodland), coastal to inland areas, early to mature successional stages [91], justifying why the taxon species richness supported the HDH.
In Curculionoidea we observed a joint contribution of drivers supporting SAR and HDH. The strong relationship between the species richness of the taxon and habitat diversity is explainable by their biology closely linked to plants [92,93], as already demonstrated for wood borers [94]. In addition, it should not be overlooked that Curculionoidea assemblages may also have been influenced by plant introduction phenomena that could have acted as vector of insect species originating from surrounding areas (islands and mainland) [32]. Curculionoidea showed also a strong insularity; this is likely due to their stricter specialization respect to Cerambycoidea. Indeed, probably two-thirds of the Curculionoidea are closely associated with plants of one genus [95,96]. In case of host plants with limited distribution, the probability of presence of these plants grows as the island size increases and, consequently, that of the associated Curculionoidea species. Overall, Curculionoidea benefit from both specific limited resources (i.e., rare genera of plants; supporting SAR), but also of diversified resources (i.e., multiple genera of plants; supporting HDH).
Euphorbiaceae and Orobanchaceae are affected by the joined effect of drivers supporting SAR and HDH, with a prevailing of habitat diversity. Despite the limited number of species present in the archipelago for both taxa, they are however extremely diversified, with species ranging from specialist to generalist, making the habitat diversity a strong driver in influencing their species richness [97]. Indeed, several Orobanchaceae are speciesspecific parasites [98]. In addition, although these two taxa showed a moderate degree of insularity, it is not surprising that the contribution of more specialized species on the overall species richness of the families increases as island size increases.
Similarly, the richness of Amaryllidaceae, Asparagaceae and Brassicaceae was affected by area and habitat diversity. However, the effect of habitat diversity was substantially low in the respect to the previous two taxa, while the effect of area was higher, being the area the major driver. Amaryllidaceae and Asparagaceae are bulbous and rhizomatous plants, and are both quite diversified taxa, that benefit from habitat diversity and heterogeneity. Several species of these two taxa are rather stenoecious, needing particular and rare habitats, thus benefit also from large-size islands. In Brassicaceae, that includes a large number of specialized and euryecious species, species richness grows in accordance with habitat diversity and island size.
Poaceae includes both specialist and generalist genera and species capable of growing in multiple types of habitat. The taxon showed a typical insularity, and this is likely due to specialized species, whose uncommon habitats are more likely to be found on large islands. In addition, we found a moderate effect of habitat diversity for Poaceae as already demonstrated by Scherber and co-authors [80].
Reptilia is the only animal taxon that did not supported any of the tested hypotheses. Probably, the geological history of the archipelago exerted a more important role in forging the assemblage of this group, especially during the Pleistocene sea regression. In this period land-bridge islands (Elba, Pianosa and Giannutri) may have been largely colonized (and possibly saturated) by different non-flying animal taxa, including Reptiles, whose assemblages result richer respect to oceanic island (Montecristo, Gorgona, Capraia and Giglio) [24]. Seven plant families (Amaranthaceae, Boraginaceae, Crassulaceae, Geraniaceae, Juncaceae, Papaveraceae and Rubiaceae) did not support any of the hypotheses tested. The reason why we observed such result remains obscure to us, and no other study can provide us some insight. However, we can speculate the fact that some families did not support any hypothesis may be due to their biological or ecological heterogeneity (e.g., different reproductive system, dispersal capability or degree of habitat specialization).

Conclusions
The multi-taxon and multi-rank analysis of plants teaches us that in order to develop an effective biogeographic model that provides useful information to understand the predisposition to extinction of some taxa or their populations, it is preferable to analyze species richness at the family level rather than higher ranks. It is necessary to pay particular attention to the analysis of large taxonomic categories because the drivers that determine their species richness can vary and differ substantially in their contribution with respect to lower taxonomic categories, thus providing potentially misleading output for the conservation of island biodiversity. An example is the Magnoliophyta that respond to EQH, but in which 10 taxa have a strong insularity defined according to SAR, EQH, SAR + HDH, EQH + HDH, four taxa have a species richness determined exclusively by habitat diversity (HDH) and only two in EQH.
These findings may have not only theoretical but also practical relevance. Specifically, lower taxonomic ranks provide more significant information on the overall effect of environmental drivers, since species tend to respond consistently as a function of their phylogenetic relatedness. This interpretation derives from the fact that closely related taxa form homogeneous groups, characterized by shared traits inherited from a common ancestor. In other cases, due to adaptive radiation, the same taxonomic rank may become rather heterogeneous, whereby each species tends to respond in a contrasting way, masking the overall effect of environmental drivers on the group. However, some traits can also be shared between phylogenetically unrelated taxa, due to a partial sharing of the ecological niche resulting from a process of evolutionary convergence. Therefore, to identify and quantify the drivers shaping island biodiversity, the challenge that Island Biogeography has to face in the future is not only to analyze the effect of the drivers on phylogenetically related and homogeneous groups, but also on those unrelated groups sharing common ecological needs.
Indeed, understanding the drivers involved in determining the dynamic equilibria that forge island biodiversity is essential for identifying the most sensible and homogenous (both phylogenetically and ecologically) groups of organisms; recognizing those ecologicalmulti-taxonomic groups would permit the implementation of conservation actions aiming at buffering those natural and anthropogenic factors that determine extinction in susceptible insular systems [99].
From a conservation perspective, to counteract the effects of biogeographical drivers induced by humans, it is mandatory: (1) To preserve small and/or rare but relevant habitats, that are more prone to disappear in small islands, imperiling localized or low-density species linked to these habitats (driver: island area); (2) To maintain insular habitat diversity (driver: habitat diversity); (3) To protect local populations, thus counteracting the "isolation" effect among islands (driver: distance).
This framework can give insights for the conservation not only of island biodiversity, but it can be effectively used to understand ecological dynamics of mainland taxa equally prone to extinction in consequence of environmental fragmentation that leads to habitat insularization.
Concluding, what emerges from our analyzes is also that, in the Tuscan Archipelago, as in many other realities, there is still a huge knowledge gap about the biodiversity of many groups of animals, and invertebrates in particular. This gap is the only real limiting factor for the conservation of biodiversity; above all, it is not safe to infer the sensitivity of an animal taxon to disturbance factors or its predisposition to extinction starting from a wider or different group. It is from this condition that we want to encourage further field research and promote the collection and sharing of data (see [100]), because we cannot adequately protect what we do not know. Funding: This research was partially supported by the Tuscan Archipelago National Park, as part of the activities funded by "BIONETPARKS-La rete delle aree protette per la tutela degli impollinatori naturali".

Data Availability Statement:
Unpublished data used in the analyses are available under request and will be published as separate articles.