Field collection

We collected Trichoptera larval specimens from the subarctic location of Churchill, MB, Canada from June 5 to August 25, 2010. This research, which did not involve any human or other vertebrate individuals, embryos, or tissues, was conducted under a permit (WB11245) issued by the Manitoba Conservation Wildlife and Ecosystem Protection Branch from Winnipeg, MB to the Churchill Northern Studies Centre (CNSC) for conducting research in the Churchill Wildlife Management Area. We sampled specimens from a variety of freshwater habitats, including 30 rocky coastal bluff ponds, 30 tundra ponds, seven creeks, five lakes, and three points along the Churchill River (see Table A in S1 Table). Using a 250 μm dip net and hand picking, we sampled each location on three dates, once per month in approximately the same order. We defined a local community of potentially interacting individuals as the entire aquatic habitat encompassing both the benthic and pelagic regions. We standardized sampling effort across sites by the area sampled. For small ponds (<20 m of shoreline), we sampled the entire shoreline using the same protocol as for large habitats, executed to the degree possible given the nature of the habitat. If the habitat was large (i.e. a lake), then 20 m of shoreline was sampled. For large habitats we sampled a 20 m transect parallel to the shore, collecting 1 m away from the shore and then again 5 m from the shore or until a depth of 1.5 m. To ensure that all the species occurring in the location were collected, we performed the sweep along the transect twice in succession and retained a minimum of 10 individuals per field morphospecies at each site on each sampling date. Field morphospecies were delimited on site as individuals with similar size, colour, markings, and case material; these morphospecies were later validated through microscopic and genetic analysis. A sub-set of sites was sampled three times during one monthly visit; a comparison of genetically confirmed species accumulation curves for these sites verified that two sweeps accurately captured the local biodiversity [25]. We preserved specimens in 95% ethanol and upon return from the field stored them at -20^o C.

Molecular analysis and species identification

We sequenced portions of one mitochondrial (COI) and three nuclear genes: cadherin (CAD), elongation factor 1 alpha (EF1-α), and RNA polymerase II (POL-II). These genes are commonly used for phylogeny reconstruction in Trichoptera [20,26,27], and we found high polymerase chain reaction (PCR) success based upon a pilot study including 2 other candidate gene regions (28S and Wingless). We sorted specimens to family based on Wiggins [28] and then, when available, selected 10 individuals per morphospecies (delineated using a microscope in the lab using the characters discussed in [23]) across all sites and samples for sequencing COI. We removed one leg from each specimen that we selected for sequencing and extracted the DNA using a standard, high-throughput invertebrate DNA extraction protocol [29]. We increased the initial DNA elution volume of 50 μl for COI to 100 μl for the nuclear genes as this increased amplification success.

To amplify the barcode region of COI, we performed PCR using standard DNA barcoding methods and a Lep/Folmer primer cocktail [30] (see Table B in S1 Table [27,31–34]). For the nuclear genes (CAD, EF1-α, and POL-II) (see specimen selection criteria below), we used a 25 μl reaction consisting of the same proportion of trehalose, 10x buffer, MgCl₂, dNTP, and Platinum Taq polymerase as used for COI, but we increased the amount of forward and reverse primer to 1.25 μl of 10 μM and added 2 μl of DNA template. Successful PCR products were bidirectionally sequenced using protocols outlined by the Canadian Centre for DNA Barcoding [35] (see Table C in S1 Table).

Using Codon Code Aligner (Codon Code Corporation, v. 3.7.1), we edited and concatenated the forward and reverse chromatograms. We then subsequently aligned sequences using Clustal W and checked the amino acid sequences for all genes for stop codons and indels in MEGA 5.0 [36]. We uploaded sequences to BOLD and identified the specimens to the species level using only the expert-identified published database of COI sequences available on BOLD that had ≥98% sequence similarity to our sequences [21,23]. We employed a 2% threshold for sequence dissimilarity since previous work has found this cut-off to correspond closely to morphological species limits, as determined using both adults and larvae, for Trichoptera species of the Churchill region [21–23]. For completing the species presence/absence matrix, we performed additional COI sequencing to obtain species-level identifications for those taxonomic groups where our prior morphospecies designations did not perfectly match the COI genetic clusters.

