Invasion disharmony in the global biogeography of native and non‐native beetle species

The concept of “island disharmony” has been widely applied to describe the systematic over‐ and under‐representation of taxa on islands compared to mainland regions. Here, we explore an extension of that concept to biological invasions. We compare biogeographical patterns in native and non‐native beetle (Coleoptera) assemblages from around the world to test whether beetle invasions represent a random sample of species or whether some families are more prone to invade than others.


| INTRODUC TI ON
"Island disharmony," originally termed by Carlquist (1965), describes the systematic over-or under-representation of certain taxonomic groups on islands compared to mainland source regions. This phenomenon has been repeatedly observed when comparing the flora and fauna of oceanic islands with that of mainland assemblages and is believed to result from selective assembly mechanisms, such as filtering based on dispersal capacity, permitting only a subset of mainland species to successfully colonize islands (Gillespie & Roderick, 2002;König et al., 2020). Traits other than dispersal capability may also act as filters selecting species and include environmental tolerances and dependencies on other species (Taylor et al., 2019;Weigelt et al., 2015). Given the potentially analogous relationship of community assembly via island colonization with assembly via biological invasions (Burns, 2015), it is logical to investigate whether a similar phenomenon of disharmony can be observed in comparing native and non-native species assemblages. Despite their analogous relationship, there are also fundamental differences between biological invasions and island colonizations so we can anticipate that filtering mechanisms driving invasion disharmony may differ from those that drive island disharmony. For example, species traits that promote associations with invasion pathways may increase the chances of species invasions (Kiritani & Yamamura, 2003;Meurisse et al., 2019). Here, we investigate the existence of invasion disharmony by analysing the biogeography of beetle (Coleoptera) invasions worldwide.
Insects are the most diverse group of animals in the world.
Approximately 1 million species have been described and estimates of total extant species (including undescribed species) range from 4 to 10 million (Chapman, 2009;Stork, 2018;Stork et al., 2015;Zhang, 2013). Among insects, the Coleoptera are the largest order.
About 400,000 beetle species have been described and the total number in the world is estimated between 850,000 and 4 million (Bouchard et al., 2017;Nielsen & Mound, 1999;Ślipiński et al., 2011;Zhang, 2013). The evolutionary success of the Coleoptera can be attributed to their diversity of life histories, which include herbivory, predation and detritivory across terrestrial and aquatic habitats (Bouchard et al., 2017).
Given the extent of their diversity, it comes as no surprise that non-native insect species outnumber those from all other animal groups; Seebens et al. (2018) report 10 times as many non-native insect species established worldwide than non-native species in any other animal taxon. Along with the Hemiptera and the Hymenoptera, the Coleoptera are one of the most species-rich orders among the non-native fauna of various world regions (Liebhold et al., 2016;Wheeler & Hoebeke, 2009). The dominance of Coleoptera among non-native insects appears to reflect their dominance in the entire global fauna; the proportional representation of Coleoptera species in non-native insect assemblages is generally similar to their representation in native species assemblages (Liebhold et al., 2016;Yamanaka et al., 2015).
While most non-native Coleoptera are relatively innocuous in their invaded ranges, a few species cause extensive damage and are considered serious pests. For example, the rice weevil, Sitophilus oryzae, is native to India but has invaded most world regions, sometimes damaging stored rice such that it has no value as food (Longstaff, 1981). The Colorado potato beetle, Leptinotarsa decemlineata, is a major invasive agricultural pest; the species had a limited native range in western North America but has now invaded virtually all regions of the Northern Hemisphere where potatoes are grown, sometimes causing massive crop damage which has occasionally resulted in famines (Clark, 2007). Some non-native beetle species have serious impacts on biodiversity and other non-market values.
For example, the emerald ash borer, Agrilus planipennis, is spreading across North America and Eastern Europe, causing widespread mortality of ash trees, Fraxinus spp. This damage is anticipated to cost billions of dollars to homeowners and municipalities and extensive loss of biodiversity in organisms associated with ash trees (Herms & McCullough, 2014).
A key question to understanding and managing invasions is the quantification of the risk that species may invade new areas in the future. Such risk analysis is critical to the implementation of biosecurity practices aimed at excluding invasions of new species (Robinson et al., 2017). Even though many species of Coleoptera are known invaders and have caused extensive damage, information is lacking on invasion probabilities for different groups within the Coleoptera.
Thus, the discovery of invasion disharmony may provide crucial information that would lead to more effective risk analysis for setting biosecurity policies targeting these insects. Here, we analyse comprehensive lists of non-native Coleoptera species from several regions in the world, comparing numbers of non-native species among families and the proportional representation of families in native and non-native assemblages. Our objective here is to use this information to quantify invasion disharmony and assess the relative invasion risk among various beetle families. Specifically, we test

