Morphological and molecular analysis of Willowsia nigromaculata (Collembola, Entomobryidae, Entomobryinae) reveals a new cryptic species from the United States

Commonly reported as a household pest throughout the northern hemisphere, Willowsia nigromaculata (Lubbock) is among the most abundant and widely distributed springtails. However, taxonomic uncertainty due to incomplete morphological descriptions based on specimens from different continents may lead to incorrect identifications and/or prevent the recognition of distinct lineages within this morphospecies. Here, we perform the first comprehensive morphological and genetic comparison between W. nigromaculata specimens collected from North America and Europe. Morphological and genetic evidence reveals that populations in the United States and France represent two distinct nigromaculata-like species, but a phylogenetic analysis indicates both species may also be present in Canada. Based on these results, we redescribe W. nigromaculata from France, provide a description for Willowsia neonigromaculata sp. nov. from the United States, and propose new diagnostic characters for their separation, including the number of inner appendages on the maxillary sublobal plate. We also highlight the need for morphological and molecular investigations of additional populations to better understand the diversity and distribution of W. nigromaculata and related species.


Introduction
Originally described as Seira nigromaculata Lubbock, 1873, and later designated as the type species for the genus Willowsia Shoebothan, 1917, Willowsia nigromaculata (Lubbock, 1873) is considered to be among the most abundant and widely distributed species of Collembola, occurring in human structures such as homes, cellars, windowsills, greenhouses, and grain silos throughout the northern hemisphere (Guthrie 1903;Maynard 1951;Scott 1991;Christiansen & Bellinger 1998;Fjellberg 2007;Zhang et al. 2011;Bellinger et al. 2020). Despite being a common household and agricultural pest, the taxonomic status of this species remains uncertain.
Until recently, W. nigromaculata was identified primarily by color pattern (Gisin 1960;Stach 1967;Fjellberg 2007), especially in North America (Guthrie 1903;Maynard 1951;Christiansen & Bellinger 1998). Color pattern has been shown to be a useful species-level diagnostic character in certain cases (Soto-Adames 2002;Katz et al. 2015a;Ding et al. 2018), but other studies also highlight the importance of including chaetotaxy and other morphological characters, in addition to color pattern, for Entomobryinae species diagnosis (Jordana 2012; Katz et al. 2015b). Szeptycki (1979), Christiansen & Bellinger (1980), and Zhang et al. (2011) provided modern descriptions of dorsal chaetotaxy for populations of W. nigromaculata in Poland, United States of America (USA), and France, respectively, but many important characters were not revealed, including head chaetotaxy, mouthparts, collophore and manubrium. Katz (2017) later provided descriptions of some of these missing characters based on specimens collected in the USA but did not examine European material for comparison.
Incomplete species descriptions of W. nigromaculata based on specimens collected from different continents may prevent the recognition of distinct species, and consequently, species identifications could be misleading and/or populations may be incorrectly assumed to be introduced on other continents. In response to this concern, we conduct the first comprehensive morphological and genetic comparisons of W. nigromaculata with specimens collected from both North America and Europe to (1) determine whether morphologically and/or genetically distinct lineages are present within W. nigromaculata, (2) determine if morphological and/or genetic differences correspond to collection locality, (3) provide taxonomic redescriptions of W. nigromaculata and any potential new species and (4) reveal new diagnostic characters for the identification of any distinct lineages uncovered within W. nigromaculata.

