Challenges associated with Neobrachinus identification led us to limit our taxon sampling to expertly identified specimens in museum collections (Suppl. material 1: table S1). We targeted material from institutions where Erwin conducted work on Neobrachinus. Specimen loans were acquired from the University of Alberta E.H. Strickland Entomological Museum (Edmonton, Alberta, CA), where Erwin deposited his vouchers after completing his PhD research with George E. Ball, culminating in several publications (Erwin 1965, 1967, 1970, 1973). We also used specimens deposited in the University of Arizona Insect Collection (Tucson, Arizona, USA), where Erwin identified Brachinus specimens to the species level as a Visiting Arthropod Systematist in 2014.

Efforts were made to sample several species from as many species groups as possible, especially within the large fumans group which we hypothesized may not be monophyletic. We also downloaded all available sequences of Neobrachinus species from the Barcode of Life Database (BOLD) and GenBank and tested their species identities against sequence data from expertly identified specimens.

Total genomic DNA was extracted from the right middle leg of specimens using the Qiagen DNeasy Blood & Tissue Kit (Valencia, CA) following the manufacturer suggested protocol. Extractions on older specimens were conducted in the Schlinger Ancient DNA Laboratory at the University of Arizona Insect Collection using the QIAamp DNA Micro Kit (Qiagen Inc., Valencia, CA) following the manufacturer suggested protocol. The concentration of total genomic DNA in extraction products was measured on a Qubit 3.0 Fluorometer (Thermo Fisher, USA). Samples with quantifiable DNA were used in subsequent PCRs.

The gene regions CAD (carbamoyl-phosphate synthetase 2, aspartate transcarbamylase, dihydroorotase) and COI (cytochrome c oxidase subunit I) have been shown to be phylogenetically informative in Neobrachinus by Schaller et al. (2018), and sequences generated in that study were downloaded from GenBank. Sequences for additional taxa to these datasets were generated following published protocols (Schaller et al. 2018).

Sequence data were also obtained for the D2-3 region of large subunit ribosomal gene (28S) from the total genomic DNA extracted for Schaller et al. (2018) as well as the new taxa added herein. We chose to add 28S to our analyses for several reasons; it has been shown to be phylogenetically informative in other genera of carabid beetles, and as a multicopy gene it is easier to amplify from older museum specimens (Sproul and Maddison 2017). Therefore, we started building a reference library of 28S sequences obtained from expertly identified specimens to facilitate future molecular work with museum material. The D2-3 region of 28S was amplified using primers LS58F and LS998R (Ober 2002) and the following PCR cycling conditions: initial denaturation at 94 °C for 2 min, followed by 35 cycles at 94 °C for 22 s, 57 °C for 22 s, 72 °C for 1 min and 10 s, and a final elongation at 72 °C for 5 min.

PCR products were quantified, normalized, and sequenced in forward and reverse directions using Sanger sequencing at the University of Arizona Genetics Core (UAGC) using an Applied Biosystems 3730 DNA Analyzer (ThermoFisher Scientific). Chromatograms were assembled into contigs, and initial base calls were made using Phred (Green and Ewing 2002) and Phrap (Green 1999) as implemented by the Chromaseq 1.52 module (Maddison and Maddison 2020) within Mesquite 3.61 (Maddison and Maddison 2019). Final base calls were made by visual inspection of the contigs.