K E Y W O R D S
Coleoptera, composition, disharmony, family, insect, invasion, native, non-native whether beetle invasions represent a random sample of the order or whether some families are more prone to invade. We hypothesize that much like island disharmony, an analogous filtering of species occurs during biological invasions that produces systematic compositional differences between native and non-native assemblages. We also test whether the family-level composition of species invading a given region is similar to the composition of native species in that region or whether they more closely resemble the composition of other invading assemblages.

| Compilation of data on numbers of native and non-native species
We assembled lists of non-native Coleoptera species from 10 regions worldwide, namely North America (excluding Mexico), Japan, the Okinawa and Ogasawara Islands, the Hawaiian Islands, South Korea, Europe (including the European part of Russia), New Zealand, Australia and the Galapagos Islands. These lists can be downloaded from Dryad at https://doi.org/10.5061/dryad.tx95x 69z1. Species lists originated from sources listed in Appendix S1 and were updated from various miscellaneous sources (e.g., recent publications and websites). A full list of species can be downloaded from Dryad (https://doi.org/10.5061/dryad.tx95x 69z1). These 10 regions were selected because they represent the few regions worldwide with existing comprehensive lists of non-native beetles. The regions vary considerably in size but cover extensive geographical diversity (e.g., the Northern and Southern Hemispheres, tropical and temperate climates, continental and island regions). We acknowledge that these regions mostly coincide with countries with highly developed economies, which may introduce unavoidable bias. Furthermore, we recognize that these lists of non-native species may actually be slightly incomplete; there typically are lags between establishment and discovery and reporting of new non-native species (Essl et al., 2011;Morimoto et al., 2019). Species that were known to have been intentionally introduced (e.g., biological control agents) were not included in our analyses.
Lists from each region were standardized to overcome duplication due to the existence of synonyms and misspellings, by performing taxonomic "cleaning" so that all species names and higher-level taxon designations were based on a single taxonomic classification system. This was performed using the GBIF taxonomic database (GBIF Secretariat, 2019) and the "taxize" package in r (Chamberlain & Szöcs, 2013). Code used for this taxonomic cleaning can be downloaded from the Zenodo repository (https://doi.org/10.5281/zenodo.4555787). A total of 91% of Coleoptera species user-supplied names were recognized (including as synonyms or misspellings) by GBIF. For the remainder (177 species), standardization was performed manually via searches of alternative databases and manual researching of names.
Like in many other taxonomic groups, there are frequent changes in the higher-level taxonomy (family groupings) within the Coleoptera. For consistency, we adopted the classification of Bouchard et al. (2011) with three exceptions. Dryophthoridae and Brachyceridae were merged into the Curculionidae, Passandridae were merged into the Cucujidae, and Megalopodidae were merged into the Chrysomelidae to align with source datasets.
The composition of species assemblages was characterized at the family level. Native and non-native assemblages from each region were summarized by the number of species per family. To accomplish this, we used our non-native species lists to calculate numbers of species per family for each region. We also compared the composition of non-native assemblages to the composition of native assemblages from the same region. Numbers of species per family for native assemblages were compiled from sources given in Appendix S1. We also compared the composition of non-native assemblages to that of all known world Coleoptera species, hereafter referred to as "world described species." The composition of the world Coleoptera was based on a published list of estimated numbers of all described Coleoptera species in the world by family (Bouchard et al., 2017).
Based on the lists of non-native species in each of the 10 regions, we calculated the total number of non-native species in each family pooled among all regions combined. To limit the stochastic effects of smaller families, we restricted our analyses to include the 25 families with at least 10 species known to have invaded at least one of the 10 regions.

| Species-area relationships
To analyse the influence of varying region sizes on species numbers, species-area relationships were explored by plotting numbers of species per regional assemblage vs. regional land area. Plots were made using log scales for both axes, and the slope and R 2 were calculated for the linear regression relating log species numbers to log land area. Species-area plots were generated using the total number of native and non-native Coleoptera species from each region.
Similarly, species-area plots were generated using numbers of native and non-native species in the Curculionidae and the Staphylinidae (the two most species-rich Coleoptera families) from each region.

| Proportional representation by families
For each of the 10 regions, we compared numbers of species per family in the non-native assemblage vs. the number of native species of the same family in that region and vs. the total number of known species per family in the world. We calculated the expected number of species per family assuming an equivalent proportion of species in each family as in the reference assemblage (i.e., the native species assemblage). To illustrate the variation in each region, we created scatterplots of numbers of non-native species per family vs. numbers of species in the same family for the reference assemblage.
We plotted a line showing expected numbers under equivalent representation and a region around this line indicating where the probability that a family was over-or under-represented (because they fell outside of the boundaries) was less than or equal to α = 0.01.
We used a Bonferroni correction for the multiple comparisons, so the boundaries of this region were calculated as the upper and lower quantiles of the binomial distribution such that (1α /m) × 100% of the distribution lay within the boundary, where m is the number of families compared.
We also investigated whether there was a tendency of families of certain feeding groups to be over-or under-represented. To accomplish this, we first classified each family into one of three feeding groups (herbivores, predators and scavengers) based on their dominant life history described in Arnett and Thomas (2001) and Arnett et al. (2002). Each family was then categorized according to whether they had more than the expected number of non-native species (in the 10 regions pooled) with reference to the world described species assemblage. Pearson's chi-squared test was used to test whether the proportion of families with more than the expected number of nonnatives was significantly dependent on their feeding group.

| Similarity between native and non-native assemblages
To quantify similarity among native and non-native assemblages, we computed a full correlation matrix between all native and nonnative assemblages using numbers of species per family. As is characteristic of community composition data, distributions of numbers of species per family were highly skewed. Thus, for the correlation and ordination analyses (described below), species richness of each family in each of the assemblages was transformed using a Hellinger transformation (square root of the number of species per family divided by the total number of species in all families) (Legendre & Gallagher, 2001).

| Ordination analysis
We employed a direct ordination, redundancy analysis (RDA), to characterize differences among all assemblages (both native and nonnative) based upon the distribution of species among families in each assemblage (ter Braak, 1986;Legendre & Legendre, 1998). Analyses were again based on Hellinger-transformed values. Following ordination, each assemblage was plotted using their scores for the first two RDA axes; the position of each assemblage in this space provided a map of compositional similarities (or dissimilarities) among assemblages. We also plotted the RDA scores of each family in order to visualize how the differences among assemblages were related to the relative dominance of different families. The effect of status (native or non-native) and of regions (including both the pooled set of all non-native species and the world described species) was evaluated by a permutation test based on 999 permutations. Neither crossed nor nested effects were considered. The RDA and permutation tests were computed using the "vegan" package in r (Oksanen et al., 2020).

| Numbers of species per region
There were 1,967 non-native beetle species recorded among all regions. North America had the most recorded non-native beetle species followed by the Hawaiian Islands, Europe and New Zealand (Table 1). The region with the fewest recorded non-native beetle species was the Galapagos Islands followed by the Ogasawara Islands and South Korea. Australia, North America and Europe each reported slightly more than 25,000 native Coleoptera species while the Galapagos and Ogasawara Islands both reported the fewest native species, with just slightly more than 300. Given this information, the Hawaiian Islands had the largest number of non-native species relative to the number of native species though the ratio was nearly as high for the Galapagos Islands. Australia, Japan and Europe had the fewest non-native species relative to the number of native species.

Region
Non-Native species Native species Ratio

| Species-area relationships
Much of the variation among regions in total numbers of both native and non-native Coleoptera species was explained by land area. Total numbers of non-native species per region followed a classic linear species-area relationship in a log-log plane (Figure 1b). A similar pattern was observed for numbers of native species (Figure 1a), as well as for both native and non-native Staphylinidae and Curculionidae

| Proportional representation by families
The numbers of species per family for the world described species, for families with at least one non-native species. Families are coloured by their principle feeding group. The black line shows the expected number of non-native species per family if they were in the same proportions as in the world described species. Grey shading shows the range outside of which families would be deemed to be over-or under-represented at the α = 0.01 level, assuming a binomial distribution and using a Bonferroni correction to account for the number of families compared. The labelled families are the most extreme-those more or less than the boundaries of the grey region Dermestidae, Bostrichidae, Ptiliidae and Latridiidae, were overrepresented ( Figure 3). There was no consistent trend for families that were predominantly herbivores, predators or scavengers being over-or under-represented in the non-native assemblage. While Figure 3 suggests a weak trend of over-representation of predominantly scavenger families, the distribution of feeding groups among those families with more non-native species than the expected was not significantly different to the distribution among families with fewer species than expected (Pearson's chi-squared test, p =.87, for families with at least one established species, where the expected number of non-native species per family is proportional to the number in the world described species).
Relationships between the composition of native and nonnative assemblages in each region (Figure 4 and Figure S2.2) were similar to the relationship described above for all regions (Figure 3). The

| Similarity between native and non-native assemblages
Correlations in numbers of species per family between assemblages describe similarities in the composition of assemblages, both native and non-native ( Figure 5). Overall, there was more similarity among native assemblages than among non-native assemblages and F I G U R E 4 Heatmap of number of species within each family for native and non-native Coleoptera assemblages from each region. Numbers are expressed as square root (number of species in family / total number of Coleoptera species in assemblage) the native assemblages were most similar to the world described species. Among native assemblages, those of North America and Europe were the most similar ( Figure 5). Native assemblages from some of the East Asian regions, such as Japan, Korea and Okinawa, were also quite similar to each other. Among non-native assemblages, Australia, New Zealand, North America and the Galapagos were generally similar to each other and the East Asian regions of Japan, Ogasawara, Okinawa and South Korea were also similar to each other.

| Ordination analysis
Ordination analysis (RDA) indicated that the family-level composition of non-native assemblages was distinct from that of native assemblages since the native and non-native assemblages generally fell on opposite sides of the ordination space, primarily defined by the first axis ( Figure 6a). The significance of native vs. non-native status was confirmed by the permutation test (Table 2). However, there was no significant effect of region which indicates that native and non-native assemblages from the same region are not more similar to each other than to other assemblages. Even though native assemblages in the East Asian regions Japan, S. Korea, Okinawa and Ogasawara Islands all had similar composition, non-native assemblages were not generally similar among these regions (Figure 6a). Likewise, the family-level composition was similar among North America, Europe, Galapagos and Hawaii native assemblages, but there was generally little similarity among the non-native assemblages from these same regions

| D ISCUSS I ON
The Coleoptera are the most species-rich taxonomic group in the world (Bouchard et al., 2017). Here, our analysis adopts the number of 385,602 described Coleoptera species from Bouchard  Grove and Stork (2000) estimate that 70% to 95% of all beetle species remain undescribed. We report a total of 1,967 species that are non-native in at least one region, but, like estimates of total numbers of beetles in the world fauna, the true total number of non-native beetle species may be substantially higher. First, the 10 regions from which comprehensive lists of non-native beetle species were available represent only a fraction of the world. In addition, there typically exists a lag between the establishment of non-native species and their discovery (Essl et al., 2011;Kiritani & Yamamura, 2003) so it is likely that undiscovered non-native species exist in these regions. Figure 1 suggests that much of the variation in numbers of species, both for non-native and for native species, among regions can be explained by differences in land area. The existence of a log-log linear species-area relationship is perhaps the most ubiquitous biogeographic phenomenon known to science, having been observed in virtually every taxon and world region (Lomolino, 2000) including both native and non-native species assemblages (Blackburn et al., 2016). Sax and Gaines (2006) reviewed studies that compared species-area relationships for both native and non-native assemblages and noted that relationships for non-native species are often different from those of native assemblages. Our data suggest that species-area relationships for non-native assemblages were generally not as steep as those for native assemblages from the same regions though none of them were significantly different ( Figure 1).
Further, lower R 2 values observed here for non-native assemblages indicate that land area explains less variation in numbers of these species than it does for native assemblages. Numbers of non-native species in a region may be influenced by a multitude of other factors that are related to propagule pressure or habitat invasibility (Lonsdale, 1999;Simberloff, 2009 (Buckland, 1981). The Curculionidae also contain the true bark beetles and ambrosia beetles (Scolytinae) which are notorious hitchhikers in wood and wood packaging that have extensively moved around the world in trade (Brockerhoff et al., 2006). these insects historically may have been moved globally with soil, rock and wood bark (Ødegaard & Tømmeras, 2000). Lindroth (1957) chronicles in detail the historical accidental transport of Carabidae from Europe to North America with soil and rock used as ballast in sailing vessels prior to 1900. Like the Carabidae, Staphylinidae are likely to have been transported in this manner. During this era of large clipper ships, trade routes were intense between Europe and North America and to a lesser extent between Britain and its colonies (e.g., Australia, New Zealand) but connectivity to east Asia was less strong (Pascali, 2017). While it remains uncertain, it is possible that these historical trade routes explain the much lower numbers of Staphylinidae in east Asian regions analysed here. However, an alternative explanation for the highly variable proportion of Staphylinidae is the existence of differences in the intensity of reporting or discovery for this group among regions. This family generally does not include pest species and therefore these species may be overlooked.
It should also be noted that numbers of native Staphylinidae species were also proportionally low in Ogasawara, Okinawa and South Korea (Figure 4).
Beetle species that were generally over-represented in nonnative assemblages included the Dermestidae, Bostrichidae, Ptiliidae and Latridiidae. The majority of these over-represented families belong to Bostrichoidea and Cucujoidea superfamilies.
The Bostrichidae mostly feed on wood so, like the Scolytinae mentioned above, have likely been transported extensively with raw wood, wood products and wood packaging material (Brockerhoff The family-level composition of non-native assemblages is completely distinct from the composition of native assemblages; all non-native assemblages were separate from all native assemblages in RDA ordination space (Figure 6a). This result resembles those of Liebhold et al. (2016) who analysed species richness by order across 44 world regions and found that the order-level composition of native assemblages was distinct from those of non-native assemblages.
Variation in family-level or species-level native species richness among regions reflects long-term (evolutionary) biogeographic radiation that drives global beta-diversity (Novotny & Weiblen, 2005).
Biogeographic factors that determine this intercontinental variation in composition include latitudinal gradients and continental drift.
However, the composition of non-native assemblages is driven by very different factors, namely the propensity of species to become associated with invasion pathways and the establishment success of small nascent populations, both of which are influenced by lifehistory traits that vary among insect orders and families (Kiritani & Yamamura, 2003;Liebhold et al., 2016;Meurisse et al., 2019).
Indeed, RDA redundancy analysis indicated that there was no significant effect of region which indicates that native and non-native assemblages from the same regions are not significantly more similar to each other than to other assemblages, and likely reflects the very different processes that determine their composition (Table 2).
For example, while the family-level species richness of native assemblages was similar among east Asian countries (all had positive values in the secondary RDA axis), the corresponding non-native assemblages were much less similar among regions (Figure 5a). Variation in the composition of non-native assemblages among regions may be driven, in part, by differences in the origins and types of goods that they import. Consequently, the composition of a non-native assemblage in any region tends to be more similar to non-native assemblages from other regions than it is to the composition of the native assemblage in the same region.
The differences that we report here between the family-level composition of native beetle assemblages and the composition of non-native assemblages are in many ways analogous to the concept of island disharmony that has been used to describe systematic overand under-representation of certain taxa on islands (Carlquist, 1965).
While our ordination analysis indicated a coherent distinction in the family-level composition of native and non-native assemblages, it showed no evidence of a consistent difference in the composition of native beetle assemblages between mainland and island regions ( Figure 6a). Thus, results here suggest that invasion disharmony is generally greater than island disharmony in the Coleoptera.
Similar to invasion disharmony observed here in the Coleoptera, there are several reports of differences in the family-level composition between native and non-native plant assemblages and these differences are attributed to selection for plant species exhibiting traits that promote invasion (Procheş et al., 2008;Pyšek, 1998).
Island disharmony is also believed to result from filtering based on species characteristics. These characteristics include dispersal capacity and disharmony may intensify as a result of radiation from early colonizers (Finston & Peck, 2004;Givnish et al., 2009). Recent work indicates that island disharmony may also arise from filtering of species based on traits other than dispersal capacity (Taylor et al., 2019;Weigelt et al., 2015). For example, colonization of islands may be more successful for species that do not depend on specialized interspecific interactions given that key species may be missing on islands. Taylor  are predominantly herbivorous, many other families (e.g., Staphylinidae, Carabidae, Tenebrionidae) are comprised of mostly non-herbivorous species. We found no evidence for a consistent effect of dominant feeding category on invasiveness (Figure 3).
While the impacts of many herbivorous non-native Coleoptera species are well known because of their impacts on agriculture and forestry, the impacts of non-herbivorous species are much more poorly understood. Though some predaceous Coleoptera species have been introduced for purposes of biological control (these were excluded from our analyses) and some information exists about their impacts, most predatory beetle species were not intentionally introduced, and almost nothing is known about their ecological impacts. Because non-native predators generally have greater impacts than native species (Paolucci et al., 2013), it is quite likely that some of these species have profound impacts on ecosystem structure and function. Similarly, Coleoptera scavenger species play key roles in nutrient cycling and other crucial ecosystem processes (Yang & Gratton, 2014) yet almost nothing is known about how the many non-native coleopteran detritivores and fungivores are altering these processes. In short, our current knowledge of both the macroecology and ecosystem ecology of beetle invasions is limited and requires extensive exploration.

| CON CLUS IONS
Here, we propose the term "invasion disharmony" to refer to systematic over-and under-representation of higher-level taxa among non-native assemblages compared to native assemblages. This phenomenon is analogous to "island disharmony" which is widely seen when comparing the composition of island and mainland faunas and floras (Gillespie & Roderick, 2002;König et al., 2020). In a manner similar to the colonization of islands, transport and establishment of non-native species filter the pool of potentially invading species and drive differences based on life history and other characteristics of species. However, invasion disharmony differs from island disharmony in that uniqueness in the higher taxa make-up of nonnative assemblages is not driven by radiation from early colonizers which may drive random differences between island and mainland communities.
We have shown that the fraction of species represented by each family varies considerably among native and non-native beetle as-

CO N FLI C T O F I NTE R E S T
None of the authors have any conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data available for download from Dryad include a table of numbers of non-native Coleoptera species per family in each of 10 world regions, a table of numbers of native Coleoptera species per family in each of 10 world regions and a list of non-native Coleoptera species from 10 world regions. These data can be downloaded at https://doi. org/10.5061/dryad.tx95x 69z1.