After completing the COI work, we selected one specimen per species for the nuclear gene sequencing; we randomly selected an individual among those having a high-quality COI sequence (658 base pairs [bp], 0 ambiguous bases) that was closely related to the other specimens of that same species (i.e. not an outlier for its species on the neighbour-joining phenogram). Sequences of the nuclear genes were verified by using BLAST on GenBank and by building single-gene phylogenies to screen for contaminants, in addition to inspecting the amino acid alignment. All of our sequences are currently stored on BOLD in the public project EBTCH, Trichoptera Larvae of Churchill 2010 and on GenBank (COI: JX681817—JX682383; CAD: KR030383—KR030415; EF1-α: KR030340—KR030382; POL-II: KR271613-KR271655) (see Table D in S1 Table).

As Lepidoptera is the well-supported sister group to Trichoptera [37], we rooted all trees using sequences for Lepidoptera species downloaded from GenBank (GenBank IDs: Coleophora serratella COI—GU828594.1, CAD—GU828096.1, and EF1-α - GU828929.1; Phalera bucephala COI- GU828607.1, CAD—GU828108.1, and EF1-α - GU828941.1; Hypenodes humidalis COI—GU828672.1). No lepidopteran sequences for POL-II were available at the time of analysis for these species.

Phylogeny construction

To assess the influence of phylogenetic reconstruction on community structure metrics, we built in total six species-level phylogenies based upon a single individual per species: (1) multi-gene Bayesian; (2) COI Bayesian; (3) COI Bayesian phylogeny with a constrained backbone topology at the family level based on Holzenthal et al. [19]; (4) COI Bayesian with no backbone but with family-level branch lengths stretched 2x their original length and (5) 5x their original length, to test whether a simplistic correction for transitional saturation may result in COI-based trees yielding similar results as obtained using multi-gene nuclear phylogenies [16]; and finally (6) COI Neighbour Joining. Finally, as a negative control to examine metric behaviour, we also built a random tree by using the COI Bayesian phylogeny for the topology and branch lengths and then randomly shuffling the taxa among the tips in Mesquite 2.75 [38].

For the construction of the Bayesian phylogenies, we selected the best model of nucleotide substitution based upon the lowest Akaike information criterion (AIC) score for each gene using MrModelTest 2.3 [39], in conjunction with PAUP 4.0 [40]. For all genes, MrModelTest found the best model to be the generalized time reversible (GTR) with a gamma distribution parameter describing among-site rate variation and a proportion of invariant sites parameter. We built all Bayesian phylogenies in MrBayes 3.2 [41] using 4 chains with 10,000,000 generations, a sampling and diagnostic frequency of 1000, and a 25% burnin.

For the construction of the COI Neighbor Joining (NJ) tree, we used MEGA 5.0 to select the best model of nucleotide substitution and build the tree. A Tamura-Nei model with a gamma distribution was applied, and a bootstrap test with 1000 replicates was performed to assess the support for the tree.

Congruence of distance matrices

Since all phylogenies are converted to species pairwise distance matrices for calculating the phylogenetic community metrics, we first evaluated how different the input matrices were from one another by calculating the congruence among distance matrices (CADM) metric in R using the package ape [42–44]. We performed a posteriori testing to assess which matrices are the most incongruent.

Phylogenetic community structure metrics

For each of the 6 phylogenies, we calculated NRI and NTI. We calculated each metric using the same observed species presence/absence matrix created for local communities determined for Trichoptera larvae for the Churchill area from the DNA-barcode-validated morphospecies (see Table E in S1 Table). Each timepoint (June, July, August) for which two or more species were detected within a site was treated as a separate community, for a total sample size of 101 communities analyzed (see Table E in S1 Table). As this study is investigating impact of molecular dataset and phylogenetic analysis methods upon community phylogenetic metrics, no temporal analysis was performed here (but see [25]). In the Picante package in R, we calculated NRI and NTI using a null model with 1000 randomizations and 1000 iterations using the independent swap algorithm [45].