Three single gene matrices (COI, CAD, and 28S) were assembled. Each matrix contained sequences generated specifically for this study as well as all homologous sequences of Neobrachinus publicly available on BOLD and GenBank (databases searched January 2021) (Suppl. material 1: table S1). The COI matrix contained 270 taxa, the CAD matrix contained 70 taxa, and the 28S matrix contained 54 taxa. Sequences in each matrix were aligned using default settings in MAFFT v. 7.474 (Katoh and Standley 2013) within Mesquite. A concatenated matrix with the data from all three genes was also assembled, which contained 282 taxa including 228 Neobrachinus (representing 9/15 Neobrachinus species groups, and 32/62 Neobrachinus species) and 54 outgroups. In the concatenated matrix and single gene matrices, COI and CAD characters were partitioned by codon position. Maximum-likelihood analyses and bootstrap analyses were conducted on single gene matrices and on the concatenated dataset using IQ-TREE v. 1.6.10 (Nguyen et al. 2015). The ModelFinder feature within IQ-TREE (Kalyaanamoorthy et al. 2017) was used to find the optimal character evolution models. The ModelFinder Plus model option was used for 28S, and the TESTMERGE option for the protein-coding genes and for the concatenated dataset. One hundred searches were conducted for the maximum-likelihood tree for each matrix in Mesquite. Bootstrap analyses for the four trees were conducted with 1000 replicates using IQ-TREE v. 1.6.10 (Nguyen et al. 2015), as orchestrated by the CIPRES Science Gateway (Miller et al. 2010). Support for and against clades were calculated for each species group and the subgenus Neobrachinus in Mesquite using its “Clade Frequencies in Trees” feature and bootstrap trees generated by IQ-TREE (Maddison et al. 2019; Kavanaugh et al. 2021).

IQ-Tree ModelFinder and ModelFinder Plus identified the following models of evolution for each character partition: COI codon 1 = K2P+I+G4; COI codon position 2, CAD codon 1, CAD codon 2 = 2K2P+I; COI codon 3 = TIM2+F+G4; CAD codon 3 = HKY+F; and 28S = GTR+F+I+G4.

The three-gene IQ-Tree analysis resulted in the phylogeny shown in Figs 3–7. The full concatenated tree, and individual gene trees for 28S, CAD, and COI, are shown in Suppl. material 1: figs S1–S4). The COI analysis identified several specimens obtained from public databases that could be misidentified (Suppl. material 1: fig. S4). In all analyses, the following species groups originally proposed by Erwin 1970 were recovered as monophyletic: cordicollis, lateralis, and texanus.

Figure 3. 

The maximum-likelihood three-gene molecular phylogeny of Neobrachinus with clades collapsed. Clades are colored by species group. Clades in solid blue were formerly placed within the fumans species group. Nodes with bootstrap values > 0.90 are denoted with grey circles. Clades present in South America are denoted with “SA.”

Figure 4. 

Phylogeny within the cinctipennis and phaeocerus species groups.

In all analyses the fumans species group proposed by Erwin (1970) was polyphyletic (Fig. 3, Suppl. material 1: figs S1–S4). B. medius fell within a highly supported clade containing Erwin’s americanus and cordicollis groups (Fig. 6, Table 1). Support for and against the subgenus Neobrachinus and each species group are shown in Fig. 8. A revised classification of Neobrachinus species that reflects molecular and morphological support for species groups is shown in Table 1.

This study used molecular data to test previous hypotheses of species group membership and phylogenetic relationships in the subgenus Neobrachinus that were proposed based on morphological data (Erwin 1970). The majority of Erwin’s species groups were recovered as monophyletic and were supported with high bootstrap values in all analyses (Figs 3–7). However, both molecular and morphological evidence support splitting the fumans species group into new species groups. Furthermore, relationships between species groups, for the most part, remain unclear.

The shape of the virga was not found to be phylogenetically informative as envisioned by Erwin (1970). For example, the revised cordicollis group now encompasses members of the former americanus group, as well as B. medius which was previously placed in the fumans group (Fig. 5). The virga of the americanus group was considered plesiomorphic among Neobrachinus, while the “H-shaped” virga of the cordicollis group was considered highly derived. Erwin (1970) also hypothesized that speciation within the subgenus Neobrachinus was largely connected to the evolution of the virga. Although molecular data largely corroborated species groups that were diagnosed by morphological characters, including the form of the virga, the polyphyly of Erwin’s fumans group indicates that molecular data is necessary to confirm the monophyly of species groups within Neobrachinus.

Figure 5. 

Phylogeny within the fumans species group, B. fumans colored blue and B. perplexus colored orange to highlight cryptic diversity and/or potential misidentifications of sequences on public databases.

Erwin’s fumans group contained 26 morphologically diverse species and was defined by a troughed virga. Considering the molecular evidence that supports splitting the fumans group, the troughed virga could be an ancestral or convergent form among Neobrachinus.

