Adult terrestrial beetle specimens were collected between 2006 and 2010 at 223 sites in the Churchill region (Figure 1) primarily using pitfall, coloured pan, and/or Malaise traps. Berlese extractions of the substrate and/or active netting occurred at some sites. Both adults and larvae of aquatic species were collected using dip nets and bottle traps from a variety of freshwater habitat types: coastal rock bluff pools, tundra ponds, fen ponds, streams, lakes, and the Churchill River. Sampling was concentrated in areas accessible by road in Churchill, bounded by the Churchill River in the west, Hudson Bay in the north, a boundary approximating 58.6°N in the south, and a boundary approximating 93.4°W in the east. A small number of specimens were also collected at 11 sites outside these boundaries, including one site at Bylot Station (58.42°N, 94.13°W) and 10 sites in Wapusk National Park, including five sites near Nestor 1 field camp (58.66°N, 93.19°W), four sites near the mouth of the Broad River (58.42°N, 92.87°W), and one isolated relict beach ridge (57.58°N, 93.87°W).

All specimens were either pinned and air dried, or preserved in 95% ethanol. Specimens from earlier years were held at room temperature following collection, but aquatic specimens collected in 2010 were held in ethanol and at −20°C following the field season. Specimens are deposited at the Biodiversity Institute of Ontario, University of Guelph, Guelph, Ontario, Canada (BIOUG); the Wallis-Roughley Museum of Entomology, University of Manitoba, Winnipeg, Manitoba, Canada (WRME); and the Canadian National Collection of Insects, Arachnids and Nematodes, Ottawa, Ontario, Canada (CNC). Locality data, specimen photographs, voucher numbers, and identifier for all specimens are available through the Barcode of Life Data Systems (BOLD; http://www.boldsystems.org) 23 under project names: Coleoptera of Churchill [COCHU] and Aquatic Coleoptera of Churchill 2010 [EBCCH]. All specimens analyzed together for this study can be accessed as a single BOLD dataset dx.doi.org/10.5883/DS-COLCHU. Also see Additional file 1 for a list of specimens, taxonomy as of October 22, 2013, and locality data.

Generally, from each site, approximately 2–5 individuals per gross morphospecies within each family were selected for DNA barcoding. For a detailed freshwater study in 2010 including both adults and larvae, at least 10 individuals per morphospecies across sites were selected by EB, whenever available. Additional specimens were selected from groups that were difficult to distinguish or which possessed many individuals in larval form (e.g. the genus Hydroporus and family Gyrinidae). For some terrestrial and aquatic families, some of the selected individuals were identified to species level on the basis of morphology by a number of experts, with the identifier listed in the specimen record for each individual specimen on BOLD. Other individuals were identified to species level on the basis of genetic matches (see below) to other specimens in our dataset or to other records on BOLD. These are designated with “Digital Photograph; BOLD Barcode Library (Date)” in the Identification Method column; these cases were verified to be reasonable matches using the photographs. For families lacking available taxonomic expertise, specimens were identified to family or genus level by TW or EB 24 25 , adjusted to the taxonomic hierarchy of Bouchard et al. 26 , and a provisional interim species name was assigned based upon the DNA barcode results (see below) to enable provisional determination of species richness. Molecular Operational Taxonomic Unit (MOTU) names were assigned using sequential codes within genera, starting at CHU1 for Churchill. For MOTUs lacking genus-level identifications, a placeholder genus name (family name plus the initials TW) was inserted into the genus field in BOLD, to permit the distance analysis calculations to be performed on the MOTUs. Similarly, single morphospecies containing two or more provisional species based on the barcode results were assigned CHU codes onto the end of the binomial Linnaean name, with barcode clusters treated as separate species for richness analysis. This dataset and the associated records in BOLD will continue to gain species-level identifications over time; all identifications and interim names as analyzed in this paper are available in Additional file 1.

Tissue samples consisting of one leg (occasionally two or more for small-bodied species) were removed from specimens and deposited into 96-well plates prefilled with 30 μl of 95% ethanol. All instruments used to remove leg tissues were cleaned in 95% ethanol and sterilized by flame between specimens. Most plates were submitted to the Canadian Centre for DNA Barcoding (CCDB) where DNA was extracted from tissue samples following standard invertebrate protocols 27 28 . Molecular processing of the 2010 aquatic samples (919 specimens) was performed using the manual version of the same protocols. The barcode region of cytochrome c oxidase subunit I (COI) was amplified using a polymerase chain reaction (PCR). Typically, most specimens were subjected to a second attempt at PCR amplification if the first attempt did not yield a full-length barcode sequence of at least 500 base pairs (bp) (see Additional file 2 for primer details). PCR amplification conditions, product checking, PCR cycle sequencing, and bidirectional sequencing followed standard protocols employed at the CCDB 29 30 31 . Forward and reverse sequences for each specimen were combined into a single consensus sequence and aligned using the CodonCode Aligner software v. 3.0.2 (CodonCode Corporation). Sequences and alignments were verified to be free of gaps and stop codons using the amino acid translation. The alignment was uploaded to BOLD and also imported into MEGA version 5.0 32 for data analysis. Sequences, PCR and sequencing primers, and trace files are available for each individual specimen through the two projects on BOLD listed above (primer sequences and references in Additional file 2).

Sequences were verified as belonging to Coleoptera using Genbank BLAST and using data scrutiny tools available through BOLD (BIN discordance report, ID engine, photographic data, and NJ trees). Only high-quality sequences with a minimum length of 500 bp and containing fewer than 1% missing nucleotides (Ns) were retained for analysis of genetic divergence patterns to reduce intraspecific variations due to sequence length variation 23 . Sequences of at least 300 bp length were retained for providing provisional species identifications, for analyzing habitat occupancy, and for estimation of total biodiversity. Shorter sequences still provide reliable matches to conspecifics 33 and for biodiversity analysis it is beneficial to include the maximum data available.

To visualize the barcode data, a neighbour-joining (NJ) tree 34 was built in MEGA 5.0 32 using the following parameters: Kimura 2-Parameter (K2P) distance model 35 with pairwise deletion of gaps/missing data and inclusion of all substitutions (transitions and transversions). Although p-distances have recently been advocated for DNA barcoding studies, at low taxonomic levels K2P distances tend to be similar to p-distances 36 37 ; thus, we use the more commonly applied former metric to facilitate comparison with other works. Sequences of at least 300 bp length were included in the NJ tree, but eleven short sequences were omitted because of lack of overlapping regions with other short sequences. These sequences were selected so as to minimize the total number of sequences omitted while also ensuring that no species was entirely deleted. A NJ analysis including all families was first performed to verify separation of the individuals morphologically assigned to different families. In order to facilitate bootstrap analysis (1000 replicates), NJ phenograms were subsequently constructed on four data partitions separately: the three largest families individually (Carabidae, Dytiscidae, Staphylinidae) and all remaining families together (see Additional file 3).

For families lacking species-level identifications, and when single morphospecies formed two or more clusters, we inspected the genetic distance matrices obtained through MEGA in order to assign provisional species codes (Additional file 4). Our Molecular Operational Taxonomic Units (MOTUs) were firstly defined based upon Barcode Index Numbers (BINs) 38 , accessed through BOLD on Nov 21, 2012. BINs are genetic groupings assigned by BOLD3 for sequences that are least 500 bp in length. The BIN algorithm is based on a 2.2% initial seed sequence divergence that is combined with an algorithm that permits higher or lower divergences within BINs on the basis of genetic distance patterns 38 . We assigned our shorter sequences to these MOTUs if they clustered within a particular BIN, and we separated clusters lacking BINs (due to sequences being <500 bp) when they showed an average divergence of 3% or greater. Although threshold-based methods can be problematic for large spatial scales 22 and with comprehensive global taxonomic sampling 39 , the early-proposed 2% threshold to separate intraspecific and interspecific divergences 40 is effective at discriminating most insect species within the Churchill region e.g. 18 19 . This value is slightly conservative in some insects, in that it underestimates the number of species in Diptera, which are hyper-diverse in Churchill 41 . We conservatively employed a 3% threshold, as short sequences may lie in more variable regions of the barcode sequence compared to the average for the entire barcode region. Therefore, our MOTU richness estimates are expected to represent minimum estimates of species-level richness. Closely related species that were morphologically identified yet show low genetic divergence were retained as separate species. While these may, in fact, represent single species with higher levels of morphological variation than known previously, the conservative approach dictates that there is insufficient evidence to overrule expert morphological identification.

Using the ≥500 bp data set, the maximum genetic distance among individuals within provisional species was calculated for all species having a sample size of at least two. Nearest interspecific distances were calculated for each species, as being the genetic distance to the nearest neighbouring sequence of a different species. These distances, both obtained using the Barcode Gap Analysis function in BOLD3, were then plotted on a histogram using R 42 .

Based on recorded field data, each specimen was assigned to one of 22 habitat classifications, which covered both aquatic and terrestrial habitats (Additional file 5). The completeness of our survey was investigated by constructing randomized accumulation curves of provisional species for: all species together, species belonging to the three most dominant families separately (Carabidae, Dytiscidae, Staphylinidae), and for four dominant habitat types separately. The habitats include Dry Forest, Wet Forest, Open Wet, and Open Dry habitat (i.e. dry tundra, thinly vegetated or unvegetated gravel or sand). Wet habitats include both aquatic habitats (ponds, lakes, streams) and wet terrestrial and semi-terrestrial habitats such as riparian zones, fens, and bogs. All other habitat types were designated as dry. Forested habitats were those for which significant tree cover was reported, and included ponds and riparian zones in forested areas (Forested Wet). Species accumulation curves were built with 1000 permutations and sub-sampling without replacement 43 in the R package Vegan 42 44 . The Chao1 45 biodiversity estimator was calculated for the entire collection (all 3203 individuals with sequences of ≥300 bp) using EstimateS 8.2 46 .

Geographic ranges of Carabidae and Dytiscidae species, two of the most speciose families in this study, were reviewed 2 3 4 7 12 47 48 to investigate possible sources of colonization (Beringia, southern open ground, or forested habitats) that could provide insight into present distribution patterns in the Churchill region. It is expected that colonization in open-ground habitats proceeded from Beringia or Arctic coastal refugia (leading to modern distributions west of Hudson Bay and north of the tree line) and from the south following glacial retreat (leading to modern distributions that include barren ground parts of Quebec and Labrador). Forest dwellers would be expected to have modern transcontinental distributions 7 49 . For each named carabid and dytiscid species with at least five specimens, a crosstab (chi-square) analysis was performed (PASW Statistics 19, SPSS, Inc., 2010, Chicago, IL, http://www.spss.com). These analyses test for differences in distribution of individuals among the four dominant habitat categories from what would be expected if the specimens were randomly distributed, based on total specimens collected from each habitat as an approximation of sampling effort.

Biodiversity of Coleoptera shows a strong tendency to decline at northern latitudes. Although our study led to a substantial increase in the known fauna from Churchill, from ca. 60 to >300 species and provisional species, and from 20 to 33 families, our results based upon six years of collecting and DNA barcoding of >3000 specimens nevertheless confirm the relatively impoverished fauna compared with that of more southerly locations 2 5 6 7 . Moreover, there is a relative paucity of phytophagous and soil-dwelling detritivorous species. Predatory families such as Carabidae, Dytiscidae, and Staphylinidae are the most abundant and diverse groups in both terrestrial and aquatic habitats, and accumulation curves suggest that the latter is even more diverse than reported here. Danks 3 suggested that northern ecosystems are more strongly regulated by abiotic conditions, such that resources such as plant litter and living plant tissue are incompletely exploited. Moreover, the diversity of food sources for phytophages also decreases with latitude 50 . Phytophagous and saprophagous Coleoptera, represented in this study by the families Chrysomelidae, Elateridae, Cerambycidae, Curculionidae, Cryptophagidae, and Leiodidae, may be largely unable to adapt to the soil and climatic conditions as well as lower plant diversity 51 and are thus proportionally less represented in Churchill in comparison with their overall Canadian diversity 2 7 .

In the analysis of carabid and dytiscid habitat distributions, approximately 76% of the specimens in these families were collected in Open Wet habitats, particularly bluff and tundra pools. Dytiscids were, as expected, almost completely absent from dry sites, while carabids in wet habitats were not collected from truly aquatic, submerged habitats. Among named species with five or more records, only two species of Dytiscidae (Ilybius erichsoni, Neoscutopterus hornii) showed a significant preference for Forested Wet habitats, which is in agreement with published modern distributions 52 . Thirteen species were found most frequently in Open Wet habitats (Table 2). The remaining 22 named species showed no significant deviation from random, suggesting that aquatic habitats in both forested and open areas were suitable. Within the Carabidae, only Notiophilus aquaticus showed a significant bias toward Forested Wet habitats, which is in agreement with published modern distributions 7 53 . Elaphrus americanus preferred Open Wet habitats, and Pterostichus caribou preferred both Open Wet and Open Dry habitats. Of the remaining species with five or more records, eight species showed a significant preference for Open Dry habitats, while nine species showed no preference.

Only those species with at least five specimens collected are considered. Canadian provincial records summarized in Bousquet et al. 2 ; other available distributional references are supplied below.

*distributional classes from Danks 7 A = Arctic, southern boundary north of the tree line; SA = southern Arctic, transcontinental north of the tree line; SAxB = southern Arctic excluding Beringia, transcontinental north of the tree line; WC = western and central, range from western North America east to Hudson Bay; S = southern, northern boundary at the tree line; SHA = southern high arctic, arctic species extending south of the tree line; SLA = southern low arctic, southern species extending north of the tree line.

** all Dytiscidae were collected in aquatic habitats.

A steep latitudinal diversity gradient is seen in the Lepidoptera [deWaard et al. unpubl.], although both the Hymenoptera 21 55 and Diptera 41 , Wang et al. unpubl.] are highly speciose in Churchill. However, for these latter two groups, it is currently challenging to assess the strength of their latitudinal diversity gradient due to the lack of studies in multiple regions employing genetic methods. DNA barcoding has often revealed higher diversity than traditional taxonomy, particularly in parasitoid groups 55 56 57 58 . Nevertheless, the striking diversity of Hymenoptera and Diptera in the north suggests that these groups may have radiated into some of the niches occupied by beetles at more southerly latitudes. Intriguingly, the Hymenoptera fauna of Churchill is highly dominated by parasitoid species 21 , which may mirror our detected pattern of a greater proportional representation of predator beetle families.

The current insect assemblage in the Churchill region is relatively young, incorporating species from several different refugia that colonized following the LGM, most notably species from Beringia, which followed the deglaciation of the Arctic coast and then south to Churchill, and species from the south that followed the ice margins northward. Dispersal ability of insects and their colonization patterns following the LGM may provide information on their potential responses to current climate warming 13 . The ongoing warming of the Churchill River resulting from damming and climate change may enhance its role as a corridor for northwards dispersal of both the terrestrial and aquatic beetles to the region 13 . Additional sampling of beetles in some habitats, particularly dry boreal forest, as well as further barcoding of beetles from other regions of Canada are both necessary to complete a phylogeographic analysis of likely species colonization routes following the LGM. Spitzer and Danks 59 also suggest that boreal forest peatlands, with their high habitat diversity and potential as refugia and habitat islands, could harbour endemic and relict species and contain significant and undocumented richness.

Danks 3 notes that 51% of boreal insect species have a transcontinental distribution, but only 4% of tundra/open ground species do. Other sources, however, indicate that many species have distributions that include habitats north of the tree line east of Hudson Bay in Quebec, Labrador, and sometimes Newfoundland, suggesting that retreating ice was followed by recolonization from the south, rather than assuming dispersal across Hudson Bay 2 52 . Danks 3 also points out that warming occurred prior to significant retreat of the ice margins, and many species persisting in southern refugia may not have survived that increase in temperature. Ten of the 13 dytiscid species that we found to have a preference for aquatic habitats in Open areas (Table 2) are reported as occurring only west of Hudson Bay by Bousquet et al. 2 , suggesting the difficulty of dispersal across Hudson Bay, although Larson 52 records many dytiscid species as transcontinental. Exceptions to this pattern are Hydroporus fuscipennis, which occurs across Canada in boreal and grassland lakes and ponds, and H. notabilis, which occurs in the low arctic across Canada, preferring peat pools 2 52 .

Distribution patterns of Carabidae species are more complex. Of the 11 species that showed a significant habitat preference (Table 2), nine prefer Open Dry habitat. Most of these species are currently distributed north of the tree line on open tundra 2 12 47 53 60 , with fossil records from areas in Yukon, Alaska, and other regions that had a similar dry tundra habitat during the LGM 48 61 . Carabus taedatus is a notable exception, with a fossil record from the Northwest Territories 48 but not Beringia 53 . Amara quenseli may be found throughout Canada, typically in open grassy habitat 2 53 , with a fossil record in Siberia that suggests a Beringian origin 48 . Notiophilus aquaticus, with a preference for Forested Wet habitat in this study, is circumpolar in distribution and found throughout Canada 2 47 53 . It is reported by Lindroth 47 as occurring on open ground, although Danks 7 indicates that it is found in both forested and tundra habitats. Garry 10 documented 65 species of Carabidae northwest of Churchill in the drainages of the Caribou and Seal Rivers (53 woodland, 34 tundra, 22 common to both major habitat types). Many species are well represented in midwestern fossil assemblages and possibly followed the retreating glaciers in appropriate habitats, in contrast with Danks’s 3 view that many Churchill species represent Beringian fauna. Garry 10 notes, however, that for at least 25 of the species he recovered, the glacial refugium and route(s) of recolonization are not clear.

Elaphrus americanus, with a preference for Open Wet habitat in this study, occurs throughout much of North America in moist habitats 2 47 53 . During the LGM it persisted in both Beringia and south of the ice sheet, with some morphological divergence of these refugial lineages observed 6 62 . Our results show a maximum intraspecific divergence of 1.4% in this species, which could indicate that multiple refugial lineages may be present in Churchill. Other named species with sufficient intraspecific divergence to suggest multiple refugial lineages or multiple species include Simplocaria metallica (Byrrhidae), Cymindis unicolor (Carabidae), Elleschus ephippiatus (Curculionidae), Agabus antennatus, A. bifarius, Hygrotus novemlineatus (Dytiscidae), Sericus incongruus (Elateridae), Gyrinus dubius (Gyrinidae), Hydrobius fuscipes (Hydrophilidae), and Philonthus boreas (Staphylinidae).

Globally, Coleoptera are the largest order of insects in terms of described species diversity 1 . Despite "the Creator’s inordinate fondness for beetles" 63 , Coleoptera have not been favoured to date by barcoders. For example, using the public data portal available through BOLD3 (accessed June 24, 2013), there were ca. 82 K public barcode records for beetles representing ca. 19 K provisional species. By comparison, the others of the top four most diverse insect orders are represented by approximately 2.5-fold (Hymenoptera), 3.5-fold (Diptera), and 8.5-fold (Lepidoptera) more public records. Moreover, several important studies on genetic variability within and between Coleoptera species have largely employed genetic regions other than the standard animal barcode region e.g. 22 64 . Thus, the DNA barcoding of Coleoptera is in its infancy, especially when considering their described 1 and projected 65 global diversity.

Although our study contributes important regional-scale data (and ~4% of the total) to the pool of barcode data for Coleoptera, further sampling from other regions is required to be able to use these data to “test” DNA barcoding for this taxon. Using a different region of the COI gene, Bergsten et al. 22 determined that the success rate for identifying individuals of a group of aquatic beetles (tribe Agabini, family Dytiscidae) to morphological species was nearly 100% at small spatial scales. However, for this taxon, this diminished to ~90% at continental scales, or lower, depending upon the match criteria employed. Further assessment across numerous beetle families will be required to further evaluate how well COI variability delineates species boundaries in the Coleoptera at broad spatial scales. Additionally, integrative methods incorporating genetic, morphological, and ecological information e.g. 56 57 58 may be needed, rather than assuming that morphospecies represent true species boundaries against which to test DNA barcodes. Our data are useful for such future tests, as the locality and sequence data are available online, and vouchers are retained for all barcoded specimens and are housed in publically accessible collections.

Being confined to the Churchill region, our study was originally expected to be a zone of a large “barcoding gap” 39 , in which there is a discontinuity between intraspecific vs. interspecific divergences for two main reasons. First, the total richness of Coleoptera species is less than that in more southerly regions. Second, there may have been pruning of intraspecific genetic variation by glaciations 66 67 . Our results, which are based upon a mixture of named species and provisional species, indeed do support a general difference in intraspecific vs. interspecific divergences for most species among the Churchill beetles as well as a projected high rate of success in future barcode-based identifications for beetles of this region. Of the 184 named species, only 4 (2.2%) displayed sharing of MOTUs under the BIN definition, and all others were readily separated into one or more MOTUs. However, we observed major overlap in the total divergence ranges (maximum intraspecific up to 4.1% and minimum interspecific down to 1.1%). Thus, this high rate of success at separating species reflects the general effectiveness of the BIN algorithm in recovering clusters/species of Churchill beetles. Meanwhile, new specimens belonging to the two pairs of closely related species could still be identified on the basis of their barcodes, since both members of these species pairs are now in the data set and they do show some genetic differentiation. Interestingly, the family Muscidae (Diptera) similarly does not show a clear barcoding gap 41 , despite the extremely high effectiveness of barcode-based species-level identification in the muscids using clustering methods. Despite a lack of a universal barcoding threshold, the strong species-level clustering pattern observed here and in other Churchill invertebrates points towards success for future barcode-based studies of beetles in the Churchill region and other northern regions. Well-populated datasets do not require a strict threshold to hold true for identification success to be achieved 39 41 . Thus, near-comprehensive surveys, such as being conducted for the Churchill region, are a valuable resource for the further study of biodiversity.