How Phylogenetics Can Elucidate the Chemical Ecology of Poison Frogs and Their Arthropod Prey

The sequestration by neotropical poison frogs (Dendrobatidae) of an amazing array of defensive alkaloids from oribatid soil mites has motivated an exciting research theme in chemical ecology, but the details of mite-to-frog transfer remain hidden. To address this, McGugan et al. (2016, Journal of Chemical Ecology 42:537–551) used the little devil poison frog (Oophaga sylvatica) and attempted to simultaneously characterize the prey mite alkaloids, the predator skin alkaloids, and identify the mites using DNA sequences. Heethoff et al. (2016, Journal of Chemical Ecology 42:841–844) argued that none of the mite families to which McGugan et al. allocated the prey was thought to possess alkaloids. Heethoff et al. concluded from analyses including additional sequences that the mite species were unlikely to be close relatives of the defended mites. We re-examine this by applying more appropriate phylogenetic methods to broader and denser taxonomic samples of mite sequences using the same gene (CO1). We found, over trees based on CO1 datasets, only weak support (except in one case) for branches critical to connecting the evolution of alkaloid sequestration with the phylogeny of mites. In contrast, a well-supported analysis of the 18S ribosomal gene suggests at least two independent evolutionary origins of oribatid alkaloids. We point out impediments in the promising research agenda, namely a paucity of genetic, chemical, and taxonomic information, and suggest how phylogenetics can elucidate at a broader level the evolution of chemical defense in prey arthropods, sequestration by predators, and the impact of alkaloids on higher-order trophic interactions.