Morphological analysis
Specimens preserved in 92-95% ethanol and identified as Willowsia nigromaculata from Paris, France and Illinois, USA, were cleared with Nesbitt's solution and then mounted on glass slides in Hoyer's medium following the procedures described by Jordana et al. (1997) in preparation for light microscopy. Specimens in ethanol gel were photographed using a stereo microscope M165C attached to a DFC420 digital camera with a dome as presented in Kawada & Buffington (2016). Photographs were digitized using Application Suite ver. 3.4.1. For scanning electron microscopy (SEM), specimens were transferred to absolute ethanol and critical point dried after sputter-coating with gold using the equipment BAL-TEC CPD 030 and BAL-TEC SPD 050, respectively. The images were made using a LEO VP 435 scanning electron microscope. The material is deposited in the Invertebrate Collection of the National Institute of Amazonian Research (INPA), Brazil, and the Illinois Natural History Survey (INHS) in Champaign, IL, United States.
Uncorrected pairwise COI distances were calculated in Geneious Prime ver. 2020.0.4 (https://www.geneious.com) and plotted in R ver. 3.6.3 (R Core Team 2020) to (1) characterize genetic variation within W. nigromaculata, (2) compare the genetic variation within W. nigromaculata with variation observed between other species of Entomobryinae, and (3) identify potential gaps in the distribution of genetic variation which may indicate the presence of both inter-and intra-specific variation within W. nigromaculata.
Bayesian and maximum likelihood approaches were employed to estimate phylogenetic relationships among Entomobryinae species to determine if distinct and divergent lineages are present within W. nigromaculata and if they correspond with differences in collection locality and/or morphology. PartitionFinder2 ver. 2.2.1 (Lanfear et al. 2017) was used to determine the appropriate models of sequence evolution and partitioning scheme. Bayesian analysis was performed using MrBayes ver. 3.2.7 (Ronquist et al. 2012) with the following parameters: codon positons 1-3 were partitioned separately with GTR+I+G, GTR+I+G, and GTR+G site models, respectively; ingroup taxa were constrained to be monophyletic; Markov chain Monte Carlo (MCMC) for 200 million generations; samplefreq = 1,000; and nruns = 2. Convergence was assessed by observation of average split frequencies values below p < 0.004. A 50% Bayesian consensus tree was generated combining both runs using the sumt command with a 25% burn-in. We also performed an additional Bayesian analysis using BEAST2 ver. 2.6.2 (Bouckaert et al. 2019) with the following parameters: codon positions 1-3 were partitioned separately for site model averaging with the bModelTest package (Bouckaert & Drummond 2017); a strict clock rate set to 1 for relative branch length estimation; Yule tree model; monophyletic constraint prior on ingroup taxa; MCMC for 200 million generations; and sampling trees and statistics every 1000 generations. After a 10% burn-in, effective sample size (ESS) for all parameters was determined to be greater than 200 with Tracer ver. 1.7.1 (Rambaut et al. 2018). A maximum clade credibility tree was estimated using TreeAnnotator ver. 2.6.0 (Bouckaert et al. 2019) with median node heights and a 10% burn-in. Maximum likelihood analysis was performed with RAxML-NG ver. 9.0 (Kozlov et al. 2019) using the all-in-one ML search and bootstrapping command (all), a partition file treating codon positions 1-3 as separate partitions with the same site models used for the MrBayes analysis, a constraint tree forcing monophyly for ingroup taxa, scaled branch lengths, and MRE-based bootstrapping (converged after 400 replicates).

Description
Length. Total length (head + trunk) of specimens 1.51-1.85 mm (n = 3). Specimens pale yellowish white with violet pigments on distal half of the Ant I to Ant IV, head laterally, margins of Th II-Abd II, Abd III with one transverse band that expands laterally, Abd IV dorsally with one irregular lateral spot and one weak central spot, irregular spots on dorsoventral plate and posterior margin; distal half of Abd V to VI, distal femur and tibiotarsus medially; eyepatches black (Fig. 8). Sometimes depigmented on proximal part of the Ant II; pigments on posterior one third of the Th II and Abd I, almost all Abd II-III medially, and posterior half of Abd IV (see Katz 2017: 551).

Genetic analysis
COI p-distances revealed remarkably high levels of genetic variation among specimens identified as W. nigromaculata, with p-distances ranging from 0% to 13.47% (Table 3). Moreover, there is a large gap in the distribution of genetic distances (0-1.18% and 9-13.47%, respectively) ( Fig. 10), corresponding to the separation of two phylogenetically distinct lineages within W. nigromaculata (Fig. 11). Distances between these two clades (Fig. 10, line C) are similar to distances between other Entomobryinae species (Fig. 10, line A) and distances between Entomobryinae sister taxa (Fig. 10, line B).
Both Bayesian and maximum likelihood phylogenetic approaches produced mostly congruent gene tree topologies, with lineages of the W. nigromaculata forming two distinct clades with high support (Fig. 11): (1) a W. neonigromaculata sp. nov. clade that includes sequences representing W. neonigromaculata sp. nov. collected from Dixon Springs, Illinois, USA and Brockville, Canada (previously identified as W. nigromaculata); and (2) a W. nigromaculata clade that includes W. nigromaculata sequences from specimens collected in Teillet, France and Toronto, Canada. Phylogenetic divergence between these two clades is similar to that observed between other Entomobryinae sister species (Fig. 11).

Discussion
Cryptic species are thought to comprise a vast majority of springtail diversity (Porco et al. 2012;Cicconardi et al. 2013), which has contributed to the rise of studies that focus on integrating morphological and molecular data to delimit cryptic species within various groups of Collembola (e.g., Zhang et al. 2014cZhang et al. , 2018Zhang et al. , 2019Katz et al. 2015a;Sun et al. 2017;Skarżyński et al. 2018;Yu et al. 2018;Carapelli et al. 2020). Likewise, we provide both morphological and molecular evidence that lineages of Willowsia nigromaculata comprise at least two species -W. neonigromaculata sp. nov. in North America and W. nigromaculata in France and Canada, but also reported throughout the Holarctic realm. Morphological analysis revealed a number of subtle, yet significant differences in chaetotaxy between the two species (Table 2), including the number of inner appendages on the maxillary sublobal plate (Figs 5E, 9B) -a character thought to be highly conserved across Entomobryinae (Katz 2017). The presence of only 2 inner appendages for W. nigromaculata from France is surprising, given that to date, only three other species of Entomobryinae are known to share this character state: Americabrya arida (Christiansen & Bellinger, 1980), Willowsia mexicana Zhang, Palacios-Vargas & Chen, 2007, and Willowsia pyrrhopygia Katz, 2017 -all of which are endemic to the New World (Katz 2017). However, until recently, this character was rarely reported in the literature, so it is possible that 2 inner appendages may be more common than previously thought.
The recovery of two highly supported and phylogenetically distinct clades (Fig. 11), separated by approximately 13% sequence divergence (Table 3; Fig. 10), provides additional support for our new species designation. However, the presence of both species in Canada indicates that the distribution of diversity in this group does not seem to correspond with geography. Extremely low levels of sequence divergence (approx. 1%) observed between W. nigromaculata specimens collected in France and Toronto, Canada (Table 3) are indicative of a recent range expansion via trans-Atlantic dispersal, possibly facilitated by W. nigromaculata's strong association with human structures. As expected, the phylogenetic analysis also supports previous observations of para-/polyphyletic genera within Entomobryinae (Zhang et al. 2014a(Zhang et al. , 2014b(Zhang et al. , 2016Katz et al. 2015bKatz et al. , 2017  Deharveng 2015; Ding et al. 2018), though the primary objective of this study was to assess genetic diversity within Willowsia, not to resolve generic relationships within Entomobryinae -a task that will likely require large, genomic data sets.
We consider the specimen identified as W. nigromaculata from Brockville, Ontario, Canada (deWaard et al. 2019), as a synonym of W. neonigromaculata sp. nov. based on genetic similarity (0.18% COI sequence divergence). Sequence divergence in COI seems to accumulate long before obvious morphological divergence in springtails, with cryptic species often having >8% sequence divergence in COI (Porco et al. 2012;Cicconardi et al. 2013;Katz et al. 2015a). Nevertheless, future morphological analysis of Canadian specimens would provide additional support for the presence of W. neonigromaculata sp. nov. in Canada.
In light of these results, we expect that specimens previously identified as W. nigromaculata, especially those collected in the USA and Canada, may actually be W. neonigromaculata sp. nov. and the reexamination of old material may reveal additional new species. Additional morphological and molecular investigations of specimens collected throughout the Holarctic region are needed to better understand the diversity and distribution of W. nigromaculata and W. neonigromaculata sp. nov. Fig. 10. COI p-distance (uncorrected) frequency histogram (in gray) and boxplots comparing p-distances between A. Entomobryinae species. B. Entomobryinae sister species. C. Willowsia nigromaculata (Lubbock, 1873) (France and Toronto, Canada) and W. neonigromaculata sp. nov (USA and Brockville, Canada) clades (inter-specific distances). D. Lineages within the same W. nigromaculata or W. neonigromaculata sp. nov. clade (intra-specific distances). Fig. 11. Bayesian consensus tree estimated with MrBayes for Entomobryinae, including the clade comprised of Willowsia nigromaculata (Lubbock, 1873) from France and Canada and W. neonigromaculata sp. nov. from the USA and Canada (highlighted). MrBayes and BEAST posterior probabilities and RAxML bootstraps are indicated at each node, respectively. The ingroup node constrained as monophyletic for all phylogenetic analyses is indicated with an asterisk (*).