Computational protocol: Mitochondrial Cox1 Sequence Data Reliably Uncover Patterns of Insect Diversity But Suffer from High Lineage-Idiosyncratic Error Rates

Similar protocols

Protocol publication

[…] DNA was extracted using the DNeasy Animal Tissue kit (Qiagen, Hilden, Germany). We sequenced the 3′ end of cox1 using the primers Jerry (F: CAA CAT TTA TTT TGA TTT TTT GG) and Pat (R: TCC AAT GCA CTA ATC TGC CAT ATT A) . Although this is not the fragment proposed as a standard barcode (5), our results can be generalized because the average evolutionary rate of both cox1 fragments is similar. Roe & Sperling sequenced and evaluated the information from the entire cox1–cox2 region, showing that “ultimately, no single optimally informative [for barcoding] 600 bp location was found within the 2.3 kb of COI–COII, and the DNA barcoding region was no better than other regions downstream in COI”. Our results may also, to some extent, apply to any mitochondrial protein-coding gene with a similar evolutionary rate, as problematic issues with mitochondrial DNA related to species delineation are linked to the mitochondrial genome per se rather than individual genes , . Sequences were edited in Sequencher 4.8 (Genecodes Corp., Ann Arbor, MI, USA) and translated into amino acid sequences for alignment control and screening for internal stop codons or other anomalies in MacClade . Finally, nucleotides were aligned using MUSCLE under default settings on CIPRES Portal v.2 (www.phylo.org). New sequences have been submitted to GenBank. [...] Parsimony searches were run in the program TNT version 1.1, which we also used to run 500 jackknife (character removal 36%) replications to assess node stability (hit best tree five times command, keeping 10,000 in memory). We ran maximum likelihood (ML) analyses with the program GARLI on CIPRES Portal v.2. We used the GTR+I+G model as selected by MrModeltest for the combined dataset and ran analyses until 10,000 generations revealed no significant improvement of the likelihood scores of the topology. Bootstrap values were based on 250 replicates using only the datasets of the different Australian radiations due to computational limitations. We also ran ML analyses in RAxML v7.0.4 , bipartitioning the data (1st +2nd versus 3rd codon sites) and implementing a GTR model with CAT approximation to incorporate rate heterogeneity across sites. Haplotype networks based on statistical parsimony were calculated using TCS 1.13 (; 95% connection limit). This approach subdivides the variation based on the level of homoplasy within the data themselves, i.e., distinguishes between long (homoplastic) and short (non-homoplastic) branches, which provides a relative measure of the divergence within a given dataset, rather than using a priori determined thresholds. Independent haplotype networks generally agree with named species or species groups , . [...] We ran a neighbor-joining analysis using uncorrected p-distances for fast distance-based clustering of the data. The SpeciesIdentifier module of TaxonDNA software v.1.6.2 was used to study the genetic divergences in our dataset and to cluster sequences at different preset thresholds using uncorrected p-distances (; http://code.google.com/p/taxondna/). SpeciesIdentifier accounts for threshold violations according to triangle inequity (i.e., when the divergence between A – B and B – C is 3% or less, but A – C exceeds 3%, then A, B and C would still be grouped into one 3% cluster by Taxon DNA). SpeciesIdentifier recognizes a priori delineated species from the sequence name, as long as the name follows the format “Genus species”, i.e., “Rhantus suturalis”, or “Rhantus australiaone MB1307”. The output summarizes the number of different species names in the dataset, the number of clusters found under the present threshold (e.g., 1%, 2%), the number of clusters that contain only one species name, and the number of perfect clusters (those that contain all individuals under one species name and only those individuals, i.e., monophyly). Thus, we can calculate the number of split clusters (one species split into more than one cluster, i.e., paraphyly) and lumped clusters (more than one species name in a cluster). SpeciesIdentifier was used for species richness estimation, with clusters taken as species surrogates. For any clustering threshold (e.g., at 1%, 2%, 3%…), two values were reported. The first of these values was the number of clusters found relative to the number of morphology-based species names in the dataset (agreement hereafter). For example, a dataset with a hundred species names and a threshold clustering at 25% divergence would likely reveal only one cluster. Thus, our species richness estimation would amount to a meager 1% (agreement) of the actually present species as delineated by morphology. Second, and more importantly, we report taxonomic accuracy, which was calculated as the number of perfect clusters (i.e., clusters containing all sequences of a morphology based species and only those sequences) relative to the number of species in the dataset. The number of perfect clusters can increase when the existing taxonomy is revised to accommodate cryptic or overlooked species. A one hundred percent accuracy means that all clusters perfectly mirror the species hypotheses based on morphology.Character-based group delineation, or population aggregation analysis (PAA) , was used to delineate geographically endemic subgroups or species within groups a priori identified by clustering and phylogenetic analyses. The sequences of the species were manually screened for diagnostic characters in the DNA sequence editor Se-Al (http://tree.bio.ed.ac.uk/software/seal/). Specifically, we applied PAA to several supposedly recent morphospecies that were lumped into paraphyletic species clusters. In a pairwise step addition, PAA defines populations based on the presence of at least one diagnostic (fixed) character in one population; otherwise, they are merged and then compared to another population. For species delineation, taxonomists traditionally use diagnostic characters from morphology or behavior, and this usage can be extended to nucleotide characters , .We also implemented the method developed by Pons et al. and explained in detail in a previous publication that delineates genetic clusters using a generalized mixed Yule coalescent (GMYC) model that represents independently evolving entities. This method uses a maximum likelihood approach to optimize the shift in the branching patterns of the gene tree from interspecific branches (Yule model) to intraspecific branches (neutral coalescent). The model optimizes the maximum likelihood value of a threshold, such that the nodes before the threshold are identified as species diversification events, while the branches beyond the threshold are clusters following coalescent processes. This method has previously been implemented in other taxonomically understudied groups , , , , . In large trees, a unique species-populations split for a particular time (single threshold) may not reflect the true diversification for all of the lineages included; therefore, we performed an analysis allowing multiple and independent thresholds over time and across the tree . Generalized mixed Yule coalescent clustering was performed here using the R package SPLITS (SPecies' LImits by Threshold Statistics), which allows single or multiple thresholds (http://r-forge.r-project.org/projects/splits/). This package provides confidence intervals (CI) in the output (solutions within two log-likelihood units of the maximum likelihood), but the GMYC entity content of these solutions cannot be retrieved in the output at present. Before running SPLITS, an ultrametric tree was made fully dichotomous, and branches with zero branch length were pruned or removed using the package ape in R . The underlying tree was derived from our above RAxML v7.0.4 analysis. Identical sequences were removed from the analysis using the reduced dataset provided by RAxML. ML tree searches were run 100 times starting from different parsimony trees, and the best one tree was finally optimized. Branch lengths were made ultrametric using PATHd8 software by arbitrarily setting the root node to 100 Ma. This age was chosen because it approximately renders the so called standard rate of nucleotide substitution of 2.3% per Ma in insects . The standard rate, which was suggested for the species and genus levels based on several arthropod examples, roughly agrees with the rates reported for different groups of Coleoptera based on calibrations using different biogeographical events , , and for the whole Coleoptera using fossils .We report the number of GMYC entities and the number of perfect GMYC entities containing all and only the members of an a priori identified species. […]

Pipeline specifications