Mitochondrial genomes of Columbicola feather lice are highly fragmented, indicating repeated evolution of minicircle-type genomes in parasitic lice

Whole genome sequence data

We obtained whole genome shotgun Illumina sequencing reads from the NCBI SRA database for Columbicola macrourae ex. Zenaida macroura (SRR3161953), C. passerinae ex. Columbina picui (SRR3161931; hereafter “C. passerinae 1” following [43]), and C. passerinae ex. Columbina cruziana (SRR3161930; hereafter “C. passerinae 2”). We also generated new sequence data for a specimen of C. columbae from the feral rock pigeon (Columba livia). These species represent three of the four major clades identified within Columbicola by (Boyd et al., 2017) and thus are a good representation of this speciose genus (90 described species; Bush, Price & Clayton, 2009). Columbicola passerinae is a known complex of cryptic species parasitizing New World ground-doves (Johnson et al., 2007). Although they are currently described as a single species based on morphology, individuals from C. passerinae form two well-defined clades structured according to mitochondrial (cox1) and genomic-level nuclear sequence data, biogeography, and host association (Sweet & Johnson, 2016; Sweet & Johnson, 2018). This complex is therefore a useful system for comparing mitochondrial architecture between two very closely related species.

The Columbicola columbae specimen was photographed as a voucher (available on the Illinois Data Bank) and total genomic DNA was extracted using a Qiagen QIAamp DNA Micro Kit (Qiagen, Valencia, CA, USA) with a 48-hour incubation period at 55 °C in proteinase K. The resulting DNA was quantified with a Qubit 3.0 Fluorometer (Invitrogen, Carlsbad, CA, USA) using the standard protocol, fragmented with sonication for a mean insert size of 400 bp, and prepared for sequencing with a Hyper Library Preparation Kit (Kapa Biosystems, Wilmington, MA, USA). 150 bp paired-end reads were then generated in a single lane on an Illumina NovaSeq 6000 using a NovaSeq S4 reagent kit.