Introduction
Systematics and taxonomy have aided our understanding of the evolution of chemical defense and its interplay with trophic ecology in diverse groups of plants and animals (Berenbaum 1995). Phylogenetic studies of the evolution of chemical defense in predators and prey, built on a foundation of well-documented biodiversity, have begun to demonstrate the extent to which ecological interactions among species are mediated by toxic chemicals. For example, the phylogenetic distribution of cardenolide toxins across distantly related angiosperms indicates that the biosynthesis of this defense has evolved in parallel and been lost repeatedly over evolutionary time (Agrawal et al. 2012). Indeed, progesterone 5β-reductase orthologs, which code for a conserved enzyme needed for cardenolide biosynthesis, are present in both cardenolide-free and cardenolide-producing species. This suggests that cardenolide biosynthesis might represent an ancient pathway with multiple evolutionary losses across groups (Agrawal et al. 2012). Both the host plants and herbivores resist cardenolides by target site insensitivity of the Na + /K + -ATPase gene due to selection for two parallel resistance-conferring amino acid substitutions in five insect orders (Agrawal et al. 2012;Ujvari et al. 2015;Yang et al. 2019). Diverse songbird species at a third trophic level have evolved Na + /K + -ATPase mutations that confer them with resistance to their co-distributed cardenolide-containing insect prey (Groen and Whiteman 2021). Similar molecular changes in the Na + /K + -ATPase gene also explain resistance in reptiles, amphibians, and mammals (Ujvari et al. 2015). The shared molecular underpinnings of resistance indicate that evolutionary responses to chemicals may converge across the animal kingdom.
Taricha salamanders are defended from predators by highly toxic tetrodotoxin (TTX) produced by symbiotic skin bacteria (Vaelli et al. 2020). Phylogenetic analysis shows that three species of garter snakes (Thamnophis) have independently evolved different amino acid substitutions in the sodium channel Na v 1.4 that provide the snakes a degree of resistance to TTX (Feldman et al. 2009).
The pitohuis of New Guinea (Aves) are defended by homobatrachotoxin alkaloids sequestered from melyrid beetles (Choresine) which repel chewing lice on the bird feathers (Dumbacher 1999;Dumbacher et al. 1992). Moreover, phylogenetic analysis showed that acquisition of defense in Pitohui and close relatives has occurred multiple times (Dumbacher and Fleischer 2001;Dumbacher et al. 2008) and that the genus Pitohui was polyphyletic (Dumbacher 2014). Interestingly, their batrachotoxin alkaloids are shared with Phyllobates poison frogs (Tokuyama et al. 1969), which have been hypothesized to sequester them from the same genus of beetles (Dumbacher et al. 2004).
Species of several frog taxa (e.g., Mantellidae, Eleutherodactylidae) sequester dietary chemical defenses. Perhaps the best-studied are the Neotropical poison frogs (Dendrobatidae), in which 102 of the 333 species (AmphibiaWeb 2021) are known or likely to be aposematic (defended and visually conspicuous) and typically diet-specialized (Darst et al. 2005;Santos et al. 2003;Vences et al. 1998). This "aposematic syndrome" has evolved multiple times (Santos and Cannatella 2011). Soil mites of the suborder Oribatida have been recognized as one of the most important sources of defensive alkaloids in poison frogs (Saporito et al. 2007), which sequester the alkaloids and store them in granular skin glands. These alkaloids are presumed to be secreted by oil (opisthonotal) glands; however, direct evidence of endogenous alkaloid production by these glands is lacking (Raspotnig et al. 2011a;Saporito et al. 2015). The frogs release the alkaloids when threatened by predators such as snakes, birds, and spiders (Alvardo et al. 2013;Darst et al. 2005;Neuwirth et al. 1979;Santos and Cannatella 2011). Arthropod predators and parasites that rely on chemoreception, such as neotropical bullet ants (Paraponera clavata), red imported fire ants (Solenopsis invicta), red-legged banana spiders (Cupiennius coccineus), and yellow-fever mosquitoes (Aedes aegypti) are effectively deterred by poison frog alkaloids as well (Fritz et al. 1981;Murray et al. 2016;Stynoski et al. 2014;Weldon et al. 2006Weldon et al. , 2013. Alkaloids are likely to be keystone molecules (Ferrer and Zimmer 2013;Zimmer and Ferrer 2007), impacting ecosystems more than would be anticipated based on their concentrations (about 1,000-10,000 times lower than those of ubiquitous compounds such as free amino acids, monosaccharides, and nucleotides) (Ferrer and Zimmer 2013). A well-studied example of a keystone molecule is TTX, which is involved in trophic interactions within riparian and montane communities (Williams et al. 2004;Zimmer et al. 2006). In addition to its defensive role in adult Taricha newts, TTX elicits escape reactions in conspecific larval prey at smaller doses (Zimmer et al. 2006). TTX also protects Thamnophis garter snakes from its predators, namely avian raptors (Williams et al. 2004;Zimmer et al. 2006). Raptors, snakes, and newts, in turn, have large impacts on community assemblages (Davic and Welsh 2004;Jones et al. 2001;Kurzava and Morin 1998;Marti et al. 1993;Smith 2006). Indeed, experimental work has suggested that the presence of newts leads to differential mortality and reverse competitive hierarchies between native and non-native tadpole species (Smith 2006). In systems with both newts and fish, unpalatability of anuran prey decreased predation by fish but was ineffective against newts (Kurzava and Morin 1998). Thus, keystone molecules can fine-tune community assemblage and prey community composition.
We illustrate a general scenario in Fig. 1. Frogs eat arthropods, some of which (ants and mites) are defended by different classes of alkaloids, such as pumiliotoxins and epibatidine. Most frogs are not aposematic (center row, left); they lack chemical defenses and warning signals. They are also typically diet-generalists and ingest a variety of prey, including certain defended species (perhaps in smaller quantities), but may also reject other defended species. They do not sequester defensive chemicals, and are not defended from predators such as birds and snakes (top row). In contrast, aposematic poison frogs (center row, right) are typically defended, conspicuous diet-specialists that have evolved resistance to prey chemical defenses, even at times preying on multiple species with different chemical defenses (ants and mites) but also consuming undefended prey. These sequestered defenses may protect the aposematic frogs from certain predators (birds), but not others (snakes) that have themselves evolved resistance to alkaloids. Thus, poison frogs, their prey, and their predators offer an intriguing system of overlain webs of species interactions. Thus, alkaloids may be keystone molecules determining the trophic structure and composition of communities.
Dissecting this complexity requires fine-grained information on predators, prey, and chemical defense. Several aposematic frog species, and a few non-aposematic ones, have been tested for skin alkaloids, but in almost all cases except Oophaga pumilio (Prates et al. 2019;Saporito et al. 2007), the specific arthropod source of the alkaloids remains a mystery. McGugan et al. (2016a;hereafter McGugan et al.) used a novel phylogenetic approach to elucidate this mystery. They attempted to identify the arthropod species (ants and mites) responsible for alkaloid diversity in the little devil poison frog, Oophaga sylvatica, by coupling molecular phylogenetic identification of mite prey with the identification of alkaloids from the mite prey and the skin of the frog predators. McGugan et al. (2016a, b) analyzed both ants and mites from the stomachs of frogs of three populations of O. sylvatica. Here we focus only on the mites. To pinpoint the source of the alkaloids, they compared matched samples of frog skin and mite prey from five frogs using LC/MS (liquid chromatography/mass spectrometry) to identify structural classes in common and GC/MS (gas chromatography/mass spectrometry) to identify the alkaloids. To identify the mites, they sequenced cytochrome oxidase 1 (CO1), a fast-evolving gene used for "barcoding" species identification (Hebert et al. 2003), from nine mites and performed BLASTn searches of these sequences against the NCBI Nucleotide database. They matched one sequence to a well-known species (Archegozetes longisetosus) based on a percent similarity of 99%. However, the other eight mite sequences had similarities of only 78-84% to GenBank records, with resolution limited to family or genus; they tentatively identified the mites as belonging to the Ceratozetidae, Eremaeidae, Oppiidae, and Trhypochthoniidae.
McGugan et al. faced swift criticism of their analysis from a group of mite experts (Heethoff et al. 2016, hereafter Heethoff et al.) Fig. 1 How keystone chemicals may govern trophic relations in a community. The red lines indicate the transfer of chemical defenses (along with energy resources) between trophic levels, and black lines indicate lack of defenses. The thickness of the lines indicates the amount of transfer (large or small). Arrows indicate energy resource intake, whereas lines terminating with a hash indicate the energy resource is not ingested; i.e., the prey is defended. Alkaloid-defended arthropods (ants and mites, bottom row) may pass on their defenses to predators that are resistant and can sequester them, such as a poison frog (middle row, right), whereas the alkaloids may deter predation by frogs that are not resistant (middle row, left). Both resistant and non-resistant frogs also eat undefended prey, although poison frogs typically specialize on prey with defenses that provide defense against predators such as birds (top row, left). However, some predators might also sequester the alkaloid defenses (snakes, top row) They found poor support (many bootstrap values < 50%) for the expanded phylogeny, and they concluded that none of the species sequenced by McGugan et al. "…was closely related to any of the known alkaloid-producing groups." (p. 844). Heethoff et al. also argued that CO1 was inadequate for identification, given the small proportion of described mite species for which CO1 sequences were available, and that morphological identification by taxonomists would have been more efficacious. McGugan et al. (2016b) responded that their identifications were tentative, the primary goal not being the determination of mite species but rather profiling the diet of diverse populations of O. sylvatica. McGugan et al. could not precisely determine which mite species provided specific alkaloids to frogs due to their LC/MS experimental design and the number of mites sequenced (nine; the other mites were analyzed for alkaloids). Nonetheless, they were able to match specific alkaloids in skins of two frogs to those found in the mite prey of those frogs.
We will not attempt to resolve the points of disagreement between the two groups. McGugan et al.'s goal of linking arthropod and frog chemistry across geography and Heethoff et al.'s call for better taxonomic practice are both worthy and complementary. Both are essential to understand the unified, higher order ecological interactions in systems mediated by alkaloids as keystone molecules. As with TTX, alkaloids in the trophic hierarchy are expected to have collective impacts on species interactions and thus community composition.
Unfortunately, the uncharted mite biodiversity and lacunae in the details of mite chemical defense have obscured the chemistry-mediated trophic web in the frog-mite system (Vences et al. 2011). Here we move beyond phylogenetic classification of chemically defended mites to explore how phylogenetics can integrate components of biodiversity and chemical ecology with community trophic structure and species assembly. We expand the analyses of McGugan et al. and Heethoff et al. with broader taxonomic sampling, better phylogenetic methods, and new analyses to frame the evolution of acquisition of alkaloid chemical defense. More generally, we identify both the weaknesses of current approaches and knowledge gaps to be addressed for the exploration of this exciting system.

Taxonomic Identification and Sequence Alignment
Given that the CO1 gene is widely sequenced in oribatid mites and was used to identify candidate alkaloid-bearing species (Heethoff et al.; McGugan et al.), we performed a phylogenetic analysis to determine whether it is useful for evolutionary placement of mites. The CO1 sequences used in the McGugan et al. analysis were downloaded from Gen-Bank. These included nine prey samples of O. sylvatica, four additional sequences selected arbitrarily from the closest GenBank matches of the nine prey samples (by percent similarity), and the southern horseshoe crab, Tachypleus gigas, as an outgroup. We note that T. gigas is distantly related to the acariform mites (which include all species in this study) and therefore may root the tree incorrectly due to long-branch attraction artifacts (Sharma et al. 2014); a nonoribatid mite species would have been preferable.
Like McGugan et al., we identified sequences to highest percent identity using Nucleotide BLAST (BLASTn). We also used Protein BLAST (BLASTx) with the invertebrate mitochondrial genetic code to match the codon structure of CO1. McGugan et al. (2016a) considered a BLASTn similarity of > 96% sufficient to assign a sequence to species or genus and < 95% to family, apparently based on Hebert et al. (2003Hebert et al. ( , 2004. We used a more general criterion and simply identified the sequences according to the scientific name and taxonomic rank of the GenBank sequence of greatest percent similarity. We found that confident identification at the species rank of most of the mite specimens from McGugan et al. was not possible, because multiple BLASTn and BLASTx matches returned closest matches at percent similarities of < 80% and < 95%, respectively, far too low for confident taxonomic inference (BLASTx percent similarities are much higher than BLASTn because the former matches by protein and not nucleotide similarity). We also found that values within 1% of each other matched sequences from different families for both BLASTn and BLASTx; however, it may be that some of these very similar matches are due to misidentifications. We report the results of our GenBank queries in Table 1.
Sequences were aligned using ClustalW with default settings (Sievers et al. 2011)

Phylogenetic Analysis
We replicated the neighbor joining (NJ) analysis of McGugan et al. and Heethoff et al. and performed Bayesian and maximum likelihood (ML) analyses on both datasets. We partitioned the data by codon position and determined the optimal model of evolution for each. We present consensus NJ trees with bootstrap support (BS) values, the best ML trees with BS values, and majority-rule consensus Bayesian trees with posterior probability (PP) values for both the McGugan and Heethoff datasets (Fig. 2). Phylogenetic analyses were carried out on a MacBook Pro and the CIP-RES Gateway (Miller et al. 2015). All tree figures were edited using FigTree v1.4.4 (Rambaut 2019) and Adobe Illustrator v25.3.1.

Neighbor Joining Analyses
McGugan et al. used the default settings in MEGA 6 (Tamura et al. 2013), which specifies uncorrected p-distance as the model of sequence evolution (McGugan et al. 2016b). Heethoff et al. used ClustalW within MEGA 7 (Kumar et al. 2016) to align the sequences and performed a NJ analysis. They repeated the analysis using the Kimura 2-parameter (K2P) model with a gamma distribution among sites, citing the limitations of uncorrected p-distance and the common   c d e f application of K2P in barcoding (Hebert et al. 2003

Maximum Likelihood Analyses
ModelFinder (Kalyaanamoorthy et al. 2017), implemented in IQ-TREE v2.1.2 (Minh et al. 2020), was used to find the preferred model of molecular evolution for each data partition under the Bayesian Information Criterion (BIC) with subsequent merging (option: -m TESTMERGE). A three-partition scheme was optimal; the TIM3e + G4, TPM3 + F + G4, and HKY + F + G4 models were applied to codon partitions one, two, and three respectively (Chernomor et al. 2016). We ran 100 ML tree searches using 10,000 ultrafast bootstrap replicates (options: -runs 100 -B 10,000); the ultrafast algorithm reduces the risk of overestimated branch support due to polytomies or severe model violations (Hoang et al. 2018). Bootstrap values were plotted on the best tree using IQ-TREE 2. Data matrices, scripts, the log output file, and other output files from the McGugan and Heethoff analyses are available in Online Resources 3 and 4, respectively.

Bayesian Analyses
Bayesian trees were generated in MrBayes v3.2.7a (Ronquist et al. 2012). Because the models of DNA evolution used in the ML analyses could not be precisely implemented in MrBayes, each partition was sampled using a reversejump MCMC (Markov Chain Monte Carlo) procedure, which samples across all substitution models (Huelsenbeck et al. 2004). Among-site rate variation was estimated by the shape parameter of the gamma distribution (options: lset nst = mixed, rates = gamma) (Ronquist and Huelsenbeck 2003). Bayesian analyses were run twice for 2,000,000 generations using four chains, with sampling every 100 generations, and 25% of samples were discarded as burn-in following examination of chain characteristics. Trace plots and the harmonic means of the marginal likelihood values were checked to ensure convergence between independent runs. Chain swap ratios were examined to ensure that values were between 20-70%. Average effective sample sizes (ESS) for all parameters were greater than 300. No parameter values showed a discernible positive or negative trend when visualized in Tracer v1.7.1 (Rambaut et al. 2018). Stationarity of runs and chains was confirmed when the final value of the average standard deviation of split frequencies reached 0.01. Combined data matrix and script files, log output files, and other output files from the Bayesian analyses are included in Online Resources 5 and 6.

Expanded Bayesian Analysis
We augmented the Heethoff et al. dataset with newly available GenBank sequences from Galumna sp. and Scheloribates sp. (from the alkaloid-bearing families Galumnidae and Scheloribatidae) to enrich the representation of alkaloid-bearing taxa. The 22 CO1 sequences from the Heethoff et al. dataset, along with 21 sequences from Scheloribatidae (Scheloribates sp., n = 11) and Galumnidae (Galumna sp., n = 10), were aligned and trimmed as before. We report the BLASTn species identification (based on percent similarity), accession number, and higher classification in Table 2. All sequences were identified to taxa within Oribatida except KM834143, which was identified only to Sarcoptiformes. This dataset was analyzed using only Bayesian analysis, which was run twice for 3,000,000 generations with four chains, with sampling every 100 generations and 25% of samples discarded as burn-in. Average ESS values for all parameters exceeded 800, and all quality checks from the above analyses were applied. We present the majority-rule consensus tree of this "expanded tree" with PP support values (Fig. 3). A combined data matrix and script file and other output files are included in Online Resource 7.

Klimov et al. dataset
Our  (Fig. 4). A combined data matrix and script file and output files can be found in Online Resource 8.

Pachl-Schaefer dataset
The fourth analysis combined the 18S sequences used by Pachl et al. (2019) and Schaefer and Caruso (2019). Following the methods outlined above, all sequences were aligned and analyzed using IQ-TREE 2. We then pruned the matrix to include only the brachypyline taxa (59), retaining only one sequence from congeners unless preliminary analyses showed that presumed congeners did not form a clade. The pruned matrix was analyzed five times to assure consistency of results, and the tree was rooted with Poroliodes based on our preliminary analyses.
To reconstruct the evolution of alkaloid presence on the 18S phylogeny estimated from the Pachl-Schaefer dataset, we used the simmap routines from the phytools package version 0.7-90 (Revell 2012) and ape package version 5.5 (Paradis et al. 2004), under the SYM model of character change. Simmap fits a continuous-time reversible Markov model for the evolution of character states of a character and then simulates stochastic character histories using that model . We used the empirically derived symmetric instantaneous rate matrix Q under a symmetric model "SYM". This ancestral state reconstruction includes four categories of alkaloid data (Fig. 5): (1) Present: Alkaloids known to be present in the genus or family of the sequenced specimen (alkaloid presence in the particular species might be unknown).
(2) Present in Pooled Sample: Alkaloids known from a pooled sample of the genus or family of the sequenced specimen. McGugan et al. (2016a) characterized alkaloids from pooled mites from stomach contents, and although each stomach sample had several mites, it may be that the mite species sequenced was not represented in the mixed sample that was analyzed for alkaloids. Similarly, data on mite family occurrence of alkaloids from Saporito et al. (2007; Tables 1 and 2) were based on mixed species samples. (3) Absent: Alkaloids not detected in the genus or family of the sequenced specimen. (4) Not Tested: Mite taxa that were not screened for alkaloids. All alkaloid data are from Takada et al. (2005) and Saporito et al. (2007Saporito et al. ( , 2015. These characterizations are very coarse, but are useful as a framework. The data matrix, log output file, and other output files for the phylogeny in Fig. 5 are available in Online Resource 9, the taxonomic and chemical data used in Figs. 4 and 5 in Online Resource 10, and the R script and code used to generate the ancestral state reconstruction in Online Resource 10.

Phylogenetics
In all cases, the model selection method preferred three partitions over the single data partition used in the McGugan et al. and Heethoff et al.

analyses, by > 1000 BIC units (Online Resources 3 and 4).
Not surprisingly, our NJ trees recovered the same topology and similar bootstrap support values as did McGugan et al. and Heethoff et al. (Fig. 2a, b). The strongly supported clades (> 95% BS) in the McGugan and Heethoff NJ trees were recovered with 50% or greater BS support in the McGugan (Fig. 2c) and Heethoff ML trees (Fig. 2d) and strong PP support in the McGugan (Fig. 2e), Heethoff (Fig. 2f), and expanded Bayesian trees (Fig. 3).
Conflicts exist between our NJ, ML, and Bayesian topologies at poorly supported nodes. For example, poorly supported clades (< 45% bootstrap support) recovered Ceratozetidae sp. (KT947987) as sister to the three Scutoverticidae (KT947980, KT947988, and KT947985) in the McGugan (Fig. 2a) and Heethoff (Fig. 2b) NJ trees. However, in the McGugan (Fig. 2c) and Heethoff (Fig. 2d) ML trees as well as the McGugan (Fig. 2e) and Heethoff (Fig. 2f) Bayesian trees, Ceratozetidae sp. (KT947987) was recovered in clades including the two Eueremaeus species and the other ceratozetid species (KM831519). The ML (Fig. 2c) and Bayesian McGugan (Fig. 2e)   For both datasets, branch lengths in the NJ trees were about an order of magnitude smaller than those in the ML and Bayesian trees. Longer branches are generally expected in phylogenies based on evolutionary models which account for multiple substitutions, as opposed to p-distances, which do not.
The Klimov et al. Bayesian tree also shows very low support at the deeper nodes. In fact, the 95% credible set contains 113,996 trees, indicating an extremely broad set of topologies acceptable with 95% confidence. Although it was not our goal to examine mite classification, we note that, as other analyses have demonstrated, the Astigmata are nested within Oribatida (but with weak support). As well, most of the families represented by multiple sequences (e.g., Eremaeidae, Haplozetidae, Oppiidae) are not monophyletic, but support for this is very low.

Phylogenetic Distribution of Alkaloid-Bearing Taxa
The position of families with alkaloid-bearing species is of particular interest. In the Heethoff NJ tree (Fig. 2b), Parakalummidae sp. (KR105190) and Berniniella hauseri (KF293512) were recovered as sister-taxa, but with low support. However, the McGugan (Fig. 2c) and Heethoff (Fig. 2d) ML trees as well as McGugan (Fig. 2e), Heethoff (Fig. 2f), and expanded (Fig. 3)   The ten Galumna sp. specimens in the expanded tree (Fig. 3) were not recovered as monophyletic and could not be confidently united with the Galumnidae individual (KR103602) (Fig. 3). Across all seven trees in Figs. 2 and 3, the two ceratozetids (KM831519 and KT947987) were recovered as monophyletic (with weak support) only in the Heethoff ML tree (Fig. 2d). Mites identified as Eremaeidae (KT947980, KT947985, KT947988, KT947984, and KM827121) were not recovered as a clade in any tree. As both the eremaeids and ceratozetids include frog prey  The ancestral state reconstruction provides the marginal probabilities of alkaloid character states at the internal nodes. The data used to characterize each species are derived from Saporito et al. (2007Saporito et al. ( , 2015 and Takada et al. (2005). The states Present or Absent mean that alkaloids have been reported (or not) in the genus or family (Online Resource 10). Present in Pooled Sample means that alkaloids were reported in the genus or family, but the mite was from a pooled multispecies sample analyzed for alkaloids; hence, this is subject to sampling error. Arrows point to likely evolutionary origins of alkaloids. Although we scored Protoribates hakonensis (family Haplozetidae) as absent (Saporito et al. 2015), Takada et al. (2005) reported that two unspecified Protoribates screened positive for alkaloids using a crude reagent method. Higher taxonomic classifications for tip labels, provided in columns, are taken from Norton and Behan-Pelletier (2009) and Schatz et al. (2011). The tree with posterior probabilities is given in Online Resource 9 and the source of all sequences used is given in Online Resource 10 Journal of Chemical Ecology (2022) 48:384-400 394 species, it is possible that BLASTn identifications are erroneous and thus that non-monophyly in these groups is real, not an artifact of CO1 divergence.
The phylogeny from the Klimov et al. dataset (Fig. 4) reinforces the patterns described above. Almost all deeper nodes are poorly supported, and the McGugan et al. and Heethoff et al. sequences are scattered throughout the tree, except for some of the non-brachypyline species. Data from Saporito et al. (2007Saporito et al. ( , 2015 and Takada et al. (2005) on alkaloid-bearing taxa indicate that most families known to have alkaloids (Scheloribatidae, Mochlozetidae, and Parakalummidae) belong to the Oropodoidea; however, this superfamily is not monophyletic in this tree, but neither is there strong support for non-monophyly.
The 18S phylogeny of Brachypylina from the Pachl-Schaefer dataset (Fig. 5) is globally more strongly supported than the CO1 phylogenies (see Online Resources 3-9 for BS and PP values of CO1 trees), particularly at shallower nodes. In contrast to the CO1 trees, the 18S analysis provides strong support for monophyly of the Galumnoidea and Oropodoidea, both of which contain all mite taxa known to have alkaloids. The ancestral state reconstruction of alkaloids on the 18S tree suggests a minimum of two independent evolutionary origins of alkaloids in brachypyline mites within Oropodoidea and at the base of Galumnoidea, displayed with arrows in Fig. 5.

Discussion
It is clear from other systems (e.g., Groen and Whiteman [2021]) that phylogenetics once harnessed is a powerful tool to understand chemically mediated eco-evolutionary interactions in highly complex multi-trophic level systems. Knowledge of broad patterns of chemical defense evolution across both prey arthropods and frog predators will be critical for a deep time picture of how trophic interactions among cointeracting species assemblages have evolved to allow keystone molecules to accumulate (Prates et al. 2019; Fig. 1).
Our phylogenetic analyses of mites and their alkaloids have demonstrated (1) The use of inappropriate models hinders robust phylogenetic inference (which was already known), (2) CO1 performs poorly for deep phylogenetics of mites, whereas other molecules such as 18S perform well, (3) a deep understanding of the role of keystone molecules in trophic ecology is limited by the taxonomic impediment in mites and the corresponding dearth of genetic data for species identification, and (4) the lack of alkaloid characterization for the vast majority of higher taxa of mites, and (5) despite points (3) and (4), some important inferences about the origins of defensive mite alkaloids can be made. We discuss each of these below.
First, the use of inappropriate models likely accounts for the topological differences between the NJ trees and the ML and Bayesian trees. Tree construction using unrealistic models has been shown to be suboptimal (Huelsenbeck and Hillis 1993;Tamura et al. 2004). When the underlying evolutionary models are oversimplified and under-parameterized, methods are less accurate and less consistent (they converge to a wrong tree with an increased amount of data) (Felsenstein 1978;Huelsenbeck and Hillis 1993;Sullivan and Swofford 2001).
Second, we found that CO1 was not useful for inferring deeper relationships (in contrast to species identification), even with better phylogenetic methods. Most deeper nodes in all trees were poorly supported. Although there are several causes of poorly supported nodes (e.g., incomplete lineage sorting, poor taxon sampling), the most obvious seems to be the high substitution rates of CO1 as reflected in the long branches of the ML and Bayesian trees. CO1 is a fastevolving molecule in animals, and coupled with the divergence times of mites (the estimated age of the acariform ancestor is 455-552 Ma [Arabi et al. 2012;Arribas et al. 2020]), this would predict large numbers of undetected substitutions. Unfortunately, a perusal of GenBank shows that sequence data from mites in alkaloid-bearing families are largely restricted to CO1. Phylogenomic analysis of many, diverse genes capable of resolving deeper phylogenetic relationships (such as we found with 18S, see below) will be critical to reconstruct the evolution of alkaloid chemistry across Oribatida.
Due to the limitations of CO1 in inferring mite phylogeny, we question the general conclusion of Heethoff (Heethoff et al. 2016;Saporito et al. 2015).
Third, identification of the mites with chemical defenses is hindered by the unfortunate dearth of mite taxonomists, a point made by McGugan et al. (2016b). Schatz et al. (2011) reported 16,197 extant species of oribatid mites, which is on a par with the roughly 14,000 species of ants. Overall estimates of mite species richness span three orders of magnitude (0.4-114 million); their densities can be > 500,000 individuals per square meter in forest soils of temperate zones (Walter and Proctor 2013). But a quick perusal of the literature would reveal far more ant than mite taxonomists. Although pleas for greater funding of taxonomy and taxonomists are constant, the situation has not improved substantially for decades. We have no specific solution, but traditional taxonomic practice alone will not remedy this shortfall (Drew 2011).
Species identification through "barcoding" (molecular identification of unknown species) has been proposed (Hebert et al. 2003(Hebert et al. , 2004) as a panacea for the taxonomic impediment. Although the CO1 gene has been argued to play a critical role in the global effort to document biodiversity, its effectiveness and desirability as a way of overcoming the taxonomic impediment (i.e., the global gaps in taxonomic knowledge and the shortage of trained taxonomists to fill these) is questionable (Buhay 2009;Moritz and Cicero 2004;Will et al. 2005). Leaving aside for the moment the barcoding controversy, the question is whether CO1 can be used to identify mite species reliably.
Species diversity by CO1 sequences in GenBank is far from adequate for species identification. Our GenBank queries failed to identify the sequences to species, and in the best-matching records the voucher specimens often were not identified to species. For example, our GenBank search for species-level taxa of Brachypylina (the oribatid group to which the alkaloid-bearing species belong) yielded 11,106 taxonomy records at the species rank including unspecified (informal) names. Of these, 10,068 (91%) were barcode records (BOLD, BIOUG) or had an associated field number but lacked a nomenclaturally valid species name. Only 300 (2.7%) of these records were valid species names. The overwhelming majority of GenBank CO1 records cannot be used for reliable species identification, despite this being the goal of barcoding.
Fourth, although much progress has been made in the last two decades, our knowledge of the phylogenetic distribution of alkaloids among oribatid mites is rudimentary. Alkaloids have been confirmed in only a few families of the brachypyline oribatid mites as well as in one family (Uropodidae) of non-oribatids (Saporito et al. 2007(Saporito et al. , 2015. Within the brachypylines, alkaloids appear to be most ubiquitous in the Oropodoidea (Saporito et al. 2015). In particular, the alkaloid-bearing families Drymobatidae, Mochlozetidae (Saporito et al. 2007), Parakalummidae (Saporito et al. 2015), and Scheloribatidae (Saporito et al. 2015;Takada et al. 2005), are oropodoids. In addition to these, mixedspecies samples in which alkaloids have been identified comprise representatives from Oribotritiidae, Haplozetidae, Oppiidae, Trhypochthoniidae, Tectocepheidae, Epactozetidae, Hypochthoniidae, Dampfiellidae, Austrachipteriidae, and Lohmanniidae (Saporito et al. 2007;compare Table 1 of sites where mite alkaloids were detected with Table 2 of taxonomic identifications of the mites). These families are all candidates as alkaloid reservoirs.
Identification of chemically defended mites is further hampered because most of the 249 oribatid families remain to be surveyed (Raspotnig et al. 2011b;Saporito et al. 2007;Schatz et al. 2011). Because many taxa have not been surveyed, attempts to infer the presence of alkaloids requires closer scrutiny. For example, Raspotnig et al. (2011b) asserted that several mite species identified in the diets of the poisonous Cuban frogs Eleutherodactylus iberia and E. orientalis, are not possible candidates for alkaloids, since they are not members of groups that possess opisthonotal glands, specifically the "glandulate oribatids" and the Astigmata (Heethoff et al. 2016;Kuwahara 2004;Raspotnig 2010;Saporito et al. 2007;Takada et al. 2005). However, oribatid mites in the Nothrina, which Raspotnig et al. (2011b) reject as possible reservoirs, have not been screened specifically for alkaloids (Raspotnig et al. 2005;Shimano et al. 2002). Moreover, Saporito et al. (2007) detected alkaloids in a two species sample of only non-glandulate oribatids, Afronothrus incisivus (Trhypochthoniidae, Nothrina) and Eohypochthonius sp cf gracilis (Hypochthoniidae, Enarthronota), both of which are taxa presumed not to have alkaloids (Raspotnig et al. 2011b).
Also, Saporito et al. (2015) identified alkaloids in mites from the family Uropodidae, which lack opisthonotal glands. It is possible that in such mites, alkaloids are produced in secretory porose areas beneath the dermal glands (Alberti et al. 1997;Saporito et al. 2015). Further investigation into how mites without opisthonotal glands obtain and secrete chemicals is required. Perhaps mites synthesize alkaloids and also sequester alkaloids from plants and saprophytic fungi (Raspotnig et al. 2011a). The oribatids Scheloribates sp. and S. azumaensis lack alkaloids at the nymph stage but possess them as adults (Takada et al. 2005), a phenomenon that is in fact common across oribatids (Raspotnig et al. 2005). The ability to synthesize alkaloids may begin at maturity, and the nymphs of some species could lack the physiological ability to sequester the alkaloids or their precursors (Saporito et al. 2011;Schneider and Maraun 2005).
That identical alkaloids occur in oribatid mites and other arthropods raises the possibility of transfer between arthropod groups via trophic interactions (Saporito et al. 2011;Takada et al. 2005). Mutualisms exist between some oribatid mites and ants (Aoki et al. 1994;Ito and Takaku 1994), and some ants are also predators of Oribatida (Masuko 1994). Pheidole ants prey on mite taxa known to have alkaloid-bearing species, namely Galumna and Scheloribates (Wilson 2005). A second non-exclusive explanation is that a widespread alkaloid-producing microsymbiont, such as a bacterium or fungus, could account for alkaloid classes shared across arthropods (Saporito et al. 2011). Testing for the presence of alkaloid-producing microsymbionts would require culturing microbes from mites (Saporito et al. 2011(Saporito et al. , 2015. Last, we have shown that despite the lack of well-resolved phylogenies of mites and associated information on alkaloid distribution, some preliminary but important inferences about the phylogenetic origins of alkaloid defense are possible. In contrast to the CO1 trees, the 18S phylogeny is well supported; Galumnoidea, Oropodoidea, and Brachypylina were found to be monophyletic (Fig. 5), whereas they were not monophyletic in the CO1 analyses (Fig. 4). Importantly, the Galumnoidea and Oropodoidea include almost all species known to bear alkaloids (Saporito et al. 2007(Saporito et al. , 2015Takada et al. 2005). The ancestral state reconstruction points to independent origins of alkaloids in these two groups. There is possibly a third evolutionary origin of brachypyline alkaloids within Oppioidea (Oppiidae has been identified in mixed species samples positively screened for alkaloids [Saporito et al. 2007]). However, more confident determination of the origin of alkaloids for both will require denser sampling of chemically screened mites with improved taxonomic precision.
We have tried to exemplify the advantages of phylogenetics for research on sequestered defensive compounds and their effects on trophic biology and community structure (Savitzky et al. 2012), beginning with a re-examination of a controversy (Heethoff et al. 2016;McGugan et al. 2016a, b), considering newer data and analyses, and continuing with a consideration of impediments to progress in the mite-poison frog system. Despite the lack of data in some critical components of the system, it appears that alkaloid defense in mites arose more than once. This will direct future research questions. For example, are certain classes of alkaloids associated with one origin but not the other? Do certain taxa of poison frogs specialize on specific taxa of mites? How do frogs and other higher order vertebrates resist alkaloids from their specialized mites?
We have expanded on how phylogenetic systematics can bridge aspects of biodiversity and ecology to clarify chemically-mediated trophic interactions. Although phylogenetics can yield insights, solutions to the lack of data on both mites and alkaloids extend beyond phylogenetics, and will require collaborations among groups with the requisite expertise, framed in questions about the proximate and ultimate causes of the origin of chemical novelty in living systems.