Species subgroups of the fumans group were also polyphyletic in the molecular phylogeny, highlighting potential convergent character states of the male genitalia. Members of Erwin’s quadripennis subgroup of the fumans group were recovered throughout the Neobrachinus tree: in the phaeocerus species group (Fig. 4), the quadripennis species group (Fig. 3), and the mexicanus species group (Fig. 7). This group was characterized by a ridge on the ventral surface of the male genitalia. The two members of Erwin’s gebhardis subgroup, B. gebhardis and B. galactoderus, were recovered in separate clades (Figs 3, 7); this group was characterized by the form of the median lobe of the male genitalia and the restriction of elytral pubescence to outer edge in the eighth depression.

Erwin’s fumans species group also contained seven monotypic species subgroups, of which four were included in this study. Two of these, B. cyanipennis and B. ovipennis, formed a clade and are now placed together in the cyanipennis species group (Fig. 7). As previously mentioned, B. medius had strong molecular support for its new placement in the revised cordicollis species group. Finally, B. tenuicollis remains a monotypic species group, as proposed by Erwin, with the support of molecular and morphological data (Fig. 6). All species of the fumans group that were in monotypic species groups not included in the molecular study, and all taxa not included in Erwin (1970) are considered incertae sedis (Table 1).

Figure 6. 

Phylogeny within the tenuicollis, americanus, medius, and cordicollis species groups. Novel additions to the cordicollis group are color-coded: americanus group (green) and B. medius (blue). Misidentifications of sequences on public databases are colored red.

Figure 7. 

Phylogeny within the mexicanus, cyanipennis, and gebhardis species groups. B. kavanaughi colored red to highlight potential misidentifications of sequences or the need to synonymize this species with B. mexicanus.

Erwin postulated all 84 species of Neobrachinus evolved from a single most recent common ancestor that crossed the Bering Land Bridge. The molecular phylogeny supports a Nearctic origin of the Neobrachinus, as predicted by Erwin (1970). Two clades, the lateralis and texanus species groups, have members that are present in South America (Figs 3–7). All other clades of Neobrachinus are only present in the Nearctic (Figs 3–7). He also hypothesized that the South American Neobrachinus species diversified from a single colonization event by an ancestral Neobrachinus lineage, giving rise to a monophyletic group containing the brunneus, grandis, lateralis, and texanus species groups (1970). However, the molecular phylogeny inferred in this study indicates otherwise. The two clades represented in this study with membership in South America, the lateralis and texanus species groups, are not sister taxa (Figs 3–7), indicating that multiple colonization events to South America must have occurred. Inclusion of additional South American taxa in a molecular phylogeny (Table 1) would illuminate their biogeographic history.

Among the previously published sequences downloaded from public databases, molecular phylogenetic analysis revealed several cases where specimens were likely misidentified. Some specimens from North Dakota were identified as B. medius (BETN1837-18, BETN9260-20, BETN9117-20, BETN9121-20), however the sequences were in a well-supported clade, separate from B. medius from the same region (Fig. 6, Suppl. material 1: fig. S1).

Other potential misidentifications exist yet are difficult to confirm. For example, within the new fumans species group, there are several clades of the species B. fumans and B. perplexus (Fig. 5). Given the lack of molecular sequence data from other members of the species group, it is impossible to determine whether these clades represent cryptic diversity or whether the specimens are misidentified members of other species of the fumans group.

Another example is the clade containing B. kavanaughi and B. mexicanus (Fig. 7). Representatives of the species B. kavanaughi do not form a separate clade from B. mexicanus. Without examining the voucher specimens it is impossible to determine whether these specimens represent misidentified members of B. mexicanus, or whether B. kavanaughi should be synonymized with B. mexicanus.

Figure 8. 

Bootstrap support for and against clades of Neobrachinus. Each column has maximum likelihood bootstrap values as percentages for or against each clade recovered in each dataset: the three-gene concatenated matrix (3G), and the single-gene datasets, 28S, CAD, and COI. Positive values indicate support while negative numbers indicate support for the contradictory clade with the highest support. Cells with bootstrap values ≥ 90 are in black, values between 75 and 89 in dark grey, and values between 50 and 74 in light grey. Cells in red have bootstrap values for the contradictory clade ≥ 50.