For both SRA and newly generated sequence data, we trimmed adapters and filtered reads following Boyd et al. (2017). We removed duplicate read pairs with fastqSplitDups.py (https://github.com/McIntyre-Lab/mcscript), removed the first 5 nt from the 5′-end of de-duplicated reads to account for potential low-quality bases using Fastx_clipper v.0.0.14 (FASTX Toolkit, http://hannonlab.cshl.edu/fastx_toolkit/), removed bases with a phred quality score <28 from the 3′-end using a 1 nt window in Fastx_trimmer v.0.014, and after trimming removed reads <75 nt. We confirmed adapter and quality trimming by running the trimmed libraries through FastQC v.0.11.7 (Babraham Bioinformatics, https://www.bioinformatics.babraham.ac.uk/projects/fastqc/).

Mitochondrial genome assembly of Columbicola

We used three methods to assemble the mitochondrial genomes of Columbicola: two approaches that combine read-mapping and de novo assembly, and a standard de novo approach. The former two methods are also complementary for assembling more challenging sequences (e.g., repetitive non-coding regions). We initially used these approaches to assemble the mitochondrial genome of C. passerinae 2 and then searched for similar patterns in the other three Columbicola taxa. Our methodology is further outlined in Fig. 1.

First, we used a combination of aTRAM v.2.0 (Allen et al., 2018) and MITObim v.1.8 (Hahn, Bachmann, & Chevreux, 2013) to assemble contigs of mitochondrial sequences. aTRAM uses BLAST (Altschul et al., 1990) to identify reads from a user-specified target and then uses third-party software to assemble those particular reads de novo. We attempted to assemble all 13 mitochondrial protein coding genes in C. passerinae 2 with aTRAM, using the amino acid sequences of each gene from the published Campanulotes compar mitochondrial genome as target sequences (AY968672; Covacin et al., 2006). For each assembly, we ran aTRAM for three iterations, used ABySS (Simpson et al., 2009) for de novo assembly, and used 10% of the read library. Following Sweet & Johnson (2018), we used a fraction of the library because mitochondrial reads are often much more prevalent than nuclear reads in a shotgun sequencing library, and using only a portion of these reads helps to avoid erroneous assemblies associated with numts, highly-repeated motifs, or contaminants (including contaminants from index- swapping (Carlsen et al., 2012; Schnell, Bohmann & Gilbert, 2015)). We then used the best- scoring contig from each gene assembly as a starting reference in MITObim. MITObim is a pipeline that uses MIRA v.4.0.2 (Chevreux & Suhai, 1999) to iteratively map reads against a starting reference (e.g., a single gene) and then against subsequent contigs to potentially reconstruct whole mitochondrial genomes. For each starting reference, we ran MITObim for a maximum of 100 iterations with the –quick option and utilizing paired-end information. Our MITObim runs resulted in a single contig for each starting reference. We then used Geneious v.11.1.5 (Biomatters Ltd., Auckland, NZ) to assemble the MITObim contigs de novo, using “Medium Sensitivity” and the default parameters.

Our second assembly method, an extension from our first approach, leverages the conserved control regions (CR) among potential minicircles. Because the CRs have highly conserved sequence blocks among minicircles in mammal lice, we expected there to be many more reads that map to the non-coding CR in a given contig compared to those which map to genic regions. However, identifying large read-coverage differences across a contig can be used for more than just evidence for presence/absence of minicircles. We can use this coverage information for a single contig to identify the boundary between the non-coding control and genic regions, and then map reads to this “boundary region” to identify reads from multiple minicircles. Song et al. (2019) used a similar approach to assemble fragmented mitochondrial chromosomes in mammal lice. To do this, we identified potential CRs by aligning the MITObim contigs produced from the first assembly approach with MAFFT (Katoh et al., 2002) in Geneious (default parameters) and identifying conserved regions among the contigs (Fig. S1). We identified coding regions using the beta version of MITOS2 to annotate the MITObim contigs (Bernt et al., 2013). Using each annotated contig as a reference, we then mapped the paired-end reads used by MITObim for the assembly of that particular contig. Read-mapping was carried out in Geneious using “Medium-Low Sensitivity” and up to 5 iterations. As predicted, we identified coverage spikes upstream and downstream from the regions with annotated genes (Table 1, Fig. S2). For the cox1 contig, we extracted a 250 bp region immediately upstream from the coding region, and mapped the entire 10% paired-end read library (i.e., not only reads that mapped to the contig in MITObim) to this region using Geneious (same parameters as above). We then identified reads that mapped to and extended beyond the 3′end of the reference. Theoretically, these reads include sequences from all minicircles. We assembled the reads de novo in Geneious and identified contigs with mean coverage >20X. We then aligned these contigs with MAFFT, trimmed the conserved 5′ends of each contig, and used the trimmed contigs (which are therefore putative coding portions of the minicircles) as starting references in MITObim (max 100 iterations, 10% read library).

Third, we used a standard de novo approach to assemble the C. passerinae 2 mitochondrial genome. We used SPAdes v.3.11.1 (Bankevich et al., 2012) to assemble the 10% read library using the –plasmid option to invoke plasmidSPAdes (Antipov et al., 2016). We assessed the assembly quality with the web version of QUAST (Gurevich et al., 2013). We then converted the resulting scaffolds to a BLAST-formatted database in Geneious, and identified scaffolds containing mitochondrial protein-coding genes using tblastn with amino acid sequences from Campanulotes compar as queries.

Assemblies with aTRAM, MITObim, and plasmidSPADES were run on a cluster with 24 cores, two 2.60 GHz Sky Lake processors, and 384 GB of memory maintained by Information Technology at Purdue (West Lafayette, IN, USA). De novo assembly, read mapping, and BLAST analyses in Geneious were run locally on a 2.9 GHz Intel i5 processor with 8 GB of memory.

To annotate the C. passerinae 2 mitochondrial genome, we initially used MITOS2 to annotate the assembled contigs. We then manually fine-tuned gene annotations (following the criteria in Cameron, 2014b) and identified the number of minicircles by focusing on annotated regions with conserved flanking regions. We identified conserved regions by using BLAST to find highly similar sequences (>80% identical, e-value <1e⁻²⁰) and searching for these sequences both within and between different contigs. We tested for evidence of circularization by identifying identical stretches of sequence >50 bp near the 3′  and 5′ ends of our contigs. We also confirmed circularization using AWA (Machado et al., 2018) with kmer length=31 and a minimum circular contig length of 1,000 bp (Table S1). After finalizing the annotations, we used Polymerase Chain Reactions (PCR) to verify our assembled minicircles (Article S1, Fig. S3, Table S2).

We also searched for minicircles in three other representatives of the Columbicola genus: C. columbae, C. macrourae, and C. passerinae 1. Evidence for minicircles in these taxa would indicate that the architecture is conserved across the genus. We used our first two approaches to assemble the mitochondrial genomes of these taxa: first generating longer contigs with aTRAM and MITObim and then assembling reads associated with CR/coding region boundaries to produce starting references for MITObim. We tested for a consistent architecture within Columbicola by using BLAST to compare our annotated C. passerinae 2 contigs to contigs from the other three species. We then compared the variation in gene-order and content between species and verified MITOS annotations following the process outlined in Cameron (2014b).

Mitochondrial genome assembly of Pediculus humanus

As another confirmation of our computational approaches, we assembled the mitochondrial genome of Pediculus humanus, a species with known minicircles verified with both PCR and Southern Blotting (Shao, Kirkness & Barker, 2009). We downloaded the P. humanus A 2013 paired-end Illumina data from the GenBank SRA database (SRR5088471), cleaned the library, and assembled the minicircles using the same three approaches used for Columbicola (MITObim from aTRAM-assembled starting seeds, MITObim from CR/PCG edge starting seeds, and de novo with SPAdes). We used published gene sequences (translated PCGs and rRNA) from P. humanus (Shao, Kirkness & Barker, 2009) as targets for aTRAM, and mapped reads from the MITObim cox1 assembly to identify a CR/protein-coding boundary. Because the Pediculus library was considerably larger than the Columbicola libraries (128,102,260 cleaned reads), we used 1% of the library (compared to 10% of the Columbicola libraries) in each assembly approach.

Genome data

Paired-end library sizes varied among the four Columbicola genomes. There was an average of 7,598,374 reads from 10% of the adapter and quality-trimmed libraries. The number of reads ranged from 3,557,924 in C. passerinae 1 to 7,617,990 in C. macrourae (Table 2). We report statistics for 10% of the libraries because these were the reads used for assembling mitochondrial genomes.

Mitochondrial genome assembly of Columbicola passerinae 2

The software aTRAM produced assembled contigs for seven protein-coding genes in C. passerinae 2 using the Campanulotes compar amino acid bait sequences: cox1 (1,530 bp), cox2 (873 bp), cox3 (801 bp), cob (1,089 bp), nad1 (783 bp), nad3 (477 bp), and nad5 (945 bp). The reported contig lengths are from the best-scoring contigs that we subsequently used as starting references in MITObim. Resulting MITObim contigs had the following lengths: cox1: 3,246 bp (5 iterations); cox2: 2,665 bp (7 iterations); cox3: 3,597 bp (10 iterations); cob: 3,673 bp (11 iterations); nad1: 2,973 bp (7 iterations); nad3: 2,465 bp (7 iterations); nad5: 3,913 bp (10 iterations). An average of 61,459 reads mapped to each contig in the final MITObim iteration (Table S3). When we mapped these reads to the MITObim contigs in Geneious, all contigs showed elevated coverage near the 5′  and 3′  ends of the contigs and lower coverage toward the centers of the contigs. However, coverages toward the centers were still high. For example, coverage at the ends (positions 1–805 and 2,378–3,246) of the contig produced with the cox1 reference was 12 times higher than coverage toward the center (positions 806–2,377) of the contig (Fig. S2, Table S3). The other contigs showed a similar pattern; among all contigs the coverage at the ends of each contig were on average 11 times higher than coverage near the centers of the contigs (Table 1, Table S3). When we annotated the contigs with MITOS, in every instance the lower-coverage regions corresponded to coding genes, whereas the regions with high-coverage peaks did not have any genes annotated with high confidence. We were unable to produce any additional contigs by assembling all the MITObim contigs together in Geneious. However, aligning each contig to itself with BLAST indicated >50 bp regions of identical sequences in the high-coverage 5′  and 3′  regions, thus suggesting each contig is a complete circle (Fig. S2).

Although our first approach indicated there are several complete minicircles in Columbicola, we were not able to assemble every mitochondrial gene with this approach. One goal of our second approach was to assemble a set of contigs that contained all of the expected mitochondrial genes. The alignment of MITObim contigs from our first approach showed similar sequences near the ends of all contigs but highly variable sequences towards the middle of the contigs (Fig. S1). The variable regions also corresponded to annotated genes and with regions of lower coverage. Together, these results suggest the conserved, high-coverage regions are part of the conserved CRs in the minicircles. We identified 6,400 reads from the 10% read pool (6,226,688 reads) that mapped to the boundary of a CR and the 5′  end of cox1. De novo assembly of these reads resulted in 19 contigs with ≥20X coverage. Alignments of the 19 contigs showed conserved sequences near the 5′  ends and variable sequences towards the 3′  ends (Fig. S4). After trimming the conserved regions, the 19 contigs had an average length of 118 bp (min 31, max 139). Using the trimmed 19 contigs as starting references, MITObim produced contigs with an average length of 2,983 bp (min 1,814; max 4,213) from an average of 10 iterations per contig. BLAST identified >50 bp regions of identical sequences on the 5′  and 3′  ends of each MITObim contig.

The plasmidSPADES de novo assembly of the 10% read pool produced 2,024 contigs (N50: 12,832) and 2,020 scaffolds. We found BLAST hits for every Campanulotes compar protein-coding and ribosomal RNA gene among the scaffolds. However, atp8, nad4l, and nad6 BLAST returns had high e-values >0.01 and low bit scores <25, indicating the homology of these sequences should be interpreted with caution. The longest scaffold with a high-scoring BLAST hit was 1,733 bp. These scaffolds were highly similar to corresponding contigs (i.e., containing the same protein-coding gene) assembled with MITObim (from both aTRAM and control edge starting references), with the only differences coming near or outside the gene boundaries identified by MITOS (Fig. S5). However, none of the de novo-assembled scaffolds covered the entire length of a MITObim contig. Relevant data for all three assembly approaches, including starting references and resulting contigs/scaffolds, are available on the Illinois Data Bank.

Annotation of the Columbicola passerinae 2 mitochondrial genome

An initial automated annotation of the control edge MITObim contigs with MITOS and subsequent manual annotations identified all 13 protein-coding genes, both rRNA genes, and 19 of the 22 tRNAs. These genes are separated on 17 separate mini-chromosomes. Although we recovered 19 contigs from MITObim, the assembly process can produce multiple contigs which possess identical coding regions, but whose CRs have sequence differences in portions of the CR outside of those conserved across all minicircles of C. passerinae 2. These ‘duplicate’ contigs were treated as assembly artefacts and removed from subsequent analyses. We determined the number of distinct chromosomes using a combination of the results described above, identifying conserved stretches of sequence and annotating the intervening regions (Fig. S6). Each protein-coding and rRNA gene is on its own circular chromosome, with two chromosomes containing only tRNA genes. The chromosomes range from 1271 bp to 2862 bp with two contigs (cox1 and nad4) forming incomplete circles, as their assemblies lacked non-coding sequence in the 3′  direction of the coding gene (Fig. 2, Tables S3–S4). AWA and PCR results were consistent with the lengths of our assemblies and further confirmed circularity of most chromosomes (all but the chromosomes containing cox1 and nad4) (Table S1, Fig. S3).

Assembly of other Columbicola species

Assemblies of the three other Columbicola species showed similar patterns to C. passerinae 2. Contigs from MITObim using the aTRAM-assembled references did not assemble together in Geneious, but when aligned showed conserved regions of sequence with highly variable intervening regions. Although there was variation in the number of mapped reads and coverage values, for each species the conserved non-coding CRs in each contig had higher (4 to 9 times) coverage than the variable coding regions (Table 1). MITOS only identified genes in the variable, lower coverage regions. De novo assembly of reads that mapped to control/coding boundaries resulted in a similar number of contigs for most species (Table 2). MITOS was able to recover most genes from the resulting extended contigs from MITObim. BLAST searches identified similar regions flanking the annotated genes. As with C. passerinae 2, we used these flanking regions to identify individual minicircles. We recovered evidence for 15-17 minicircles in each species, with lengths of recovered complete chromosomes ranging from 1207-3118 bp for C. passerinae 1, 1262-3173 bp for C. macourae, and 1119-2657 bp for C. columbae (Fig. 2, Table S4).

Assembly of Pediculus humanus

Assemblies of the P. humanus mitochondrial chromosomes showed similar patterns to the Columbicola assemblies. aTRAM assembled 11/13 PCGs (all excluding nad3 and nad4l). Using these assemblies as starting references in MITObim produced contigs with an average length of 3,014 bp. Read-mapping to these MITObim contigs showed large relative coverage difference between coding and non-coding regions (average of 12 times higher coverage in non-coding regions; Table 1 and Table S3). Using the contigs assembled from the CR/cox1 edge reads, we recovered 17/18 chromosomes (we were unable to assemble the chromosome that contains rrnS) that showed high sequence similarity (>99% in coding regions) to the published P. humanus mitochondrial sequences from Shao, Kirkness & Barker (2009) (Fig. S7 ; contigs and alignments of published sequences are available on the Illinois Data Bank). There was lower similarity between the non-coding CRs, but most of these contigs still had >90% total similarity to the published sequences. Four contigs had particularly lower similarity (70–80%) with published sequences in the CRs, and five contigs were shorter than their comparable published chromosomes. We were able to assemble and identify 11 PCGs and both rRNA genes from de novo assembly. However, as with Columbicola, these contigs were shorter than contigs produced from the other two assembly approaches; most assemblies seemed to fail at the CRs (Fig. S8). There was also a chimeric contig containing portions of nad5 and rrnL.

Comparative mitochondrial genome structure among Columbicola

Overall mitochondrial genome structure was consistent across Columbicola species, with each protein-coding and rRNA gene on a separate chromosome and two minicircles with only tRNAs. All four genomes showed AT compositional bias (average 67.2 AT% across all coding regions of the minicircles, ranging from 63.4–69.2 total AT% among the four species; Table 2, Table S4), which is consistent with some insect mitochondrial genomes and comparable to other louse mitochondrial genomes (Cameron, Johnson & Whiting, 2007; Cameron, 2014a). Additionally, protein coding genes aligned to the corresponding genes from the C. passerinae 2 assembly and had expected uncorrected distances between overlapping regions (13% between C. passerinae 1, and 22%–28% between the more distantly related taxa). None of the completely assembled genes were significantly longer or shorter than homologous copies in other Columbicola species or other louse species, which suggests the assemblies were reliable. However, nad5 was incompletely assembled in C. passerinae 1 and C. columbicola. Assembled contigs from C. macrourae did not include minicircles with either atp8 or nad6; however, this is likely an assembly or annotation failure rather than evolutionary loss of these genes within this species. Both genes are comparatively short and have higher rates of amino-acid substitution than other louse mitochondrial genes creating challenges for BLAST. Additional aTRAM and BLAST searches using the C. passerinae 2 amino acid sequences of atp8 and nad6 as references likewise failed to identify these two genes in C. macrourae, further highlighting the difficulty of assembling or annotating them from short-read genomic data. Similarly, each Columbicola species assembly was missing between two (C. columbae) and five (C. macrourae) tRNA isotypes and had one or more duplicated tRNA isotypes. None of the canonical 22 tRNA isotypes were missing across all four Columbicola species, suggesting that missing tRNAs are likely present on other tRNA-only minicircles which failed to assemble or be adequately annotated.

Although the overall genome structure was consistent across the four Columbicola mitochondrial genomes, there was substantial variation in gene arrangement between Columbicola species. Of the 17 minicircles identified, only one (trnI-cox1) had a common gene arrangement across all four Columbicola species. Each of the other 14 minicircles containing a PCG or rRNA gene differed between at least two species due to the presence, absence, or translocation within the minicircle of accompanying tRNAs. These differences varied among comparisons: C. passerinae 2 and C. passerinae 1 shared 10 minicircles, C. passerinae 2 and C. columbae shared 7,  C. passerinae 2 and C. macrourae shared 6, C. passerinae 1 and C. columbae shared 7, C. passerinae 1 and C. macrourae shared 4, and C. macrourae and C. columbae shared 9. The assembled and annotated chromosomes for all four Columbicola species are available on GenBank (Table 2). Relevant data generated from the various analyses, including assembled contigs and read-mapping references, are available on the Illinois Data Bank.

Assembling fragmented mitochondrial genomes with a fully bioinformatic approach

Assembling fragmented mitochondrial genomes can be a challenging and labor-intensive process. In this study, we were able to assemble the mitochondrial minicircles of lice using a fully bioinformatic approach from raw genomic Illumina reads. In particular, we used an approach that assembled short, variable contigs from reads that mapped to the edge of the control region (CR) and coding regions. These contigs served as starting references for assembling longer contigs, essentially acting as “bioinformatic PCR primers.” By identifying conserved non-coding CRs that flank adjacent coding regions, we could then detect separate circular chromosomes. Subsequent annotations recovered most of the 37 mitochondrial genes present in metazoans.

We used two other methods for assembling mitochondrial genomes: de novo and read-mapping extension. De novo assembly is fast (∼12 min on five cores for SPAdes) and has the advantage of not requiring a starting reference. However, finding mitochondrial scaffolds is reliant on BLAST searches with a closely related target. In our case, we were not able to reliably identify every PCG or rRNA gene in each Columbicola species, presumably because the initial target was too distantly related (Campanulotes), the missing genes were particularly divergent (e.g., nad6), or the genes simply failed to assemble. It is also challenging to identify tRNA-only minicircles with a de novo approach, because tRNA genes are short (<100 bp) and can be highly divergent from even a closely-related target. Our de novo assemblies also had trouble assembling through the CR because of the high similarity of this region among minicircles, and as a result the de novo scaffolds are considerably shorter than contigs produced from other approaches (Fig. S5). Assemblies with aTRAM were also fast (assembled PCGs in ∼3.5 min total) and using the assembled genes as starting seeds with read-mapping in MITObim allowed us to extend the contigs into the CR. However, we were limited by the starting seeds we could assemble with aTRAM, which as with standard de novo assembly tended to exclude faster-evolving PCGs and tRNA genes. Our approach based on the coding/non-coding boundary combines the benefits of de novo and read-mapping approaches, while also enabling us to assemble all mitochondrial minicircles (perhaps excluding small chromosomes that only contain tRNAs).

Our approach is an important progression that builds on previous work that primarily used wet laboratory methods or a combination of bioinformatics and laboratory methods to identify and confirm mitochondrial minicircles. The first minicircular mitochondrial genome from insects, that of Pediculus humanus, was assembled in 18 separate contigs from Sanger shotgun reads, with non-coding regions sequenced from clone libraries of PCR products generated by designing primers within the coding regions of assembled contigs (Shao, Kirkness & Barker, 2009). Circularization was also confirmed by restriction digest and Southern hybridization of purified mitochondrial DNA extractions from Pediculus (Shao, Kirkness & Barker, 2009). We were able to recover these same chromosomes from Pediculus using our bioinformatic approaches. Subsequently, a procedure was developed for sequencing minicircles by first sequencing several partial coding gene sequences, designing outward facing primers which would amplify the non-coding portion of several minicircles, then designing primers in conserved regions of the non-coding region to amplify all minicircles present in a given species. This mixed population of amplicons would then be sequenced using NGS methods to yield clean, minicircle contigs (Shao et al., 2012). Recently, Song et al. (2019) simplified this approach to minimize primer design and PCR steps by using initial partial coding sequences (generated by universal PCR primers) as seeds to directly assemble minicircles from genomic extracts. Similar to our approach, conserved CRs were identified for select mitochondrial genes (initially sequenced with Sanger) from whole genome assemblies and subsequently used to identify all contigs derived from mitochondrial minicircles. Finally, Song et al. (2019) confirmed putative minicircle findings with PCR amplification and gel electrophoresis.

Approaches that rely on laboratory methods to either identify or confirm fragmented mitochondrial genomes are certainly useful and should continue as an option for these types of studies, particularly because there are some limitations to assembling these genomes with Next Generation Sequencing (NGS) data. For example, because it is repeated across multiple minicircles, it is challenging to assemble through long stretches of relatively conserved non-coding regions (e.g., CRs) with short-read data. The entire CR may not assemble for a particular minicircle, which prohibits confirmation of circularization. This was the case for a few of our assembled fragments (incomplete circles in Fig. 2). Due to their similarity, reads from the CRs could also misassemble among different minicircles, so the exact sequences of the CR should be interpreted with caution or reported with ambiguity codes (as we have done in this study). Additionally, a bioinformatic approach may fail to recover all genes present in the organism. This could particularly be an issue for tRNA-only minicircles or highly divergent PCGs such as atp8 and nad6. However, laboratory-based methods can also fail to identify minicircles or genes, especially tRNAs. Additionally, the most rigorous wet-lab methods for confirming the presence of minicircles, such as Southern blotting, require large volumes of starting tissue, thus limiting their utility to species that are either physically large (unlike lice) or are maintained in culture, which is rarely feasible for the vast majority of louse species.

It should also be noted that bioinformatics methods do not appear to erroneously assemble fragmented mitochondrial genomes for taxa where they are not present. Many hundreds of insect mitochondrial genomes have been assembled from NGS datasets (e.g., Linard et al., 2018) using bioinformatics assembly procedures such as those presented here (e.g., de novo or MITObim baited assembly). Complete circularized mitochondrial genomes are routinely recovered, or when incomplete, assemblies typically only lack portions of the CR (e.g., Gómez-Rodriguez et al., 2015). In these assemblies there is no evidence for chimeras between incompletely assembled portions of genomes. Within insects, evidence for mitochondrial genome fragmentation derived from NGS sequencing and bioinformatics analysis is confined to taxa for which prior evidence for fragmentation had been collected (such as lice). Conversely, even within taxa such as lice where aberrant mitochondrial genomes appear to be common, NGS-only approaches have successfully assembled complete, canonical, single-chromosome mitochondrial genomes when such an architecture was present in a given taxon (e.g., for Falcolipeurus (Song et al., 2019).

Thus despite potential limitations, with an ever-increasing reliance on NGS data it is important to demonstrate that fragmented genomes can be assembled from sequence read libraries. In this regard there is capacity for the study of “aberrant” mitochondrial genome architectures to follow the path of “conventional” mitochondrial genome studies over the past decade where the expansion of NGS techniques has led to a sharply reduced reliance on PCR-centric methods (see review Cameron (2014a)). Avoiding PCR is especially relevant for rare specimens or for incorporating publicly available data (NCBI SRA), where it may not be possible to confirm genomic architecture with PCR due to a lack of DNA extracts. We assembled mitochondrial genome fragments from whole genome shotgun sequencing, but our approach could also be used with reads generated from target capture methods (e.g., ultra-conserved elements (UCEs)). These genome reduction methods target nuclear regions, but mitochondrial reads are also captured at relatively high depth as “by-catch” in the process (Amaral et al., 2015; Ströher et al., 2017). Utilizing the ever-growing diversity of genomic data, most of which is publicly available through databases such as GenBank, allows for further exploration of mitochondrial genome architecture and arrangements on a broad taxonomic scale. The bioinformatic methods described here accurately reconstruct aberrant mitochondrial genomes whose structures have been verified by wet-lab procedures (i.e., Pediculus) and extend our capacity to search for fragmented mitochondrial genomes into taxa which cannot be verified by wet-lab methods (i.e., Columbicola and other feather-lice).

Evolutionary implications for highly fragmented mitochondrial genomes in Columbicola

Highly fragmented mitochondrial genomes in lice were first identified in human lice (Pediculus humanus) (Shao, Kirkness & Barker, 2009), and minicircle-type mitochondrial genomes have since been confirmed in several other genera of sucking lice (Anoplura (Jiang, Barker & Shao, 2013; Dong et al., 2014a; Dong et al., 2014b; Herd, Barker & Shao, 2015; Shao et al., 2017), elephant and warthog lice (Rhyncophthirina Shao et al., 2015), and chewing lice in the family Trichodectidae (Cameron et al., 2011; Song et al., 2019). Phylogenetic analyses have indicated that these three lineages of mammalian lice form a monophyletic group, which suggests highly fragmented mitochondrial genomes evolved once in the common ancestor of this group, as suggested by Song et al. (2019). Given the rarity of this architecture in other groups of organisms and the implications of fragmentation for disease, it is very reasonable to conclude there have been few transitions from full to fragmented mitochondrial genomes. However, our study implies a second such transition to minicircle-type mitochondrial genomes in parasitic lice. The mitochondrial genomes of Columbicola species are the first recorded instance of a fragmented mitochondrial genome from within the family Philopteridae for which the complete genome has been sequenced. Other published complete mitochondrial genomes of Philopteridae include species of Bothriometopus (Cameron, Johnson & Whiting, 2007), Ibidoecus (Cameron et al., 2011), Falcolipeurus (Song et al., 2019), Campanulotes (Covacin et al., 2006), and Coloceras (Cameron et al., 2011). All five of these species have a single circular mitochondrial chromosome, although Coloceras has a full-sized chromosome and at least one additional heteroplasmic fragment missing approximately one third of the full-sized genome. Recent phylogenomic analysis shows that Columbicola is nested within these other three taxa (Johnson et al., 2018), so the minicircularized genomes of Columbicola likely represent an independent transition from a single chromosome (Fig. 3). It is possible that fragmentation occurred multiple times within the Columbicola genus, which could explain why we recovered fewer minicircles in C. columbae. However, given the similar levels of fragmentation and conserved gene order in some minicircles among the four Columbicola species, fragmentation is at most an infrequent phenomenon within the genus and most parsimoniously represents a single independent origin. A more densely-sampled analysis is needed to confirm the history of fragmentation within the genus.

Cameron et al. (2011) identified partial fragmented mitochondrial genomes in three other genera of Philopteridae (Anaticola, Quadraceps, and Philopterus), so it is possible other transitions to fragmented genomes from single chromosomes have occurred in these feather lice. However, none of these taxa possess the minicircle-type genome architecture (i.e., few genes per chromosome and a relatively long CR with conserved motifs) reported here for Columbicola. The independent origin of mitochondrial minicircles in Columbicola and mammal-associated lice allows for phylogenetically independent examination of the life-history traits proposed to be related to this phenomenon. As Columbicola are feather-feeding and associated with birds, the proposals that mitochondrial minicircularization in lice is related to blood- or mammal feeding are unlikely (c.f. Dong et al., 2014b). Similarly, the generation-time (24 days) of Columbicola, which is the most fragmented mitochondrial genome yet sequenced, falls within the range of generation times (Martin, 1934) for other lice whose mitochondrial genomes are complete, or at least less fragmented, suggesting that genome fragmentation is not causally related to short generation times. Multiple independent transitions from single to fragmented mitochondrial chromosomes in one group of organisms (parasitic lice) do, however, suggest that there are underlying potentially shared mechanisms for the formation of minicircles, such as the lack of a functional mtSSB gene, and/or similar selective pressures retained across relatively deep evolutionary time (∼65 my). Future work can address these questions by building on the genome architectural framework established in this study.