We used our multi-gene Bayesian phylogeny as the default against which to compare the community structure metrics generated using the other phylogenies. Although our goal is to assess the impact of input phylogeny upon conclusions from a real field study, not to assess which phylogeny is most likely to represent the truth, we assumed for our purposes here that the multi-gene Bayesian tree would be the strongest phylogenetic hypothesis, and therefore the most accurate for estimating metrics of phylogenetic community structure. Increased phylogenetic accuracy is usually achieved by increasing the number of independent genes and by utilizing more complex models of evolution [46–48]. We therefore performed a linear regression of the multi-gene Bayesian-calculated NRI and NTI values against those from all other phylogenetic hypotheses in sequence, with the origin forced through zero, as in Swenson [10], using R. A slope close to 1 would indicate that the values estimated by the other phylogenies are very similar to those estimated by the multi-gene Bayesian phylogeny, while a high r² would suggest that there is little variation between the values calculated by the other phylogenies vs. the multi-gene Bayesian phylogeny.

Phylogenetic signal metrics for traits

To test the influence of phylogeny on phylogenetic signal metrics for traits, we calculated two metrics prevalent in the literature, Blomberg et al.’s K statistic [17] and Pagel’s λ [18], for maximum body length and maximum case length (see Table F in S1 Table). Body size of Trichoptera larvae has been linked to habitat preferences [49,50] and is a commonly used ecological trait to test for phylogenetic signal in a variety of taxonomic groups [51]. Using a Nikon AZ100M microscope and NIS elements BR 3.0, we measured body length (mm) from the tip of the mandible to the end of the anal claw on specimens, while case length (mm) was the maximum distance spanned. After an initial visual sort of the largest specimens of each COI-identified species, microscopic measurements were taken on one to eight specimens, with the sample size dependent on abundance, enabling the measurement of maximum size for each species among the specimens we processed. Four species of Trichoptera (Hydropsyche alternans, Neureclipsis crepuscularis, Polycentropus aureoles, and Rhyacophila angelita) were excluded from the case length analysis, as they do not build portable cases, but either make fixed retreats or are free living [28].

We utilized the maximum length measured for each species and calculated Blomberg et al.’s [17] K statistic with 1000 replicates using the package Picante in R [45]. In addition, we calculated Pagel’s λ [18] using Geiger in R [52] and used a likelihood ratio test and a chi-squared distribution to compare the estimated λ to a tree with no phylogenetic signal (λ = 0), thus testing for significant phylogenetic structure in the trait data.

Distance matrix similarity

The distance matrices from all of the phylogenies (exempting the random phylogeny) were congruent. A posteriori testing showed that the mean of the Mantel correlations, computed on rank-transformed distances, between the multi-gene Bayesian distance matrix and all others were high (>0.9, p = 0.010 for all trees; Table 1). The trees with the highest to lowest congruence score, as compared with the multi-gene Bayesian tree, were: COI Bayesian + backbone > COI Bayesian > COI Bayesian 2x > COI Bayesian 5x > COI NJ.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Mantel correlations and probabilities of genetic distance matrices of phylogenies in comparison to the multi-gene Bayesian phylogeny.

https://doi.org/10.1371/journal.pone.0126662.t001

Estimation of phylogenetic community structure metrics

NRI and NTI values calculated using the multi-gene Bayesian phylogeny distance matrix were well estimated using the other phylogenies, excluding the random phylogeny (Table 2). We observed a strong association between the NRI and NTI values of the multi-gene phylogeny and NRI and NTI values for all other phylogenetic hypotheses (r²>0.75), and the slopes were slightly below 1.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Linear regression of phylogenetic community metrics using the multi-gene Bayesian phylogeny against those based on other phylogenetic hypotheses.

https://doi.org/10.1371/journal.pone.0126662.t002

Both the NRI and the NTI values were well estimated by the COI Bayesian + backbone phylogeny. The COI Bayesian + backbone phylogeny yielded NRI and NTI values that displayed the highest r² when regressed against the metrics as obtained using the multi-gene Bayesian phylogeny (0.898 and 0.878, respectively); slopes were close to 1 (0.945 and 0.893, respectively). The COI Bayesian 5x performed best in terms of the closest slope to 1 (NRI = 0.987, NTI = 0.903); however, it had a lower r² value than several other phylogenies (NRI = 0.848, NTI = 0.787).

Presence of phylogenetic signal

For body length, K showed significant phylogenetic signal for all of the phylogenies except the multi-gene Bayesian phylogeny (Table 3). Case length was shown to display phylogenetic signal based on K for only the NJ tree (p = 0.018). By contrast, λ detected significant phylogenetic signal across all of the phylogenies for both body length and case length (p<0.05). For both traits, the lowest λ values were calculated from the multi-gene Bayesian phylogeny (0.316 for body length, 0.308 for case length), while λ was somewhat higher for the other phylogenies and highest using the NJ tree (0.532 for body length, 0.581 for case length).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Phylogenetic signal metrics for Trichoptera maximum body length and maximum case length, using Blomberg et al.’s K metric [17] and Pagel’s λ [18].

https://doi.org/10.1371/journal.pone.0126662.t003

Best approaches for estimating phylogenetic community structure

Our study investigated how different phylogenetic reconstruction methods applied to COI data can approximate NRI and NTI calculated from a multi-gene tree, presumed here to provide the more robust phylogenetic hypothesis and branch lengths across the entire depth of the phylogeny. We found that both NRI and NTI values calculated from COI Bayesian phylogenies were generally concordant with those generated using the multi-gene tree, therefore not supporting our original hypothesis. Since the COI distance matrices were congruent with the multi-gene genetic distance matrix, and the phylogeny was generally well supported at the deeper nodes, it is not surprising that NRI and NTI values were also well estimated. Therefore, the rather surprising finding from our study was how well the COI data estimated the relative genetic distances between pairs of co-occurring species. All COI Bayesian phylogenies and the COI NJ tree had slopes slightly less than 1 for both NRI and NTI, indicating a slight bias towards overestimating the probable true value, i.e. increased detection of non-random phylogenetic community structure (Type I error). It is also interesting to note that the estimated NRI and NTI have a higher variance for values indicative of overdispersed phylogenetic community structure than for clustered values (S1 Fig and S2 Fig). This variance may be the result of the branch lengths being more poorly estimated for deeper nodes of the phylogeny for COI.

In terms of which phylogenetic reconstruction method most accurately estimated the multi-gene phylogenetic community structure metrics, we found that a Bayesian approach using COI while enforcing a backbone for the deeper relationships in the phylogeny performed well. The COI Bayesian + backbone phylogeny had the highest r² values for both NRI and NTI when comparing against the multi-gene Bayesian tree, indicating that this approach yielded values that varied little from the multi-gene Bayesian NRI and NTI values. In addition, the COI Bayesian + backbone phylogeny for NRI and NTI had slopes close to 1, suggesting that the values estimated using the COI Bayesian + backbone were very similar to the multi-gene Bayesian phylogeny. It is intuitive that this approach produced the most accurate NRI and NTI values, as it also produced the most congruent distance matrix with the multi-gene phylogeny. In addition, it was expected that a stronger phylogenetic reconstruction method including more biological information and a more realistic model of molecular evolution would produce a more accurate phylogeny and thus estimates of phylogenetic community structure metrics. As such, we would recommend that future community phylogenetic studies use realistic phylogeny reconstruction methods, as well as include information about the supported relationships among taxa from prior studies, particularly when using animal DNA barcode data for community phylogenetics.

Discrepancy of metrics measuring phylogenetic signal

Our study revealed substantial variation between the two metrics measuring phylogenetic signal in trait data but limited variability in conclusions across the phylogenetic reconstruction method employed. Similar to Münkemüller et al. [53], we found that Blomberg et al.’s K and Pagel’s λ suggested different conclusions about the phylogenetic signal of traits in our phylogenies. For instance, K suggested that both traits are not significantly conserved in the multi-gene Bayesian phylogeny; however, λ implied that they are significantly conserved. For all of the COI Bayesian phylogenies for case length, K found no support for phylogenetic signal, while λ found significant phylogenetic signal.

Münkemüller et al. [53] used simulated data to investigate the sensitivity of these metrics to phylogenetic structure and found that λ had similar values with repeated simulations, was less sensitive to variation in the number of species in the phylogeny, and was less prone to missing branch length information than K. In addition, λ had a smaller Type I statistical error rate and was able to detect phylogenetic signal in traits evolving under Brownian motion better than K, which was prone to Type II error [53]. Freckleton et al. [51] also reported appropriate Type I error rates for λ in simulated data and found that λ was able to detect phylogenetic signal in >90% of simulations for trees with 20 species and close to 100% for trees with 40 species. In addition, Freckleton et al. [51] found that 88% of published phylogenies display phylogenetic signal for at least one morphological or ecological trait with λ, and overall, 60% of traits displayed phylogenetic signal.

In our study, the K metric is failing to detect phylogenetic signal in the traits we measured, but λ suggests that phylogenetic dependence is present. Our phylogenies contain a reasonable sample size of species for this type of study (n > 40), as indicated by simulations and a review of λ values from empirical studies of morphological and ecological traits; therefore, our study should have strong power to detect significant phylogenetic signal, if present, using λ [51]. Since we do not explicitly know the process of trait evolution for our Trichoptera phylogenies, we cannot state which metric is superior to the other. However, we tend to favour Pagel’s λ for assessing phylogenetic dependence of trait data due to prior evidence of this metric’s behaviour [53], as well as the frequent finding of phylogenetic signal in body size measures across diverse taxa [51]. Focusing on our λ results, conclusions about the significance of phylogenetic signal did not vary across input trees; however, all of the COI phylogenies slightly overestimate phylogenetic signal when compared with the multi-gene tree. In sum, we conclude that choice of input tree was of modest impact in our study, when using λ, but that the Bayesian COI trees recovered λ values closer to those for the multi-gene phylogeny than did COI NJ trees.

Utility of COI

Our study has highlighted the utility of COI for estimating phylogenetic community structure for local communities within a small regional source pool within the order Trichoptera. Of course, when more genetic information is available, it is expected to be useful for better estimating NRI and NTI; nonetheless, our results suggest that using only COI is a reasonable approach at the taxonomic and geographic scale examined. We stress that for calculations of phylogenetic community structure, branching order is not as important as relative branch lengths within a dataset, as the tree is converted to a distance matrix for calculations. Short internal branches may be associated with rapid successive diversification events; different phylogenetic reconstructions may yield different branching orders yet similar relative values for the total sums of the branch lengths connecting pairs of tips. Therefore, the main consideration for community phylogenetics is being able to accurately estimate the relative genetic distances between species, which COI appears to do quite well at the geographic and taxonomic scale examined. However, for phylogenies reconstructed using only COI, theoretically we expected NTI to be more accurately estimated than NRI, since NTI focuses on the tips of the phylogeny where there is expected to be more support for a COI phylogeny. By contrast, our results indicated that NRI had higher r² values and slopes closer to 1 (in comparison with NRI from the the multi-gene tree) than NTI values. However, this finding may be an artefact of the communities in Churchill being principally phylogenetically clustered; the performance of NRI may be reduced when the communities are predominantly overdispersed, and this is a subject requiring further investigation. The Churchill area that we sampled is classified as a region according to Webb et al. (10–1,000 km) [2], and while our study suggests that COI is suitable at this scale for estimating NRI and NTI, difficulties may arise when examining larger geographic areas and broader taxonomic groups. It might be expected that a wider geographic area would include more taxa (species-area relationship; [54]), which is expected to increase phylogenetic accuracy. However, at a broader geographic scale, more taxonomic families may be sampled, yielding more deep branches requiring resolution, which could reduce the accuracy of a COI phylogeny and thus the NRI and NTI of the community. Nevertheless, improved phylogeny reconstruction methods are allowing researchers to build accurate phylogenies amongst insect orders with mitochondrial genomes [55]; with improved analytical methods, together with the availability of more backbone phylogenies, these problems may be alleviated. In addition, the accuracy of specimen identifications using COI may be reduced at broad vs. local geographic scales in some taxa [56]; hence, this should be another consideration for future studies. Despite these challenges, the large quantities of geo-referenced DNA barcodes being generated from initiatives such as the International Barcode of Life (iBOL) project (www.ibol.org) could be a vast, largely untapped resource of community data at varying geographic scales.