Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Conflicting Evolutionary Histories of the Mitochondrial and Nuclear Genomes in New World Myotis Bats

Conflicting Evolutionary Histories of the Mitochondrial and Nuclear Genomes in New World Myotis Bats Downloaded from https://academic.oup.com/sysbio/article/67/2/236/4095671 by DeepDyve user on 13 July 2022 Syst. Biol. 67(2):236–249, 2018 © The Author(s) 2017. Published by Oxford University Press, on behalf of the Society of Systematic Biologists. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. For Permissions, please email: journals.permissions@oup.com DOI:10.1093/sysbio/syx070 Advance Access publication August 26, 2017 Conflicting Evolutionary Histories of the Mitochondrial and Nuclear Genomes in New World Myotis Bats 1 2 1 3 3 ROY N. PLATT II ,BRANT C. FAIRCLOTH ,KEVIN A. M. SULLIVAN ,TROY J. KIERAN ,TRAVIS C. GLENN ,MICHAEL W. 1 4 1 5 1,∗ VANDEWEGE ,THOMAS E. LEE Jr. ,ROBERT J. BAKER ,RICHARD D. STEVENS , AND DAVID A. RAY 1 2 Department of Biological Sciences, Texas Tech University, 2901 Main St, Lubbock, TX, USA; Department of Biological Sciences and Museum of Natural Science, Louisiana State University, 202 Life Science Building, Baton Rouge, LA, USA; Department of Environmental Health Science, University of Georgia, 206 Environmental Health Sciences Building, Athens, GA, USA; Department of Biology, Abilene Christian University, 1600 Campus Ct. Abilene, TX, USA; and Natural Resource Management, Texas Tech University, 2901 Main St, Lubbock, TX, USA Correspondence to be sent to: Department of Biological Sciences, Texas Tech University, Lubbock, TX, USA; E-mail: david.4.ray@gmail.com. Received 17 March 2017; reviews returned 31 July 2017; accepted 15 August 2017 Associate Editor: Matthew Hahn Abstract.—The rapid diversification of Myotis bats into more than 100 species is one of the most extensive mammalian radiations available for study. Efforts to understand relationships within Myotis have primarily utilized mitochondrial markers and trees inferred from nuclear markers lacked resolution. Our current understanding of relationships within Myotis is therefore biased towards a set of phylogenetic markers that may not reflect the history of the nuclear genome. To resolve this, we sequenced the full mitochondrial genomes of 37 representative Myotis, primarily from the New World, in conjunction with targeted sequencing of 3648 ultraconserved elements (UCEs). We inferred the phylogeny and explored the effects of concatenation and summary phylogenetic methods, as well as combinations of markers based on informativeness or levels of missing data, on our results. Of the 294 phylogenies generated from the nuclear UCE data, all are significantly different from phylogenies inferred using mitochondrial genomes. Even within the nuclear data, quartet frequencies indicate that around half of all UCE loci conflict with the estimated species tree. Several factors can drive such conflict, including incomplete lineage sorting, introgressive hybridization, or even phylogenetic error. Despite the degree of discordance between nuclear UCE loci and the mitochondrial genome and among UCE loci themselves, the most common nuclear topology is recovered in one quarter of all analyses with strong nodal support. Based on these results, we re-examine the evolutionary history of Myotis to better understand the phenomena driving their unique nuclear, mitochondrial, and biogeographic histories. [Summary tree methods; concatenation; vespertilionidae; phylogenomics; UCE; ultraconserved elements.] The bat genus Myotis (Order Chiroptera, Family Ves- single-genetic unit that is susceptible to evolutionary pertilionidae) comprises more than 100 species that ori- processes that may cause its history to diverge from ginated during the last 10–15 million years (Stadelmann the history of the species (Edwards and Bensch 2009). et al. 2007), making it one of the most successful extant The most widely accepted phylogenies of Myotis rely mammalian radiations. Myotis are distributed world- heavily on mitochondrial data and even phylogenies wide, excluding polar regions, and generally share a containing nuclear data demonstrate an over reliance on similar ecological niche: aerial insectivory. Myotis species mitochondrial markers for resolution, likely swamping often exhibit little morphological differentiation and, as out potentially conflicting signals from the nuclear a result, the rate of cryptic speciation within the genus is genome. For example alignments of the nuclear RAG2 thought to be high. For example, specimens identified and mitochondrial cytb contained 162 and 560 variable as M. nigricans and M. albescens form multiple para- characters, respectively (Stadelmann et al. 2007). Despite phyletic lineages distributed throughout the phylogeny these concerns, we find ourselves in a situation where of Neotropical Myotis (Larsen et al. 2012). Confounding our understanding of evolutionary relationships, our matters, the morphological variation that exists is often basis for conservation strategies, and our biogeographic a poor indicator of species-level relationships. Early clas- hypotheses are all founded, almost entirely, on mito- sifications of Myotis identified three major morphotypes chondrial phylogenies. (Findley 1972). Subsequent phylogenetic analyses of the Targeted sequencing of ultraconserved elements mitochondrial cytochrome-b (cytb) gene demonstrated (UCEs; Faircloth et al. 2012) to collect sequence data paraphyletic origins of each morphotype, suggesting from thousands of loci across the nuclear genome has frequent convergent evolution (Ruedi and Mayer 2001). resolved a number of difficult phylogenetic problems These same analyses demonstrated that geography was (e.g. see Faircloth et al. 2012; McCormack et al. 2013; a better predictor of phylogenetic relationship than mor- Green et al. 2014; Faircloth et al. 2015; McGee et al. 2016). phology (Ruedi and Mayer 2001; Stadelmann et al. 2007). Here we used UCE baits to collect ~1.4 Mbp from >3600 The ability of mitochondrial markers to resolve a nuclear loci in addition to random shotgun sequencing to well-supported topology does not guarantee that the collect full mitochondrial genomes in 37 taxa, primarily mitochondrial tree represents the species tree (e.g. see representing New World (NW) Myotis. Rather than Willis et al. 2014; Li et al. 2016; Leavitt et al. 2017). The analyzing each sequence data set once, or a handful of lack of recombination and uniparental inheritance of times, we chose to use a range of sampling, partitioning, the mitochondrion means that it is transmitted as a and inference methods to fully explore phylogenetic tree [08:41 13/7/2018 Sysbio-OP-SYSB170072.tex] Page: 236 236–250 Downloaded from https://academic.oup.com/sysbio/article/67/2/236/4095671 by DeepDyve user on 13 July 2022 2018 PLATT ET AL.—PHYLOGENOMICS OF MYOTIS BATS 237 space of the nuclear and mitochondrial genomes. We recovered six times, all by analysis of the translated recovered 294 trees representing 175 distinct topologies protein coding genes (Fig. 1a). Despite being the most for the nuclear UCE data and 28 trees representing common topology, bipartition frequencies and clade 14 distinct mitochondrial topologies. Our results show probability values were lower than nucleotide based data that, despite the range of trees recovered from each sets, as would be expected with more slowly evolving marker type, nuclear and mitochondrial markers occupy amino acid sequences. Analysis of RNA coding genes distinct regions in tree space. Given that the nuclear and recovered two topologies depending on the method of mitochondrial trees are distinct from one another it is phylogenetic inference. The remainder of the analyses necessary to reevaluate hypotheses made based solely recovered various trees without any pattern relating to on the mitochondrial phylogeny. partitioning scheme or inference method. A majority rule consensus tree for all 28 analyses is shown in Figure 1c. The mitochondrial consensus tree reflects pre- viously reported mitochondrial phylogenies for Myotis. RESULTS The New World Myotis are generally monophyletic We used targeted sequencing of UCEs in 37 individu- splitting into Nearctic and Neotropical clades with the als (Table 1) to collect sequence data from 3648 nuclear Old World taxon, M. brandtii, sister to the Neotropical loci which we assembled into concatenated alignments clade. Ambiguity, represented by polytomies, is found as large as 1.37 Mb. In addition, we assembled mito- only in terminal relationships. chondrial genomes for most taxa from random shotgun sequencing data. We then used the data to infer the phylogenetic history of New World Myotis using a range UCE Assembly and Alignment, and Phylogenetic Inference of locus sampling strategies, alignment partitioning We averaged 3.29 million reads per sample after methods, and phylogenetic inference methods. demultiplexing reads from the UCE-enriched, sequen- cing pool. These reads were assembled into an aver- age of 5778 contigs per sample (min=1562 M. nyctor, Mitochondrial Assembly, Alignment, and max=11,784 M. nigricans ). Recovery rates for UCE loci Phylogenetic Inference varied across taxa. Of the 5500 loci in the Amniote We generated an average of 5,498,440 paired reads per probe set, we successfully recovered 3898 UCE loci, sample (min = 2,2370,68, max = 9,961,348) of random 3648 loci from five or more samples and 212 loci in all DNA sequence data for mitochondrial genome assembly. 37 samples (Table 2). On average, 3332 UCE loci were Assemblies varied in quality. Some were almost entirely recovered per sample, ranging from 1106 (M. nyctor) complete while others were missing small regions. We to 4008 (M. keaysi). Repetitive sequences, identified via were not able to assemble the entire mitochondrial RepeatMasker searches, were minimal, occupying less control region for any of the samples. We found three than 0.02% of sites across all UCE alignments. premature stop codons in mtDNA protein coding genes. In all, 294 different combinations of alignment, par- Subsequent manual alignment and validation suggested titioning, and inference method were used to identify that these regions were miscalled by MitoBim, and we 175 unique nuclear topologies. The most common tree corrected the errors prior to analysis. Sequence coverage was recovered 45 times (Fig. 1b). All nodes were highly- of the mitochondrial genomes averaged 58× (range supported based on clade probability values and all >1×–297×). but three nodes were highly supported with maximum The 37 mitochondrial protein coding, rRNA and tRNA likelihood bootstrap replicates. A consensus tree (Fig. 1d) genes were concatenated into an alignment of 15,520 of the 294 nuclear topologies resolves New World Myotis bp containing 5007 informative characters. rRNA and (and M. brandtii) as a single polytomy that excludes tRNA genes were concatenated into an alignment con- M. volans and represents a departure from previous taining 4157 bp containing 862 parsimony informative mitochondrial phylogenies, including the one recovered characters. Protein coding genes were concatenated into herein. Comparison of the mitochondrial and nuclear an alignment containing 11,363 bp or 3770 amino acids consensus trees (Fig. 1c,d) reveals substantial differences (aa) and containing 4145 or 509 parsimony informative in the relationships within the genus. characters, respectively. For the alignment containing As noted above, we were concerned that phylogenetic all protein coding, rRNA and tRNA genes, 30 samples error, sampling bias or other methodological factors were  90% complete, and alignments for five samples may be driving the conflict between the mitochondrial were 68–84% complete. Only 21% and 50% of nucleotide and nuclear trees. To address this, we investigated positions were present in the M. albescens (TK 61766) the effects of matrix composition (or completeness) and M. levis alignments, respectively. on phylogenetic inference by generating 10 alignments When considering alignment type, partitioning having varying levels of matrix completeness. Nine strategy, and phylogenetic inference method we ana- alignments were composed of loci with 15 to 95% (at 10% lyzed the mitochondrial data 28 ways as described intervals) of taxa. The tenth alignment was composed in the Methods section. These analyses recovered 14 of all loci which were recovered in 100% of specimens distinct topologies, the most common of which was examined. For example, all loci that were represented [08:41 13/7/2018 Sysbio-OP-SYSB170072.tex] Page: 237 236–250 Downloaded from https://academic.oup.com/sysbio/article/67/2/236/4095671 by DeepDyve user on 13 July 2022 238 SYSTEMATIC BIOLOGY VOL. 67 TABLE 1. Specimens examined Species Name herein Museum identification Number of Short read Mitochondrial genome or accession number UCE loci archive accession accession Eptesicus fuscus fuscus GCA_000308155.1 2467 NA KF111725.1 Eptesicus fuscus fuscus TK 178736 2849 SAMN06141755 MF143474 Myotis albescens albescens RDS 7889 2764 SAMN06141756 MF143471 Myotis albescens albescens QCAZ 9157 1185 SAMN06141757 MF143497 Myotis albescens albescens TK 61766 2872 SAMN06141758 MF143470 Myotis ater ater M4430 2774 SAMN06141760 MF143493 Myotis auriculus auriculus MSB 40883 2229 SAMN06141761 MF143486 Myotis brandtii brandtii GCA_000412655.1 2446 NA KT210199.1 Myotis californicus californicus UMMZ 175828 2948 SAMN06141763 MF143469 Myotis davidii davidii GCA_000327345.1 2450 NA KM233172.1 Myotis diminutus diminutus QCAZ 9168 3078 SAMN06141777 MF143481 Myotis dominicensis dominicensis TK 15624 2576 SAMN06141764 MF143467 Myotis evotis evotis MSB 47323 2586 SAMN06141765 MF143468 Myotis fortidens fortidens MSB 54941 2791 SAMN06141766 MF143483 Myotis horsfieldii horsfieldii MHNG 1926.039 3017 SAMN06141767 MF143494 Myotis keaysi keaysi TK 13525 3195 SAMN06141768 MF143477 Myotis keenii keenii UAM 113849 2723 SAMN06141769 MF143472 Myotis leibii leibii TK 24872 3119 SAMN06141770 MF143488 Myotis levis levis RDS 7781 2538 SAMN06141771 MF143482 Myotis lucifugus lucifugus GCA_000147115.1 2429 NA KP273591.1 Myotis lucifugus lucifugus MSB 46679 2736 SAMN06141772 MF143491 Myotis melanorhinus melanorhinus TK 193888 3177 SAMN06141774 MF143489 Myotis nigricans nigricans QCAZ 9601 2854 SAMN06141775 MF143495 Myotis nigricans nigricans TK 194145 3159 SAMN06141776 MF143484 Myotis nyctor nyctor TK 151413 856 SAMN06141773 MF143498 Myotis occultus occultus MSB 121995 2957 SAMN06141778 MF143490 Myotis oxyotus oxyotus UMMZ RCO1013 3106 SAMN06141779 MF143479 Myotis riparius riparius TK 22688 2924 SAMN06141782 MF143473 Myotis riparius riparius TK 101723 2990 SAMN06141759 MF143475 Myotis riparius riparius TK 145199 2890 SAMN06141780 MF143480 Myotis ruber ruber MVZ 185692 2757 SAMN06141781 MF143478 Myotis septentrionalis septentrionalis TK 194420 2916 SAMN06141762 MF143487 Myotis thysanodes thysanodes 07LEP 2821 SAMN06141783 MF143492 Myotis velifer velifer MSB 70877 2704 SAMN06141784 MF143499 Myotis vivesi vivesi MSB 42658 2469 SAMN06141785 MF143476 Myotis volans volans MSB 40886 2819 SAMN06141786 MF143496 Myotis yumanensis yumanensis TK 194144 2589 SAMN06141787 MF143485 Samples beginning with “GCA_” represent genome assemblies available through NCBI. Species with more than one sample are designated with a superscript. Specimens derived from whole genome alignments are designated with a superscript “G.” MSB = Museum of Southwestern Biology; MVZ = Museum of Vertebrate Zoology; MHNG = Natural History Museum of Geneva; QCAZ = Pontificia Universidad Catolica del Ecuador Museo de Zoologia; TK = Texas Tech University Natural Science Research Laboratory; UAM = University of Alaska Museum of the North; UMMZ = University of Michigan Museum of Zoology. in at least 15% of species (n  5) were included generations. Unfortunately, computational limits forced taxa in the alignment (n = 3648). In the “100%” matrix, us to abandon the individually partitioned, Bayesian loci only loci which were identified in all species (n = analyses. In general, the same alignment produced taxa 37) were included (n = 212). Alignments were parti- the same topology regardless of inference method or loci tioned using three schemes: unpartitioned, individually partitioning scheme with the only exception being the partitioned by locus, or combined partitions. Align- terminal relationships of the M. levis/M. albescens clade ment lengths, numbers of informative characters and in the combined versus unpartitioned Bayesian analysis number of partitions identified by PartitionFinder are of the 100% complete data matrix. available in Table 2. Resulting alignments were analyzed We also tested how locus length could affect phylo- with Bayesian, maximum likelihood, and coalescent genetic discordance. To do so, trees were generated from methods using RaxML, ExaBayes, ASTRID, ASTRAL-II, data matrices incorporating UCE loci of differing lengths and SVDquartets as described in the Methods section. (Hosner et al. 2016). All 3648 loci were grouped into Bootstrap topologies stabilized in each alignment and 10 separate bins based on locus length. The number partition combination within 150 replicates and all of informative characters per bin ranged from 1115 to Bayesian runs converged in less than ten thousand 6995 and the number of informative characters was [08:41 13/7/2018 Sysbio-OP-SYSB170072.tex] Page: 238 236–250 Downloaded from https://academic.oup.com/sysbio/article/67/2/236/4095671 by DeepDyve user on 13 July 2022 2018 PLATT ET AL.—PHYLOGENOMICS OF MYOTIS BATS 239 FIGURE 1. Comparison of nuclear and mitochondrial phylogenetic trees of Myotis. The most common mitochondrial (a) and nuclear (b) topologies. The mitochondrial topology arises from Bayesian and maximum likelihood analysis of the amino acid sequences portioned based on the PartitonFinder recommendations. The nuclear topology is from the 55% complete data matrix partitioned based on the recommendations of PartitionFinder. For clarity, bootstrap values greater than 95 and clade probability values greater than 0.95 have been omitted from the nuclear consensus tree. Majority-rule consensus trees from 28 mitochondrial (c) and 294 nuclear (d) topologies. Values above the branches in the consensus trees (c,d) are bipartition frequencies for that clade across all nuclear or mitochondrial topologies. Conflicting tips between data types (consensus nuclear vs. consensus mitochondrial) are indicated with red lines between the topologies. Biogeographic regions are color coded, as are subgeneric classifications based on morphotypes, as defined by Findley (1972). Species with more than one sample are designated with a superscript that is referenced in Table 1. Specimens derived from whole genome alignments are designated with a superscript “G.” [08:41 13/7/2018 Sysbio-OP-SYSB170072.tex] Page: 239 236–250 Downloaded from https://academic.oup.com/sysbio/article/67/2/236/4095671 by DeepDyve user on 13 July 2022 240 SYSTEMATIC BIOLOGY VOL. 67 TABLE 2. Matrix composition analysis Matrix ID (% taxa Min. num. Number of loci Parsimony- Variable Alignment Optimum with data from locus) of taxa characters informative uninformative length partitions 15 5 3648 38,718 62,588 13,77,262 31 25 9 3379 37,894 57,672 12,60,248 37 35 12 3232 37,259 56,453 12,27,093 34 45 16 3064 36,539 54,782 11,87,492 31 55 20 2890 35,284 53,200 11,44,471 27 65 24 2668 34,373 51,148 10,91,620 27 75 27 2481 33,031 48,732 10,41,099 20 85 31 2034 29,179 41,711 903,903 16 95 35 1193 18,288 24,189 575,321 15 100 37 212 3480 4778 112,125 6 General alignment information. For a subset of analyses a series of alignments were generated based on the number of taxa represented by each locus. Thirty-seven taxa were examined so an alignment with all 37 taxa was considered 100% complete. Parsimony-informative characters make up a small portion of the total alignment. The optimum partitioning scheme was calculated with PartitionFinder. correlated with average locus length (Supplementary concatenation methods only recovered the same tree in Material Fig. 1 available on Dryad at http://dx.doi.org/ 1 of 100 attempts. 10.5061/dryad.5g205). On average, only 2.6% of char- Finally, we used weighted and unweighted statist- acters in each bin were parsimony-informative (Sup- ical binning to combine individual gene trees into plementary Material Fig. 1 available on Dryad at supergenes, estimate the supergene phylogeny, and then http://dx.doi.org/10.5061/dryad.5g205). Each of the infer the species tree from the supergene trees. The ten length-based alignments recovered slightly different 3648 loci were combined into 528 binned loci with 480 topologies. Terminal relationships were generally stable bins containing 7 loci each and 48 bins containing 6 across analyses with the majority of differences between loci each. Binning methods increases the normalized topologies found in the early bifurcations of Myotis. quartet scores from an average of 0.547 with gene trees Generally, longer alignments produced well resolved to 0.672 for the binned-unweighted and 0.673 for the and similar/identical topologies with significant nodal binned-weighted supergene trees. Given the relatively support regardless of the phylogenetic method or parti- even distribution of loci into bins the negligible differ- tioning scheme used. In contrast, smaller data sets were ence in quartet/species tree discordance between the more likely to yield unique topologies. unweighted and weighted supergene trees is expected. Given the goal of fully exploring nuclear tree space Binning methods have been shown to recover incorrect and generating as many reasonable nuclear UCE based species trees under coalescent methods by altering gene topologies for comparison with the mitochondrial tree, tree frequencies. In addition, it is possible that support we decided to randomly sample a limited number for incorrect trees can actually increase with the number of nuclear UCE loci to create small alignments that of input gene trees (Liu and Edwards 2015). With these would be more likely to result in unique topologies. We concerns in mind, both binning methods recovered the therefore randomly subsampled UCE loci to create 100 same topology which was the most common nuclear unique data sets of 365 loci. Loci were concatenated in topology observed across all analyses (Fig. 1c). each replicate data set and analyzed using maximum likelihood in RAxML. Of the 100 alignments analyzed, 80 unique topologies were generated (mean Robinson– Foulds distance = 4.3). Topological Comparisons In addition to concatenation methods we explored phylogenetic tree space of the UCE data set using All mitochondrial and nuclear trees are available species tree methods. Normalized quartet scores from in Supplementary Material File 1 available on Dryad ASTRAL-II (Mirarab and Warnow 2015) analyses were (http://dx.doi.org/10.5061/dryad.5g205). When visu- consistent among all analyses, with scores ranging from alizing all topologies in tree space, nuclear trees co- 0.540 to 0.553, and the number of induced quartet localized and were distinct from mitochondrial topolo- gene trees ranging from 7,745,739 (100% complete 212 gies (Fig. 2a). Comparisons of Robinson–Foulds sym- gene trees) and 63,042,410 (15% 3648 loci). SVDquartets metrical differences show 26 symmetrical differences (Chifman and Kubatko 2014) sampled all 66,045 quartets. between the most similar nuclear and mitochondrial On average, the total weight of incompatible quartets trees (Fig. 2b). For comparison only 9 of 43,070 and 0 of was 2.84%. Similar to the concatenated analysis, we 378 pairwise nuclear versus nuclear and mitochondrial inferred coalescent-based species from the same 100 versus mitochondrial tree comparisons, Figure 2c and subsamples of 365 loci described above. Despite being d, respectively, exhibited more than 26 symmetrical generated from the same underlying data, summary, and differences. [08:41 13/7/2018 Sysbio-OP-SYSB170072.tex] Page: 240 236–250 Downloaded from https://academic.oup.com/sysbio/article/67/2/236/4095671 by DeepDyve user on 13 July 2022 2018 PLATT ET AL.—PHYLOGENOMICS OF MYOTIS BATS 241 FIGURE 2. Comparison of mitochondrial and nuclear topologies (Fig. 1a). All trees recovered from mitochondrial and nuclear analyses were visualized in tree space using multidimensional scaling of Robinson–Foulds distances. Nuclear trees (green outlines) were clustered separately from mitochondrial trees (blue outlines). Robinson–Foulds distances between nuclear and mitochondrial trees (Fig. 1b), nuclear versus nuclear trees, (Fig. 1c) and mitochondrial versus mitochondrial trees (Fig. 1d) indicate the overall distinctions between the nuclear and mitochondrial topologies. Of the 45 analyses that recovered the most frequently parameters to produce a range of reasonable nuclear and observed topology (Fig. 1b), 38 were Bayesian and mitochondrial topologies that may be more useful than RAxML searches that varied by matrix completeness overreliance on a tree resulting from one or two analyses. and partitioning scheme. The fact that these analyses Rather than rejecting concordance between the two data recovered the same topology is expected given that they types from a single analysis we took steps to analyze are not independent. For example, a RAxML analysis both data sets to generate many viable topologies to of the 25% complete data set uses 1.26 Mb of the potentially reveal hidden overlap between marker types 1.38 Mb of data present in the 15% complete matrix, masked by phylogenetic error. sharing 91% of the data. Analyses that directly varied We recovered 175 and 14 unique nuclear and mito- the alignments and/or sampled less data (e.g. randomly chondrial topologies from 294 and 28 analyses using sampling loci) were more likely to generate unique multiple methodologies, sampling strategies, and phylo- topologies than the nested analyses described above. genetic parameters. By varying our phylogenetic meth- Of the 200 analyses that randomly sampled UCE loci, ods we show that recovering a consistent topology is 164 unique topologies were observed. This implies that difficult even when using the same marker type, thus when analyses of large data sets produce well-resolved making it difficult to determine the “best” topology. trees with significant nodal support, sampling smaller Rather than choosing single trees to represent the portions of the data, may provide a mechanism for nuclear and mitochondrial data, we use consensus trees creating phylogenetic uncertainty not represented by to represent the ambiguity that is present through typical tree scoring metrics. all analyses, in addition to phylogenies conforming to the most common topology (Fig. 1). Despite our efforts to identify as many plausible, distinct nuclear and mitochondrial trees as possible we were unable to DISCUSSION recover overlapping topologies between marker types We analyzed 3648 UCE loci and mitochondrial gen- suggesting that nuclear UCE loci and the mitochondrial omes of 37 Myotis bats and outgroups. The mitochon- genomes of Myotis have distinct evolutionary histories. drial and nuclear UCE phylogenies recovered distinct A comparison of topologies in tree space illustrates the topologies (Fig. 1), whether comparing the most com- differenced in trees generated for each data set (Fig. 2a). monly recovered or consensus topologies. Previous Intramarker discordance, as measured by differences in work with UCE loci demonstrated that support for nuclear versus nuclear or mitochondrial versus mito- deep divergences varied based on the number of loci chondrial topologies is low. Pairwise comparisons of examined (McCormack et al. 2013). Further, bootstrap nuclear topologies contain an average of 8.3 symmetrical replicates and clade probability values can be inaccurate differences, mitochondrial comparisons contain an aver- metrics of nodal support (Douady et al. 2003, Hedtke age of 10.8 differences. By contrast, pairwise differences et al. 2006). We varied the input data and phylogenetic between nuclear and mitochondrial topologies contain [08:41 13/7/2018 Sysbio-OP-SYSB170072.tex] Page: 241 236–250 Downloaded from https://academic.oup.com/sysbio/article/67/2/236/4095671 by DeepDyve user on 13 July 2022 242 SYSTEMATIC BIOLOGY VOL. 67 an average of 33.3 symmetrical differences. The two most agrees with Stadelmann et al. (2007)exceptinthe similar nuclear and mitochondrial topologies contain placement of M. septentrionalis and M. auriculus. Our 26 symmetrical differences. Of the 43,072 pairwise nuclear UCE topology places M. septentrionalis and M. comparisons of nuclear topologies, only 27 contain more auriculus sister to a clade containing M. californicus, M. than 26 symmetrical differences. Likewise of the 378 leibii, and M. melanorhinus. Even though this placement pairwise comparisons of mitochondrial topologies only of M. septentrionalis and M. auriculus disagrees with one contains more than 26 symmetrical differences. Stadelmann et al. (2007) it is the same relationship These observations indicate that discordance between recovered by Haynie et al. (2016). Thus, even though marker types (Fig. 2b) is much greater than within there were some concerns and miscalled nucleotides marker type (Fig. 2c,d). All but the most divergent in the MITObim assemblies, these inaccuracies have nuclear topologies are more similar to each other than not caused the mitochondrial tree to deviate from any nuclear versus mitochondrial comparison. the expectations derived from previous work. Given Conflict between mitochondrial and nuclear data may the consistency of our mitochondrial phylogeny with be driven by error in phylogenetic estimation or may previously published and independently recovered reflect genuine discordance between the two marker topologies, we are confident that the mitochondrial types (Degnan and Rosenberg 2009; Huang et al. 2010). phylogeny we recovered here, and by others, reflects We relied on multiple tree-inference methods (e.g. the true mitochondrial tree. However, the mitochondrial summary vs. concatenation), manipulated phylogenetic topology may not adequately reflect the species history, parameters (e.g. partitioning strategy), and sampling particularly when considering the factors that cause criteria (e.g. loci sampled in all taxa) to minimize the incongruence between nuclear and mitochondrial gene impacts of phylogenetic error on the data sets. In most trees. Causes of conflicting gene trees include horizontal cases, varying parameters or methodologies generated transfer, gene duplication, introgressive hybridization, unique topologies with minor differences. In most cases and incomplete lineage sorting. The rapid radiation of unique topologies resulted from rearrangements of a this clade as well as other biological factors suggests few terminal taxa. For example, placement of M. volans that some of these phenomena are more likely to have and M. brandtii was often either sister to the remaining influenced the Myotis radiation than others. New World Myotis or as an early bifurcation between Horizontal transfer of genes is thought to be rare in the Nearctic and Neotropical clades. Myotis vivesi was at eukaryotes, but, vespertilionids in general (Thomas et times found as sister to the clade containing M. lucifugus, al. 2011; Platt et al. 2014; Platt et al. 2016), and Myotis M. occultus, and M. fortidens or as sister to the clade in particular (Pritham and Feschotte 2007; Ray et al. containing the neotropical Myotis. 2007; Ray et al. 2008; Pagan et al. 2010), have experi- Summary methods failed to recover the most com- enced horizontal transfer of DNA transposons. These mon nuclear topology presented in Figure 1b except events would not be reflected in our phylogeny because when loci were binned prior to gene tree estimation. repetitive sequences were removed prior to phylogenetic Summary methods also tended to recover more unique analyses. More generally, gene duplications could create topologies than concatenation methods when analyzing conflicting signal among individual UCE markers (e.g. data from the same gene(s). For example, the random comparing nonorthologous UCE loci), but the number sample analyses recovered 80 unique topologies using of gene duplication events would have to be very high concatenation (RAxML) and 89 unique topologies with to impact enough of the 3648 UCE loci to confound the summary methods (ASTRAL-II). This likely has to do mitochondrial and nuclear phylogenies. Further ruling with the limited number of informative characters per out gene duplication events as the dominant cause of locus and by extension limited phylogenetic signal per conflicting phylogenetic signal is the fact that such events gene tree (Supplementary Material Fig. 1 available on are likely depressed in Myotis as evidenced by their Dryad at http://dx.doi.org/10.5061/dryad.5g205). In smaller genome size (~2.2 Gb) and trend towards DNA these instances, limited phylogenetic signal per gene loss (Kapusta et al. 2017) combined with low rates of would likely lead to increased opportunity for phylogen- paralogy in UCEs in general (Derti et al. 2006). etic error in gene tree estimation. Further supporting this Introgressive hybridization and reticulation could idea, binning of compatible UCE loci may have indirectly significantly influence the phylogenies of Myotis in a increased phylogenetic signal resulting in the same topo- way that leads to conflicting signal between the nuclear logy that many of the concatenation analyses recovered. and mitochondrial genomes (Sota 2002; Good et al. No other summary/coalescent method recovered this 2015). Hybridization in bats may be relatively common topology. given their propensity to swarm at cave entrances for Previous studies that used primarily mitochondrial breeding purposes. In European Myotis, swarming has data recovered effectively the same relationships among allowed for high degrees of hybridization between M. Myotis as our mitochondrial analyses (Ruedi and Mayer brandtii, M. mystacinus, and M. alcathoe (Bogdanowicz et 2001; Stadelmann et al. 2007; Roehrs et al. 2010; Larsen al. 2012). Further, M. evotis, M. thysanodes, and M. keeni all et al. 2012; Ruedi et al. 2013; Haynie et al. 2016). experienced historical gene flow during their divergence Differences between our tree and any single previ- (Carstens and Dewey 2010; Morales et al. 2016). We could ously published trees are minor. For example, our tree also explain the differences between the mitochondrial [08:41 13/7/2018 Sysbio-OP-SYSB170072.tex] Page: 242 236–250 Downloaded from https://academic.oup.com/sysbio/article/67/2/236/4095671 by DeepDyve user on 13 July 2022 2018 PLATT ET AL.—PHYLOGENOMICS OF MYOTIS BATS 243 and nuclear UCE phylogenies if Myotis experienced relationships link species with different morphotypes. extensive incomplete lineage sorting during their radi- Based on these results, it appears that the three major ation. Two factors influence the rate of lineage sorting, morphotypes in Myotis are indeed a result of convergent the fixation rate and the speciation rate (Hudson et al. evolution, as suggested by previous work (Ruedi and 2002). Increasing the time to fixation and/or decreasing Mayer 2001; Stadelmann et al. 2007) despite the conflict- the amount of time between cladogenic events will ing mitochondrial and nuclear phylogenies recovered. increase the likelihood of incomplete lineage sorting. Among the more dramatic differences between the Myotis are generally long-lived species (Dzeverin 2008) nuclear and mitochondrial topologies is the placement of and underwent a rapid radiation between 5 and 10 M. volans and M. brandtii as sister to all New World taxa MYA (Lack et al. 2010), suggesting that Myotis species by the nuclear data. Our mitochondrial analyses place are likely to experience higher levels of lineage sorting. M. volans within a Nearctic clade and M. brandtii directly The importance of these events—introgressive hybrid- inbetween the Nearctic and Neotropical bifurcations ization and incomplete lineage sorting—in driving the (Fig. 1a,c) as has been previously reported (Stadelmann differences between the mitochondrial and nuclear et al. 2007). In contrast the consensus UCE tree (Fig. 1d) phylogenies should be further explored. places M. brandtii, a species distributed throughout the Old World, within a polytomy at the base of the New World Myotis radiation and M. volans sister to all New Evolutionary History of Myotis World Myotis (including M. brandtii). The placement of M. brandtii in the nuclear UCE consensus tree does Our understanding of relationships within Myotis is not necessarily contradict the mitochondrial phylogeny heavily biased by mitochondrial data because nuclear since it is an unresolved polytomy. However, it is markers were harder to collect and produced fewer worth noting that the most commonly recovered nuclear informative sites (Ruedi and Mayer 2001; Stadelmann et al. 2007; Lack et al. 2010; Roehrs et al. 2010; Larsen topology (Fig. 1b) places M. brandtii sister to all New et al. 2012; Ruedi et al. 2013; Haynie et al. 2016)or World Myotis (excluding M. volans). This relationship nuclear phylogenies with large numbers of makers would more closely affiliate M. brandtii with other Old contained limited numbers of taxa (Platt et al. 2015). World taxa. The placement of M. volans as sister to all Our UCE-based results indicate that nuclear trees vary New World taxa (including M. brandtii) in the most substantially from the mitochondrial tree. Given that the common nuclear tree is a significant departure from nuclear and mitochondrial trees are different, we find previous work and, at first glance, does not make as it necessary to re-evaluate Myotis in the context of the much sense in a biogeographic framework. Myotis volans nuclear data. is distributed across western and northwestern North Paraphyly of M. nigricans and M. albescens was inferred America as far as far north as Alaska. Myotis brandtii is from previous mitochondrial phylogenies. Larsen et al. distributed across much of Northern Europe and into (2012) identified a minimum of 4 and potentially 12 the extreme eastern regions of Siberia. The key may lie lineages in M. albescens and M. nigricans. Our sampling in understanding the biogeography and phylogenetics included three M. albescens and two M. nigricans, of a third species, M. gracilis, a species that we were compared with Larsen’s 17 and 29 samples. Despite unfortunately unable to include. different mitochondrial and nuclear topologies overall, Myotis gracilis, along with M. brandtii,are theonlytwo our mitochondrial and nuclear phylogeny recovered the Myotis geographically distributed in the Old World, but same paraphyletic clade of three M. albescens samples phylogenetically affiliated with the New World Myotis and M. levis. Close relationships between these taxa (Stadelmann et al. 2007). If future work verifies the sister were found in previous work and are expected. More relationship between M. brandtii and M. gracilis, then we can envision a scenario where M. gracilis, M. brandtii, importantly we did not find that M. albescens was and M. volans are the result of cladegenic events that paraphyletic across much of neotropical Myotis. We also occurred during the transition of Myotis from the Old found that M. nigricans is monphlyletic in the nuclear World to the New World. It is important to remember that tree, but paraphyletic in the mitochondrial tree. These results from M. nigricans and M. albescens are interesting this interpretation relies on a fairly dramatic departure but further inference is limited due to low sample sizes from the currently accepted mitochondrial relationships of these taxa. of M. volans (represented here by a single sample) to The original subgeneric taxonomy of Myotis was based other Myotis species and abandons the nuclear UCE on three morphotypes that were later shown to be the consensus tree for the most frequently recovered UCE result of convergent evolution (Ruedi and Mayer 2001). topology. In addition, our taxonomic sampling excludes If lineage-sorting affected the mitochondrial phylogeny, Myotis species from the Old World, Ethiopian clade, the it is possible that the morphotypes truly are mono- sister clade to New World Myotis (Ruedi and Mayer 2001; phyletic. However, superimposing the previous sub- Stadelmann et al. 2007; Lack et al. 2010; Ruedi et al. 2013). generic/morphological classification onto the species Taxa from the Ethiopian clade are needed to properly tree shows an interspersed distribution of morphotypes root the New World Myotis clade and understand its throughout the most common and consensus nuclear biogeographic origins. With these caveats in mind, the topologies (Fig. 1a,b). Many strongly supported terminal hypothesis presented here should be viewed as highly [08:41 13/7/2018 Sysbio-OP-SYSB170072.tex] Page: 243 236–250 Downloaded from https://academic.oup.com/sysbio/article/67/2/236/4095671 by DeepDyve user on 13 July 2022 244 SYSTEMATIC BIOLOGY VOL. 67 speculative. Increasing the number of sampled Myotis to reflect the true species tree as any single-UCE locus lineages will shed additional light on this hypothesis. chosen at random. Phenomena such as lineage sorting, Other taxa with conflicting positions between marker reticulation, and introgression have likely influenced types include the M. lucifugus + M. occultus clade, M. the genomes of Myotis and should be accounted for in fortidens, and M. vivesi. In general, these relationships subsequent work. It is possible that the Myotis radiation are characterized by very short branches (Fig. 1d) and is more accurately reflected as a hard polytomy or a are the most likely to be affected by incomplete lineage phylogenetic network rather than a strictly bifurcating sorting or limited phylogenetic information. This could phylogeny. explain the strong support with the mitochondrial tree compared with the nuclear species tree, while allowing MATERIALS AND METHODS for a number of nuclear loci to disagree with the species tree, as well. Taxon Selection There are a number of monophyletic groups iden- Taxa were selected to span the major phylogenetic tified with nuclear data (Fig 1a) that share general break points with emphasis on the Nearctic and Neo- morphologies. For example, all of the long-eared bats tropical bifurcation as recovered in previous mitochon- (M. septentrionalis, M. auriculus, M. evotis, M. thysanodes, drial phylogenies (Stadelmann et al. 2007; Ruedi et al. and M. keenii) represent a monophyletic group of higher 2013). In addition, multiple individuals morphologically elevation, forest-dwelling species that glean insects off identified as M. nigricans and M. albescens were included of surfaces (Fitch and Shump 1979; O’Farrell and Studier to test for paraphyly as suggested by Larsen et al. 1980; Warner 1982; Manning and Jones 1989; Caceres and (2012). Three Old World species of Myotis and the Barclay. 2000). The group represented by M. fortidens, outgroup, Eptesicus fuscus, were included to root phylo- M. lucifugus, and M. occultus represent a relatively genetic analyses. All field identifications were confirmed long-haired form of Myotis. While having a distinct from museum-voucher specimens. Information for all dental formula, fortidens was historically described as a specimens examined is available in Table 1. subspecies of M. lucifugus (Miller and Allen 1928), and M. occultus has alternately represented its own species or been considered a subspecies of lucifugus (Hollister Library preparation, sequencing, and processing 1909; Valdez et al. 1999; Piaggio et al. 2002). The clade Genomic DNA was extracted from 33 samples using consisting of M. keaysi, M. oxyotus, M. ruber, M. riparius, either a Qiagen DNEasy extraction kit or a phenol- and M. diminutus represents a neotropical group of chloroform/ethanol precipitation. DNA was fragmen- primarily woolly haired bats (LaVal 1973). None of these ted using the Bioruptor UCD-300 sonication device relationships are monophyletic in the consensus or most (Diagenode, Denville, NJ, USA). Libraries were prepared common mitochondrial topologies. If the mitochondrial using the Kapa Library Preparation Kit KR0453-v2.13 genome has been subjected to phenomena that obscure (Kapa Biosystems, Wilmington, MA, USA) following the true species tree then these species groups, along the manufacturer’s instructions with five minor modi- with their synapomorphic morphological features, can fications. First, we used half volume reactions. Second, be reevaluated. subsequent to end repair, we added Sera-Mag Speed- beads (Thermo-Scientific, Waltham, MA, USA; prepared Conclusions according to Glenn et al. 2016) at a ratio of 2.86:1 for Relationships within Myotis, which until now have end repair cleanup. Third, we ligated universal iTru relied heavily on mitochondrial data, served as the y-yoke adapters (Glenn et al. 2016) onto the genomic basis for species identification (Puechmaille et al. 2012), DNA. Fourth, following adapter ligation, we performed evolutionary hypotheses (Simões et al. 2007), and even one postligation cleanup followed by Dual-SPRI size conservation recommendations (Boyles and Storm 2007). selection using 55 μL of speedbead buffer (22.5mM Previous studies using nuclear data have largely been PEG, 1M NaCl) and 25 μL of Speedbeads. Finally, we uninformative or utilized too few samples to draw performed a PCR at 95°C for 45 s, then 14 cycles of 9 °C definitive conclusions. Trees estimated from ~3650 for 30 s, 60°C for 30 sec, 72°C for 30 s, then 72°C for a nuclear loci and 295 phylogenetic analyses recovered 5 min final extension and a 12°C hold using iTru5 and 175 topologies, none of which are congruent with the iTru7 primers to produce Illumina TruSeqHT compatible mitochondrial phylogeny of Myotis. Conflict between libraries (Glenn et al. 2016). the mitochondrial and nuclear trees as well as among Libraries were quantified on a Qubit 2.0 (Life Techno- individual nuclear loci suggest that the Myotis radiation logies) and 83 ng from each library was added to create may have been accompanied by high levels of incomplete 5 pools of 6 or 7 libraries each. We then split the pools in lineage sorting and possible hybridization. Rather than two. One subsample was enriched for UCE loci, the other placing emphasis on the mitochondrial tree, it may be was not. UCE loci in the enriched library pools were cap- more appropriate to consider it for what it really is: a tured using Tetrapods 5K version 1 baits from MYcroar- single gene on par with a single-UCE locus, albeit one ray (Ann Arbor, MI, USA) following their MYbaits with many more phylogenetically informative charac- protocol v. 2.3.1 with overnight incubations (Faircloth ters. If true, then the mitochondrial genome is as likely et al. 2012). Enriched libraries were quantified with a [08:41 13/7/2018 Sysbio-OP-SYSB170072.tex] Page: 244 236–250 Downloaded from https://academic.oup.com/sysbio/article/67/2/236/4095671 by DeepDyve user on 13 July 2022 2018 PLATT ET AL.—PHYLOGENOMICS OF MYOTIS BATS 245 Qubit and pooled with other unrelated samples prior partition for gene B, and a partition for the overlapping to sequencing on an Illumina HiSeq 3000 to produce nucleotides of gene A and B. paired-end reads of 151 bases. The unenriched samples were sequenced on a separate run using a single lane of Assembly and Phylogenetic Analysis of Nuclear UCEs Illumina HiSeq 2500. All samples were demultiplexed Quality filtered raw sequence reads from the UCE- with Illumina’s fastq2bcl software. Reads were quality enriched libraries were assembled into contigs using the filtered by removing any potential adapter sequence and Trinity assembler (Grabherr et al. 2011) with a minimum trimming read ends once the average Phred quality over kmer coverage of 2. We used Phyluce to identify a four-base window score dropped below 20 using the those assembled contigs that were UCE loci. We also Fastx toolkit (Gordon and Hannon 2010). harvested UCE loci from E. fuscus (GCA_000308155.1), Myotis brandtii (GCA_000412655.1), M. davidii (GCA_000327345.1), and M. lucifugus (GCF_000147115.1) Assembly, Annotation, and Phylogenetic Analysis of the genome assemblies using the Phyluce package (Faircloth Mitochondrial Genome 2016). Once extracted from Trinity and genome Raw reads from the unenriched libraries were used assemblies, we aligned all UCE loci using MAFFT (Katoh to generate mitochondrial genomes via MitoBim (Hahn et al. 2002) and trimmed the aligned data with gBlocks et al. 2013). This program used MIRA (Chevreux et al. (Castresana 2000). Repetitive sequences (i.e. transpos- 1999) to map reads to a M. brandtii reference mitochon- able elements) in each alignment were identified with drial genome (Genbank accession number KT210199.1). RepeatMasker and trimmed where found. Alternative methods of mitochondrial genome assembly Each alignment was analyzed using three different were used when MitoBim assembly failed. These taxa partitioning schemes. Unpartitioned alignments were include M. albescens (TK61766), M. albescens (TK 101723), simply concatenated UCE loci treated as a single unit. M. albescens (RDS 7889), M. fortidens,M. keeni, M. melan- These alignments are referred to herein as “unparti- orhinus, M. nigricans (QCAZ 9601), M. septentrionalis, tioned.” Fully partitioned alignments were concatenated M. simus, M. velifer, and M. volans. For these samples, alignments mitochondrial genes that were partitioned by we first identified reads that were mitochondrial in locus. These alignments are referred to herein as “locus origin using BLAST searches against the M. brandtii partitioned.” Finally, PartitionFinder v1.1.1 (Lanfear et mitochondrial genome (KT210199.1). Those reads were al. 2012) was used to combine individual loci into an assembled using Trinity v2.2.0 with the “–single” option optimal partitioning scheme. The combined partition- to treat reads as unpaired. For taxa where we used ing schemes for each alignment were identified with NCBI genome assemblies to recover UCE loci in silico PartitionFinder v1.1.1 (Lanfear et al. 2012). Rather than mitochondrial genomes from genbank were used to searching for best-fit substitution models for each locus in the mitochondrial analyses GenBank as follows: M. or partition, the GTR+ model of sequence evolution brandtii (KT210199.1), E. fuscus (KF111725.1), M. lucifugus was assigned to all loci (Darriba and Posada 2015) (KP273591.1), and M. davidii (KM233172.1). except in the case of amino acid alignments where the Once assembled, each mitogenome was annotated MtMam model was used. Initial trees for PartitionFinder via MITOS (Bernt et al. 2013). Annotated genes were were generated using RAxML v7.4.1 (Stamatakis 2006) manually validated via BLAST to confirm sequence with linked branch lengths. Partitioning schemes were identity and length. Protein coding genes were checked heuristically searched using the hcluster algorithm. for stop codons using EMBOSS’s transeq program (Rice Partitioning schemes were chosen using the Bayeisan et al. 2000). When a stop codon was found, we used information criterion. the raw reads to verify the sequence. We used BWA Finally, each alignment and partitioning scheme v0.7.12 (Li and Durbin 2009) to align the reads to the was analyzed using Bayesian inference and maximum Mitobim assembled mitogenome to verify base calls likelihood phylogenetic methods. Bayesian phylogenies from Mitobim. The protein coding rRNA and tRNA were generated with the MPI version of ExaBayes genes from each assembly were aligned using MUSCLE (v1.4.1) using 4 independent runs of 4 chains each. and concatenated into three different alignments con- ExaBayes runs were terminated after 10 million gen- taining only protein coding genes, only rRNA and tRNA erations only if the average standard deviation of genes, or all protein coding and RNA genes. A fourth split frequencies was less than 0.01. The first 25% alignment of all protein coding genes was translated to of samples were discarded after which every 1000th amino acids and concatenated into a single alignment. generation was sampled. Proper sampling, post burn-in Several different partitioning schemes were examined was inspected via Tracer v1.6. (Rambaut and Drum- for each of the mitochondrial alignments. Alignments mond 2014) and effective sample sizes greater than 200 were either partitioned by gene, by codon, or by gene were considered acceptable. Posterior probability values type (rRNA and tRNA vs. protein coding). Genes were greater than 95% were considered to be significant. partitioned individually except in the instances where RaxML (v8.1.3) was used to estimate and score the two genes overlapped. These regions were partitioned maximum likelihood phylogeny with the rapid boot- separately from the individual genes resulting in three strapping option and 10,000 bootstrap replicates. We partitions for the two genes: a partition for gene A, a define strongly supported bipartitions as those present [08:41 13/7/2018 Sysbio-OP-SYSB170072.tex] Page: 245 236–250 Downloaded from https://academic.oup.com/sysbio/article/67/2/236/4095671 by DeepDyve user on 13 July 2022 246 SYSTEMATIC BIOLOGY VOL. 67 in 95–100% of bootstrap replicates and moderately sup- Partitioning strategies.—Alignments were analyzed ported bipartitions are present in 85–95% of bipartitions using three different partitioning schemes—single, (Wiens et al. 2008). locus, and combined—similar to the mitochondrial Alignment types.—Different combinations of partitioning schemes described above. Unpartitioned UCE loci were used to create unique matrices alignments were simply concatenated UCE loci treated based on 1) locus distribution, 2) locus length, as a single-genetic unit. Rather than searching for best-fit or 3) random locus sampling. The first set of substitution models for each UCE locus or partition, the alignments (locus distribution) was determined by GTR+ model of sequence evolution was assigned to all the number of taxa represented in the UCE alignments loci (Darriba et al. 2015). Initial trees for PartitionFinder (phyluce_align_get_only_loci_with_min_taxa; Faircloth were generated using RAxML v7.4.1 (Stamatakis 2006) 2016), or degree of completeness. Matrices were with linked branch lengths. Partitioning schemes were constructed using loci which were present in 100% heuristically searched using the hcluster algorithm (number of specimens, n = 37), 95% (n = 35), 85% and the best scheme was chosen using the Bayesian (n = 31), 75% (n = 27), 65% (n = 24), 55% (n = 20), information criterion. 45% (n = 16), 35% (n = 12), 25% (n = 9), and 15% Inference methods.—Phylogenetic trees were generated (n = 5) of specimens examined. These 10 groups were with three different phylogenetic inference methods nonexclusive, so a locus that was assembled in all across five different inference implementations includ- specimens (100% complete) would also be included ing concatenation and summary tree methods. with loci present in only 55% of specimens. On the Maximum likelihood trees were inferred for each other hand, a locus found in only 55% of specimens alignment and partitioning combination using RAxML would not be included in the 100% complete data v8.1.3 (Aberer et al. 2014). The best scoring (lowest - set. Each set of UCE alignments was concatenated lnL) tree from each data set was identified from 100 using phyluce_align_format_nexus_files_for_raxml random starting trees and bootstrapped 100 times using and a nexus character block was created using the the GTR+ in both cases. The autoMRE function in phyluce_align_format_nexus_files_for_raxml—charsets RAxML v8.1.3 was used to determine the need for option. These data sets then served as the basis for additional bootstrap replicates beyond the initial 100 downstream phylogenetic analyses. For example, when (Pattengale et al. 2009). A stopping criterion was set a a partitioning methodology (discussed below) was priori if the weighted Robinson–Foulds distance was less tested, it was performed for each of the 100%, 95%, 85%, than 5% in 95% of random permutations of computed etc. alignments. In addition to partitioning schemes, the bootstrap replicates (Pattengale et al. 2009). If necessary, effect of missing data was examined using Bayesian and an additional 100 bootstrap replicates were computed maximum likelihood methods. until the convergence stopping criteria were met. Finally, The second alignment criterion combined UCE loci bipartition frequencies of bootstrap replicates were of similar sizes. Pervious coalescent analyses of UCE drawn onto the best scoring tree from the initial data showed that subsampling the most informative RAxML searches for each of the respective data sets. loci can result in different topologies (Meiklejohn et Alignments based on UCE length or randomly sampled al. 2016). Rather than using coalescent based analyses, were analyzed using slightly different methods. In both we used concatenation of UCE loci to identify different cases UCE loci were partitioned individually and nodal topologies based on length. Under these assumptions, support was calculated using 100 bootstrap replicates UCE loci were sorted into 10 groups based on their using the RAxML fast-bootstrapping option. For all max- length and the predicted correlation between length imum likelihood analyses, we define strongly supported and number of informative characters was confirmed bipartitions as those present in 95–100% of bootstrap (Supplementary Material Fig. 1 available on Dryad at replicates and moderately supported bipartitions are http://dx.doi.org/10.5061/dryad.5g205). UCE loci in present in 85–95% of bipartitions (Wiens et al. 2008). the same size cohort were combined into a single Bayesian analyses were conducted using ExaBayes alignment, partitioned by locus, and analyzed with v1.4.1 (Aberer et al. 2014). For all Bayesian analyses RAxML using the methods described below. four independent runs of four chains each were run The final set of alignments was generated by random in parallel for a minimum of one million generations sampling of UCE loci. In large phylogenetic analyses, sampling every thousandth generation and applying systematic error can result in highly supported, but aGTR+ substitution model for each partition. After incorrect topologies due to compounding of nonphylo- one million generations, analyses continued until the genetic signal (Rodríguez-Ezpeleta et al. 2007). By standard deviation of the split frequency between chains randomly reducing the data set and replicating the max- was less than 0.01. The “-M 3” option was used to reduce imum likelihood analyses, we can reduce the potential the memory footprint of all ExaBayes runs. Proper effects of compounding error. Roughly 10% of the data sampling, post burn-in, was inspected via Tracer v1.6. set, 365 loci, were randomly sampled, concatenated, and (Rambaut et al. 2014). Effective sample sizes greater than partitioned by locus to create 100 new alignments, which 200 were considered acceptable. Posterior probability were then analyzed with RAxML using the methods values greater than 95% were considered significant. described below. An extended majority rule consensus tree was created [08:41 13/7/2018 Sysbio-OP-SYSB170072.tex] Page: 246 236–250 Downloaded from https://academic.oup.com/sysbio/article/67/2/236/4095671 by DeepDyve user on 13 July 2022 2018 PLATT ET AL.—PHYLOGENOMICS OF MYOTIS BATS 247 from all trees after the first 25% of trees were discarded analysis. For example, ASTRAL-II branch lengths are using TreeAnnotator v1.7.0 (Rambaut et al. 2013) and representative of coalescent units and ASTRID does not parameter estimates across all runs were calculated with calculate branch lengths. For accurate tree comparisons, Tracer v1.6 (Rambaut et al. 2014). branch lengths were stripped from all trees. Pairwise Species trees were calculated from gene trees for unweighted Robinson–Foulds distances were calcu- individual UCE loci recovered in five or more taxa using lated among all trees. Robinson–Foulds distances were the GTR+ substitution model and fast bootstrapping transformed into two dimensions using the stochastic (-f a) option in RAxML and 1000 bootstrap replicates. In CCA algorithm for nonlinear dimension reduction in general, gene trees were classified based on the degree TreeScaper v1.09 (Huang et al. 2016). Coordinates were of completeness (i.e. number of taxa represented) similar then visualized in R using hexagonal binning in the to the way we treated individuals as described above. hexbin library v1.27.1 (Lewin-Koh 2011). Nuclear and Species trees were estimated and bootstrapped using mitochondrial 50% majority rule consensus trees were three different programs. ASTRAL-II v4.10 (Mirarab et generated with PAUP v4.0a150 (Swofford 2003). al. 2015) was used to build a summary tree. Support values for bipartitions in the tree were generated from SUPPLEMENTARY MATERIAL 100 bootstrap replicates using site as well as site and locus resampling (Seo 2008). Species trees were estimated from Data available from the Dryad Digital Repository: ASTRID v1.4 (Vachaspati et al. 2015) using bionj and http://dx.doi.org/10.5061/dryad.5g205. bootstrapped for 100 replicates. SVDquartets (Chifman et al. 2014), as implemented in PAUP v4.0a150 (Swofford FUNDING 2003), was used to estimate a species trees from a random subset of 200,000 quartets and 1000 bootstrap replicates. This work was supported by the National Science UCE loci are relatively short markers with few Foundation [DEB-1355176]. Additional support was informative characters from which to build gene trees. provided by College of Arts and Sciences at Texas Tech Errors in gene tree estimation may reduce the accur- University. acy of summary methods and phylogenetic inference (Liu et al. 2009, Leaché Rannala 2011, DeGiorgio and ACKNOWLEDGMENTS Degnan 2014, Mirarab et al. 2016). We used weighted and unweighted statistical binning to combine genes We would like to thank the following museums, into compatible supergenes, increasing the number of collection managers, and collaborators for the tis- informative characters per “locus” and reducing the sue loans necessary to complete this work: Joseph phylogenetic error (Mirarab et al. 2014, Bayzid et al. Cook (Museum of Southwestern Biology), Museum of 2015). The gene trees used for the summary tree methods Vertebrate Zoology, Manuel Ruedi (Natural History described above were used rather than re-estimating Museum of Geneva), Santiago Burneo (Pontificia Uni- trees. Bifurcations supported by more than 50% of the versidad Catolica del Ecuador Museo de Zoologia), bootstrap replicates were retained for each gene tree. Heath Garner, Robert Bradley and Caleb Phillips (Texas Alignments from compatible trees were concatenated Tech University Natural Science Research Laboratory), into a single-supergene alignment. Trees for supergenes Link Olson (University of Alaska Museum of the were estimated using RAxML. The best trees for each North), and Cody Thompson and Priscilla Tucker supergene, as defined by log likelihood score, were (University of Michigan Museum of Zoology). In addi- retained from 500 searches. Bipartition support was tion, we would like to thank the Texas Tech HPCC estimated from 500 bootstrap replicates. For all analyses, (http://www.depts.ttu.edu/hpcc/) for providing the the GTR+ model of substitution was used and each computational resources necessary to complete this gene in the supergene alignment was partitioned separ- project. Finally we would like to thank the reviewers. ately. The resulting supertrees were then used for species Their efforts greatly improved this manuscript. tree estimation using ASTRAL-II. For the unweighted analysis, all supertrees were included in the pool of trees. REFERENCES For the weighted analysis, supertrees were weighted according to the number of genes combined in the Aberer A.J., Kobert K., Stamatakis A. 2014. ExaBayes: massively parallel supergene alignment. For example, if a supergene was Bayesian tree inference for the whole-genome era. Mol. Biol. Evol. a composite of six genes, the supertree was present six 31:2553–2556. Bayzid M.S., Mirarab S., Boussau B., Warnow T. 2015. Weighted times compared with a composite of five genes which statistical binning: enabling statistically consistent genome-scale would be represented only five times. Support for the phylogenetic analyses. PLoS One. 10:e0129183. weighted and unweighted species trees was estimated Bernt M., Donath A., Jühling F., Externbrink F., Florentz C., Fritzsch G., by site and locus resampling (Seo 2008) for 100 bootstrap Pütz J., Middendorf M., Stadler P.F. 2013. MITOS: improved de novo metazoan mitochondrial genome annotation. Mol. Phylogenet. replicates in ASTRAL-II. Evol. 69:313–319. Topological comparisons.—Trees recovered from all Bogdanowicz W., Piksa K., Tereba A. 2012. Hybridization hotspots at analyses were compared with each other in order to bat swarming sites. PLoS One 7:e53334. quantify the differences between topologies. Branch Boyles J.G., Storm J.J. 2007. The perils of picky eating: dietary breadth lengths have different meanings based on the type of is related to extinction risk in insectivorous bats. PLoS One 2:e672. [08:41 13/7/2018 Sysbio-OP-SYSB170072.tex] Page: 247 236–250 Downloaded from https://academic.oup.com/sysbio/article/67/2/236/4095671 by DeepDyve user on 13 July 2022 248 SYSTEMATIC BIOLOGY VOL. 67 Caceres M.C., Barclay R.M.R. 2000. Myotis septentrionalis. Mamm. Gresham C.R., Peterson D.G., Schmitz J., Pollock D.D., Haussler D., Species 634:1–4. Triplett E.W., Zhang G., Irie N., Jarvis E.D., Brochu C.A., Schmidt Carstens B.C., Dewey T.A. 2010. Species delimitation using a combined C.J., McCarthy F.M., Faircloth B.C., Hoffmann F.G., Glenn T.C., coalescent and information-theoretic approach: an example from Gabaldón T., Paten B., Ray D.A. 2014. Three crocodilian genomes North American Myotis bats. Syst Biol. 59:400–414. reveal ancestral patterns of evolution among archosaurs. Science Castresana J. 2000. Selection of conserved blocks from multiple 346:1254449. alignments for their use in phylogenetic analysis. Mol. Biol. Evol. Hahn C., Bachmann L., Chevreux B. 2013. Reconstructing mitochon- 17:540–552. drial genomes directly from genomic next-generation sequencing Chevreux B., Wetter T., Suhai S. 1999. Genome sequence assembly using reads—a baiting and iterative mapping approach. Nucleic Acids trace signals and additional sequence information. Proceedings of Res. 41:e129. the German Conference on Bioinformatics 99:45–56. Haynie M.L., Tsuchiya M.T.N., Ospina-Garcés S.M., Arroyo-Cabrales Chifman J., Kubatko L. 2014. Quartet inference from SNP data under J., Medellín R.A., Polaco O.J., Maldonado J.E. 2016. Placement of the the coalescent model. Bioinformatics 30:3317–3324. rediscovered Myotis planiceps (Chiroptera: Vespertilionidae) within Darriba D., Posada D. 2015. The impact of partitioning on the Myotis phylogeny. J. Mammal. 97:701–712. phylogenomic accuracy. bioRxiv. Available from: URL Hedtke S.M., Townsend T.M., Hillis D.M. 2006. Resolution of phylo- https://doi.org/10.1101/023978. genetic conflict in large data sets by increased taxon sampling. Syst. DeGiorgio M., Degnan J.H. 2014. Robustness to divergence time Biol. 55:522–529. underestimation when inferring species trees from estimated gene Hollister N. 1909. Two new bats from the southwestern United States. trees. Syst. Biol. 63:66–82. Proc. Biol. Soc. Wash. 22:43–44. Degnan J.H., Rosenberg N.A. 2009. Gene tree discordance, phylogen- Hosner P.A., Faircloth B.C., Glenn T.C., Braun E.L., Kimball R.T. etic inference and the multispecies coalescent. Trends Ecol. Evol. 2016. Avoiding missing data biases in phylogenomic inference: An 24:332–340. empirical study in the landfowl (Aves: Galliformes). Mol. Biol. Evol. Derti A., Roth F.P., Church G.M., Wu C.T. 2006. Mammalian 33:1110–1125. ultraconserved elements are strongly depleted among segmental Huang H., He Q., Kubatko L.S., Knowles L.L. 2010. Sources of duplications and copy number variants. Nat. Genet. 38:1216–1220. error inherent in species-tree estimation: impact of mutational and Douady C.J., Delsuc F., Boucher Y., Doolittle W.F., Douzery E.J.P. coalescent effects on accuracy and implications for choosing among 2003. Comparison of Bayesian and maximum likelihood bootstrap different methods. Syst. Biol. 59:573–583. measures of phylogenetic reliability. Mol. Biol. Evol. 20:248–254. Huang W., Zhou G., Marchand M., Ash J.R., Morris D., Van Dooren Dzeverin I. 2008. The stasis and possible patterns of selection in P., Brown J.M., Gallivan K. A., Wilgenbusch J.C. 2016. TreeScaper: evolution of a group of related species from the bat genus Myotis visualizing and extracting phylogenetic signal from sets of trees. (Chiroptera, Vespertilionidae). J. Mamm. Evol. 15:123–142. Mol. Biol. Evol. 33:3314–3316. Edwards S., Bensch S. 2009. Looking forwards or looking backwards Hudson R.R., Coyne J.A., Huelsenbeck J. 2002. Mathematical con- in avian phylogeography? A comment on Zink and Carrowclough sequences of the genealogical species concept. Evolution 56:1557– 2008. Mol. Ecol. 18:2930–2933. 1565. Faircloth B.C. 2016. PHYLUCE is a software package for the analysis of Kapusta A., Suh A., Feschotte C. 2017. Dynamics of genome size conserved genomic loci. Bioinformatics 32:786–788. evolution in birds and mammals. Proc. Natl. Acad. Sci. U. S. A. Faircloth B.C., Branstetter M.G., White N.D., Brady S.G. 2015. Target 114:E1460–E1469. enrichment of ultraconserved elements from arthropods provides Katoh K., Misawa K., Kuma K.I., Miyata T. 2002. MAFFT: a novel a genomic perspective on relationships among Hymenoptera. Mol. method for rapid multiple sequence alignment based on fast Fourier Ecol. Resour. 15:489–501. transform. Nucleic Acids Res. 30:3059–3066. Faircloth B.C., McCormack J.E., Crawford N.G., Harvey M.G., Lack J.B., Roehrs Z.P., Stanley C.E., Ruedi M., Van Den Bussche Brumfield R.T., Glenn T.C. 2012. Ultraconserved elements anchor R.A. 2010. Molecular phylogenetics of Myotis indicate familial- thousands of genetic markers spanning multiple evolutionary level divergence for the genus Cistugo (Chiroptera). J. Mammal. timescales. Syst. Biol. 61:717–726. 91:976–992. Findley J.S. 1972. Phenetic relationships among bats of the genus Lanfear R., Calcott B., Ho S.Y.W., Guindon S. 2012. Partition- Myotis. Syst. Biol. 21:31–52. Finder: combined selection of partitioning schemes and substi- Fitch J.H., Shump K.A. 1979. Myotis keenii. Mamm. Species Archive tution models for phylogenetic analyses. Mol. Biol. Evol. 29: 121:1–3. 1695–1701. Glenn T.C., Nilsen R., Kieran T.J., Finger J.W., Pierson T.W., Bentley Larsen R.J., Knapp M.C., Genoways H.H., Khan F.A.A., Larsen P.A., K.E., Hoffberg S., Louha S., Garcia-De-Leon F.J., Angel del Rio Wilson D.E., Baker R.J. 2012. Genetic diversity of neotropical Myotis Portilla M., Reed K., Anderson J.L., Meece J.K., Aggery S., Rekaya (Chiroptera: Vespertilionidae) with an emphasis on South American R., Alabady M., Belanger M., Winker K., Faircloth B.C. 2016. species. PLoS One 7:e46578. Adapterama I: universal stubs and primers for thousands of dual- LaVal R.K. 1973. A revision of the Neotropical bats of the genus Myotis indexed Illumina libraries (iTru & iNext). bioRxiv. Available from: Los Angeles County: Natural History Museum. URL https://doi.org/10.1101/049114. Leaché A.D., Rannala B. 2011. The accuracy of species tree estima- Good J.M., Vanderpool D., Keeble S., Bi K. 2015. Negligible nuclear tion under simulation: a comparison of methods. Syst. Biol. 60: introgression despite complete mitochondrial capture between two 126–137. species of chipmunks. Evolution 69:1961–1972. Leavitt D.H., Marion A.B., Hollingsworth B.D., Reeder T.W. 2017. Gordon A., Hannon G. 2010. Fastx-toolkit. Available from: URL Multilocus phylogeny of alligator lizards (Elgaria, Anguidae): http://hannonlab.cshl.edu/fastx_toolkit. testing mtDNA introgression as the source of discordant Grabherr M.G., Haas B.J., Yassour M., Levin J.Z., Thompson D.A., molecular phylogenetic hypotheses. Mol. Phylogenet. Evol. 110: Amit I., Adiconis X., Fan L., Raychowdhury R., Zeng Q., Chen Z., 104–121. Mauceli E., Hacohen N., Gnirke A., Rhind N., di Palma F., Birren Lewin-Koh N. 2011. Hexagon binning: an overview. Available from: B.W., Nusbaum C., Lindblad-Toh K., Friedman N., Regev A. 2011. URL ftp://ftp.naist.jp/pub/lang/R/CRAN/web/packages/ Full-length transcriptome assembly from RNA-Seq data without a hexbin/vignettes/hexagon_binning.pdf. reference genome. Nat. Biotech. 29:644–652. Li G., Davis B.W., Eizirik E., Murphy W.J. 2016. Phylogenomic evidence Green R.E., Braun E.L., Armstrong J., Earl D., Nguyen N., Hickey G., for ancient hybridization in the genomes of living cats (Felidae). Vandewege M.W., St. John J.A., Capella-Gutiérrez S., Castoe T.A., Genome Res. 26:1–11. Kern C., Fujita M.K., Opazo J.C., Jurka J., Kojima K.K., Caballero Li H., Durbin R. 2009. Fast and accurate short read alignment with J., Hubley R.M., Smit A.F., Platt R.N., Lavoie C.A., Ramakodi M.P., Burrows–Wheeler transform. Bioinformatics 25:1754–1760. Finger J.W., Suh A., Isberg S.R., Miles L., Chong A.Y., Jaratlerdsiri Liu L., Edwards S.V. 2015. Comment on “Statistical binning enables W., Gongora J., Moran C., Iriarte A., McCormack J., Burgess an accurate coalescent-based estimation of the avian tree”. Science S.C., Edwards S.V., Lyons E., Williams C., Breen M., Howard J.T., 350:171. [08:41 13/7/2018 Sysbio-OP-SYSB170072.tex] Page: 248 236–250 Downloaded from https://academic.oup.com/sysbio/article/67/2/236/4095671 by DeepDyve user on 13 July 2022 2018 PLATT ET AL.—PHYLOGENOMICS OF MYOTIS BATS 249 Liu L., Yu L., Kubatko L., Pearl D.K., Edwards S.V. 2009. Coalescent Rambaut A., Drummond A. 2013. TreeAnnotator v1. 7.0. Available methods for estimating phylogenetic trees. Mol. Phylogenet. Evol. from: URL http://beast.bio.ed.ac.uk. 53:320–328. Rambaut A., Suchard M., Xie D., Drummond A. 2014. Tracer v1.6. Manning R.W., Jones J.K. 1989. Myotis evotis. Mamm. Species Archive Available from: URL http://beast.bio.ed.ac.uk. 329:1–5. Ray D.A., Feschotte C., Pagan H.J.T., Smith J.D., Pritham E.J., McCormack J.E., Harvey M.G., Faircloth B.C., Crawford N.G., Glenn Arensburger P., Atkinson P.W., Craig N.L. 2008. Multiple waves of T.C., Brumfield R.T. 2013. A Phylogeny of birds based on over recent DNA transposon activity in the bat, Myotis lucifugus. Genome 1,500 loci collected by target enrichment and high-throughput Res. 18:717–728. sequencing. PLoS One 8:e54848. Ray D.A., Pagan H.J.T., Thompson M.L., Stevens R.D. 2007. Bats with McGee M.D., Faircloth B.C., Borstein S.R., Zheng J., Darrin Hulsey C., hATs: evidence for recent DNA transposon activity in genus Myotis. Wainwright P.C., Alfaro M.E. 2016. Replicated divergence in cichlid Mol. Biol. Evol. 24:632–639. radiations mirrors a major vertebrate innovation. Proc. R. Soc. Lond. Rice P., Longden I., Bleasby A. 2000. EMBOSS: the European molecular B Biol. Sci. 283:20151413. biology open software suite. Curr. Trends 16:276–277. Meiklejohn K.A., Faircloth B.C., Glenn T.C., Kimball R.T., Braun Rodríguez-Ezpeleta N., Brinkmann H., Roure B., Lartillot N., Lang B.F., E.L. 2016. Analysis of a rapid evolutionary radiation using Philippe H. 2007. Detecting and overcoming systematic errors in Ultraconserved Elements: Evidence for a bias in some multispecies genome-scale phylogenies. Syst. Biol. 56:389–399. coalescent methods. Syst. Biol. 65:612–627. Roehrs Z.P., Lack J.B., Van Den Bussche R.A. 2010. Tribal phylogenetic Miller G.S. Jr., Allen G. 1928. The American bats of the genus Myotis relationships within Vespertilioninae (Chiroptera: Vespertilionidae) and Pixonyx. Bull. US Nat. Mus. 144:1–218. based on mitochondrial and nuclear sequence data. J. Mammal. Mirarab S., Bayzid M.S., Boussau B., Warnow T. 2014. Statistical binning 91:1073–1092. enables an accurate coalescent-based estimation of the avian tree. Ruedi M., Mayer F. 2001. Molecular systematics of bats of the genus Science 346:1250463. Myotis (Vespertilionidae) suggests deterministic ecomorphological Mirarab S., Bayzid M.S., Warnow T. 2016. Evaluating summary convergences. Mol. Phylogenet. Evol. 21:436–448. methods for multilocus species tree estimation in the presence of Ruedi M., Stadelmann B., Gager Y., Douzery E.J.P., Francis C.M., Lin incomplete lineage sorting. Syst. Biol. 65:366–380. L.-K., Guillén-Servent A., Cibois A. 2013. Molecular phylogenetic Mirarab S., Warnow T. 2015. ASTRAL-II: coalescent-based species tree reconstructions identify East Asia as the cradle for the evolution estimation with many hundreds of taxa and thousands of genes. of the cosmopolitan genus Myotis (Mammalia, Chiroptera). Mol. Bioinformatics 31:i44–i52. Phylogenet. Evol. 69:437–449. Morales A.E., Jackson N.D., Dewey T.A., O’Meara B.C., Carstens B.C. Seo T.-K. 2008. Calculating bootstrap probabilities of phylogeny using 2016. Speciation with gene flow in North American Myotis bats. Syst. multilocus sequence data. Mol. Biol. Evol. 25:960–971. Biol. 66:440–452. Simões B.F., Rebelo H., Lopes R.J., Alves P.C., Harris D.J. 2007. Patterns O’Farrell M.J., Studier E.H. 1980. Myotis thysanodes. Mamm. Species of genetic diversity within and between Myotis d. daubentonii and 137:1–5. M. d. nathalinae derived from cytochrome b mtDNA sequence data. Pagan H.J.T., Smith J.D., Hubley R.M., Ray D.A. 2010. PiggyBac-ing on Acta Chiropt. 9:379–389. a Primate genome: novel elements, recent activity and horizontal Sota T. 2002. Radiation and reticulation: extensive introgressive transfer. Genome Biol. Evol. 2:293–303. hybridization in the carabid beetles Ohomopterus inferred Pattengale N.D., Alipour M., Bininda-Emonds O.R.P., Moret B.M.E., from mitochondrial gene genealogy. Population Ecol. 44: Stamatakis A. 2009. How many bootstrap replicates are necessary? 0145–0156. In: Batzoglou S, editor. Research in Computational Molecular Stadelmann B., Lin L.K., Kunz T.H., Ruedi M. 2007. Molecular phylo- Biology: 13th Annual International Conference, RECOMB 2009, geny of New World Myotis (Chiroptera, Vespertilionidae) inferred Tucson, AZ, USA, May 18-21, 2009. Proceedings. Berlin, Heidelberg: from mitochondrial and nuclear DNA genes. Mol. Phylogenet. Evol. Springer Berlin Heidelberg, p. 184–200. 43:32–48. Piaggio A.J., Valdez E.W., Bogan M.A., Spicer G.S. 2002. Systematics Stamatakis A. 2006. RAxML-VI-HPC: maximum likelihood-based of Myotis occultus (Chiroptera: Vespertilionidae) inferred from phylogenetic analyses with thousands of taxa and mixed models. sequences of two mitochondrial genes. J. Mammal. 83:386–395. Bioinformatics 22:2688–2690. Platt R.N., Mangum S.F., Ray D.A. 2016. Pinpointing the vesper bat Swofford D.L. 2003. PAUP*. Phylogenetic analysis using parsimony transposon revolution using the Miniopterus natalensis genome. (* and other methods). Version 4. http://paup.sc.fsu.edu/ Mob. DNA 7:12. Thomas J., Sorourian M., Ray D., Baker R.J., Pritham E.J. 2011. The Platt R.N., Vandewege M.W., Kern C., Schmidt C.J., Hoffmann F.G., limited distribution of Helitrons to vesper bats supports horizontal Ray D.A. 2014. Large numbers of novel miRNAs originate from DNA transfer. Gene 474:52–58. transposons and are coincident with a large species radiation in bats. Vachaspati P., Warnow T. 2015. ASTRID: Accurate Species TRees from Mol. Biol. Evol. 31:1536–1545. Internode Distances. BMC Genomics 16:S3. Platt R.N., Zhang Y., Witherspoon D.J., Xing J., Suh A., Keith M.S., Valdez E.W., Choate J.R., Bogan M.A., Yates T.L. 1999. Taxonomic status Jorde L.B., Stevens R.D., Ray D.A. 2015. Targeted capture of of Myotis occultus. J. Mammal. 80:545–552. phylogenetically informative ves SINE insertions in genus myotis. Warner R.M. 1982. Myotis auriculus. Mamm. Species 191:1–3. Genome Biol. Evol. 7:1664–1675. Wiens J.J., Kuczynski C.A., Smith S.A., Mulcahy D.G., Sites J.J.W., Pritham E.J., Feschotte C. 2007. Massive amplification of rolling-circle transposons in the lineage of the bat Myotis lucifugus. Proc. Natl. Townsend T.M., Reeder T.W. 2008. Branch lengths, support, and Acad. Sci. U. S. A. 104:1895–1900. congruence: testing the phylogenomic approach with 20 nuclear Puechmaille S.J., Allegrini B., Boston E.S.M., Dubourg-Savage M.-J., loci in snakes. Syst. Biol. 57:420–431. Evin A., Knochel A., Le Bris Y., Lecoq V., Lemaire M., Rist D., Teeling Willis S.C., Farias I.P., Ortí G. 2014. Testing mitochondrial capture and E.C. 2012. Genetic analyses reveal further cryptic lineages within the deep coalescence in Amazonian cichlid fishes (Cichlidae: Cichla). Myotis nattereri species complex. Mammal. Biol. 77:224–228. Evolution 68:256–268. [08:41 13/7/2018 Sysbio-OP-SYSB170072.tex] Page: 249 236–250 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Systematic Biology Oxford University Press

Conflicting Evolutionary Histories of the Mitochondrial and Nuclear Genomes in New World Myotis Bats

Loading next page...
1
 
/lp/ou_press/conflicting-evolutionary-histories-of-the-mitochondrial-and-nuclear-CamUO8d0TV

References (88)

Publisher
Oxford University Press
Copyright
Copyright © 2022 Society of Systematic Biologists
ISSN
1063-5157
eISSN
1076-836X
DOI
10.1093/sysbio/syx070
pmid
28945862
Publisher site
See Article on Publisher Site

Abstract

Downloaded from https://academic.oup.com/sysbio/article/67/2/236/4095671 by DeepDyve user on 13 July 2022 Syst. Biol. 67(2):236–249, 2018 © The Author(s) 2017. Published by Oxford University Press, on behalf of the Society of Systematic Biologists. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. For Permissions, please email: journals.permissions@oup.com DOI:10.1093/sysbio/syx070 Advance Access publication August 26, 2017 Conflicting Evolutionary Histories of the Mitochondrial and Nuclear Genomes in New World Myotis Bats 1 2 1 3 3 ROY N. PLATT II ,BRANT C. FAIRCLOTH ,KEVIN A. M. SULLIVAN ,TROY J. KIERAN ,TRAVIS C. GLENN ,MICHAEL W. 1 4 1 5 1,∗ VANDEWEGE ,THOMAS E. LEE Jr. ,ROBERT J. BAKER ,RICHARD D. STEVENS , AND DAVID A. RAY 1 2 Department of Biological Sciences, Texas Tech University, 2901 Main St, Lubbock, TX, USA; Department of Biological Sciences and Museum of Natural Science, Louisiana State University, 202 Life Science Building, Baton Rouge, LA, USA; Department of Environmental Health Science, University of Georgia, 206 Environmental Health Sciences Building, Athens, GA, USA; Department of Biology, Abilene Christian University, 1600 Campus Ct. Abilene, TX, USA; and Natural Resource Management, Texas Tech University, 2901 Main St, Lubbock, TX, USA Correspondence to be sent to: Department of Biological Sciences, Texas Tech University, Lubbock, TX, USA; E-mail: david.4.ray@gmail.com. Received 17 March 2017; reviews returned 31 July 2017; accepted 15 August 2017 Associate Editor: Matthew Hahn Abstract.—The rapid diversification of Myotis bats into more than 100 species is one of the most extensive mammalian radiations available for study. Efforts to understand relationships within Myotis have primarily utilized mitochondrial markers and trees inferred from nuclear markers lacked resolution. Our current understanding of relationships within Myotis is therefore biased towards a set of phylogenetic markers that may not reflect the history of the nuclear genome. To resolve this, we sequenced the full mitochondrial genomes of 37 representative Myotis, primarily from the New World, in conjunction with targeted sequencing of 3648 ultraconserved elements (UCEs). We inferred the phylogeny and explored the effects of concatenation and summary phylogenetic methods, as well as combinations of markers based on informativeness or levels of missing data, on our results. Of the 294 phylogenies generated from the nuclear UCE data, all are significantly different from phylogenies inferred using mitochondrial genomes. Even within the nuclear data, quartet frequencies indicate that around half of all UCE loci conflict with the estimated species tree. Several factors can drive such conflict, including incomplete lineage sorting, introgressive hybridization, or even phylogenetic error. Despite the degree of discordance between nuclear UCE loci and the mitochondrial genome and among UCE loci themselves, the most common nuclear topology is recovered in one quarter of all analyses with strong nodal support. Based on these results, we re-examine the evolutionary history of Myotis to better understand the phenomena driving their unique nuclear, mitochondrial, and biogeographic histories. [Summary tree methods; concatenation; vespertilionidae; phylogenomics; UCE; ultraconserved elements.] The bat genus Myotis (Order Chiroptera, Family Ves- single-genetic unit that is susceptible to evolutionary pertilionidae) comprises more than 100 species that ori- processes that may cause its history to diverge from ginated during the last 10–15 million years (Stadelmann the history of the species (Edwards and Bensch 2009). et al. 2007), making it one of the most successful extant The most widely accepted phylogenies of Myotis rely mammalian radiations. Myotis are distributed world- heavily on mitochondrial data and even phylogenies wide, excluding polar regions, and generally share a containing nuclear data demonstrate an over reliance on similar ecological niche: aerial insectivory. Myotis species mitochondrial markers for resolution, likely swamping often exhibit little morphological differentiation and, as out potentially conflicting signals from the nuclear a result, the rate of cryptic speciation within the genus is genome. For example alignments of the nuclear RAG2 thought to be high. For example, specimens identified and mitochondrial cytb contained 162 and 560 variable as M. nigricans and M. albescens form multiple para- characters, respectively (Stadelmann et al. 2007). Despite phyletic lineages distributed throughout the phylogeny these concerns, we find ourselves in a situation where of Neotropical Myotis (Larsen et al. 2012). Confounding our understanding of evolutionary relationships, our matters, the morphological variation that exists is often basis for conservation strategies, and our biogeographic a poor indicator of species-level relationships. Early clas- hypotheses are all founded, almost entirely, on mito- sifications of Myotis identified three major morphotypes chondrial phylogenies. (Findley 1972). Subsequent phylogenetic analyses of the Targeted sequencing of ultraconserved elements mitochondrial cytochrome-b (cytb) gene demonstrated (UCEs; Faircloth et al. 2012) to collect sequence data paraphyletic origins of each morphotype, suggesting from thousands of loci across the nuclear genome has frequent convergent evolution (Ruedi and Mayer 2001). resolved a number of difficult phylogenetic problems These same analyses demonstrated that geography was (e.g. see Faircloth et al. 2012; McCormack et al. 2013; a better predictor of phylogenetic relationship than mor- Green et al. 2014; Faircloth et al. 2015; McGee et al. 2016). phology (Ruedi and Mayer 2001; Stadelmann et al. 2007). Here we used UCE baits to collect ~1.4 Mbp from >3600 The ability of mitochondrial markers to resolve a nuclear loci in addition to random shotgun sequencing to well-supported topology does not guarantee that the collect full mitochondrial genomes in 37 taxa, primarily mitochondrial tree represents the species tree (e.g. see representing New World (NW) Myotis. Rather than Willis et al. 2014; Li et al. 2016; Leavitt et al. 2017). The analyzing each sequence data set once, or a handful of lack of recombination and uniparental inheritance of times, we chose to use a range of sampling, partitioning, the mitochondrion means that it is transmitted as a and inference methods to fully explore phylogenetic tree [08:41 13/7/2018 Sysbio-OP-SYSB170072.tex] Page: 236 236–250 Downloaded from https://academic.oup.com/sysbio/article/67/2/236/4095671 by DeepDyve user on 13 July 2022 2018 PLATT ET AL.—PHYLOGENOMICS OF MYOTIS BATS 237 space of the nuclear and mitochondrial genomes. We recovered six times, all by analysis of the translated recovered 294 trees representing 175 distinct topologies protein coding genes (Fig. 1a). Despite being the most for the nuclear UCE data and 28 trees representing common topology, bipartition frequencies and clade 14 distinct mitochondrial topologies. Our results show probability values were lower than nucleotide based data that, despite the range of trees recovered from each sets, as would be expected with more slowly evolving marker type, nuclear and mitochondrial markers occupy amino acid sequences. Analysis of RNA coding genes distinct regions in tree space. Given that the nuclear and recovered two topologies depending on the method of mitochondrial trees are distinct from one another it is phylogenetic inference. The remainder of the analyses necessary to reevaluate hypotheses made based solely recovered various trees without any pattern relating to on the mitochondrial phylogeny. partitioning scheme or inference method. A majority rule consensus tree for all 28 analyses is shown in Figure 1c. The mitochondrial consensus tree reflects pre- viously reported mitochondrial phylogenies for Myotis. RESULTS The New World Myotis are generally monophyletic We used targeted sequencing of UCEs in 37 individu- splitting into Nearctic and Neotropical clades with the als (Table 1) to collect sequence data from 3648 nuclear Old World taxon, M. brandtii, sister to the Neotropical loci which we assembled into concatenated alignments clade. Ambiguity, represented by polytomies, is found as large as 1.37 Mb. In addition, we assembled mito- only in terminal relationships. chondrial genomes for most taxa from random shotgun sequencing data. We then used the data to infer the phylogenetic history of New World Myotis using a range UCE Assembly and Alignment, and Phylogenetic Inference of locus sampling strategies, alignment partitioning We averaged 3.29 million reads per sample after methods, and phylogenetic inference methods. demultiplexing reads from the UCE-enriched, sequen- cing pool. These reads were assembled into an aver- age of 5778 contigs per sample (min=1562 M. nyctor, Mitochondrial Assembly, Alignment, and max=11,784 M. nigricans ). Recovery rates for UCE loci Phylogenetic Inference varied across taxa. Of the 5500 loci in the Amniote We generated an average of 5,498,440 paired reads per probe set, we successfully recovered 3898 UCE loci, sample (min = 2,2370,68, max = 9,961,348) of random 3648 loci from five or more samples and 212 loci in all DNA sequence data for mitochondrial genome assembly. 37 samples (Table 2). On average, 3332 UCE loci were Assemblies varied in quality. Some were almost entirely recovered per sample, ranging from 1106 (M. nyctor) complete while others were missing small regions. We to 4008 (M. keaysi). Repetitive sequences, identified via were not able to assemble the entire mitochondrial RepeatMasker searches, were minimal, occupying less control region for any of the samples. We found three than 0.02% of sites across all UCE alignments. premature stop codons in mtDNA protein coding genes. In all, 294 different combinations of alignment, par- Subsequent manual alignment and validation suggested titioning, and inference method were used to identify that these regions were miscalled by MitoBim, and we 175 unique nuclear topologies. The most common tree corrected the errors prior to analysis. Sequence coverage was recovered 45 times (Fig. 1b). All nodes were highly- of the mitochondrial genomes averaged 58× (range supported based on clade probability values and all >1×–297×). but three nodes were highly supported with maximum The 37 mitochondrial protein coding, rRNA and tRNA likelihood bootstrap replicates. A consensus tree (Fig. 1d) genes were concatenated into an alignment of 15,520 of the 294 nuclear topologies resolves New World Myotis bp containing 5007 informative characters. rRNA and (and M. brandtii) as a single polytomy that excludes tRNA genes were concatenated into an alignment con- M. volans and represents a departure from previous taining 4157 bp containing 862 parsimony informative mitochondrial phylogenies, including the one recovered characters. Protein coding genes were concatenated into herein. Comparison of the mitochondrial and nuclear an alignment containing 11,363 bp or 3770 amino acids consensus trees (Fig. 1c,d) reveals substantial differences (aa) and containing 4145 or 509 parsimony informative in the relationships within the genus. characters, respectively. For the alignment containing As noted above, we were concerned that phylogenetic all protein coding, rRNA and tRNA genes, 30 samples error, sampling bias or other methodological factors were  90% complete, and alignments for five samples may be driving the conflict between the mitochondrial were 68–84% complete. Only 21% and 50% of nucleotide and nuclear trees. To address this, we investigated positions were present in the M. albescens (TK 61766) the effects of matrix composition (or completeness) and M. levis alignments, respectively. on phylogenetic inference by generating 10 alignments When considering alignment type, partitioning having varying levels of matrix completeness. Nine strategy, and phylogenetic inference method we ana- alignments were composed of loci with 15 to 95% (at 10% lyzed the mitochondrial data 28 ways as described intervals) of taxa. The tenth alignment was composed in the Methods section. These analyses recovered 14 of all loci which were recovered in 100% of specimens distinct topologies, the most common of which was examined. For example, all loci that were represented [08:41 13/7/2018 Sysbio-OP-SYSB170072.tex] Page: 237 236–250 Downloaded from https://academic.oup.com/sysbio/article/67/2/236/4095671 by DeepDyve user on 13 July 2022 238 SYSTEMATIC BIOLOGY VOL. 67 TABLE 1. Specimens examined Species Name herein Museum identification Number of Short read Mitochondrial genome or accession number UCE loci archive accession accession Eptesicus fuscus fuscus GCA_000308155.1 2467 NA KF111725.1 Eptesicus fuscus fuscus TK 178736 2849 SAMN06141755 MF143474 Myotis albescens albescens RDS 7889 2764 SAMN06141756 MF143471 Myotis albescens albescens QCAZ 9157 1185 SAMN06141757 MF143497 Myotis albescens albescens TK 61766 2872 SAMN06141758 MF143470 Myotis ater ater M4430 2774 SAMN06141760 MF143493 Myotis auriculus auriculus MSB 40883 2229 SAMN06141761 MF143486 Myotis brandtii brandtii GCA_000412655.1 2446 NA KT210199.1 Myotis californicus californicus UMMZ 175828 2948 SAMN06141763 MF143469 Myotis davidii davidii GCA_000327345.1 2450 NA KM233172.1 Myotis diminutus diminutus QCAZ 9168 3078 SAMN06141777 MF143481 Myotis dominicensis dominicensis TK 15624 2576 SAMN06141764 MF143467 Myotis evotis evotis MSB 47323 2586 SAMN06141765 MF143468 Myotis fortidens fortidens MSB 54941 2791 SAMN06141766 MF143483 Myotis horsfieldii horsfieldii MHNG 1926.039 3017 SAMN06141767 MF143494 Myotis keaysi keaysi TK 13525 3195 SAMN06141768 MF143477 Myotis keenii keenii UAM 113849 2723 SAMN06141769 MF143472 Myotis leibii leibii TK 24872 3119 SAMN06141770 MF143488 Myotis levis levis RDS 7781 2538 SAMN06141771 MF143482 Myotis lucifugus lucifugus GCA_000147115.1 2429 NA KP273591.1 Myotis lucifugus lucifugus MSB 46679 2736 SAMN06141772 MF143491 Myotis melanorhinus melanorhinus TK 193888 3177 SAMN06141774 MF143489 Myotis nigricans nigricans QCAZ 9601 2854 SAMN06141775 MF143495 Myotis nigricans nigricans TK 194145 3159 SAMN06141776 MF143484 Myotis nyctor nyctor TK 151413 856 SAMN06141773 MF143498 Myotis occultus occultus MSB 121995 2957 SAMN06141778 MF143490 Myotis oxyotus oxyotus UMMZ RCO1013 3106 SAMN06141779 MF143479 Myotis riparius riparius TK 22688 2924 SAMN06141782 MF143473 Myotis riparius riparius TK 101723 2990 SAMN06141759 MF143475 Myotis riparius riparius TK 145199 2890 SAMN06141780 MF143480 Myotis ruber ruber MVZ 185692 2757 SAMN06141781 MF143478 Myotis septentrionalis septentrionalis TK 194420 2916 SAMN06141762 MF143487 Myotis thysanodes thysanodes 07LEP 2821 SAMN06141783 MF143492 Myotis velifer velifer MSB 70877 2704 SAMN06141784 MF143499 Myotis vivesi vivesi MSB 42658 2469 SAMN06141785 MF143476 Myotis volans volans MSB 40886 2819 SAMN06141786 MF143496 Myotis yumanensis yumanensis TK 194144 2589 SAMN06141787 MF143485 Samples beginning with “GCA_” represent genome assemblies available through NCBI. Species with more than one sample are designated with a superscript. Specimens derived from whole genome alignments are designated with a superscript “G.” MSB = Museum of Southwestern Biology; MVZ = Museum of Vertebrate Zoology; MHNG = Natural History Museum of Geneva; QCAZ = Pontificia Universidad Catolica del Ecuador Museo de Zoologia; TK = Texas Tech University Natural Science Research Laboratory; UAM = University of Alaska Museum of the North; UMMZ = University of Michigan Museum of Zoology. in at least 15% of species (n  5) were included generations. Unfortunately, computational limits forced taxa in the alignment (n = 3648). In the “100%” matrix, us to abandon the individually partitioned, Bayesian loci only loci which were identified in all species (n = analyses. In general, the same alignment produced taxa 37) were included (n = 212). Alignments were parti- the same topology regardless of inference method or loci tioned using three schemes: unpartitioned, individually partitioning scheme with the only exception being the partitioned by locus, or combined partitions. Align- terminal relationships of the M. levis/M. albescens clade ment lengths, numbers of informative characters and in the combined versus unpartitioned Bayesian analysis number of partitions identified by PartitionFinder are of the 100% complete data matrix. available in Table 2. Resulting alignments were analyzed We also tested how locus length could affect phylo- with Bayesian, maximum likelihood, and coalescent genetic discordance. To do so, trees were generated from methods using RaxML, ExaBayes, ASTRID, ASTRAL-II, data matrices incorporating UCE loci of differing lengths and SVDquartets as described in the Methods section. (Hosner et al. 2016). All 3648 loci were grouped into Bootstrap topologies stabilized in each alignment and 10 separate bins based on locus length. The number partition combination within 150 replicates and all of informative characters per bin ranged from 1115 to Bayesian runs converged in less than ten thousand 6995 and the number of informative characters was [08:41 13/7/2018 Sysbio-OP-SYSB170072.tex] Page: 238 236–250 Downloaded from https://academic.oup.com/sysbio/article/67/2/236/4095671 by DeepDyve user on 13 July 2022 2018 PLATT ET AL.—PHYLOGENOMICS OF MYOTIS BATS 239 FIGURE 1. Comparison of nuclear and mitochondrial phylogenetic trees of Myotis. The most common mitochondrial (a) and nuclear (b) topologies. The mitochondrial topology arises from Bayesian and maximum likelihood analysis of the amino acid sequences portioned based on the PartitonFinder recommendations. The nuclear topology is from the 55% complete data matrix partitioned based on the recommendations of PartitionFinder. For clarity, bootstrap values greater than 95 and clade probability values greater than 0.95 have been omitted from the nuclear consensus tree. Majority-rule consensus trees from 28 mitochondrial (c) and 294 nuclear (d) topologies. Values above the branches in the consensus trees (c,d) are bipartition frequencies for that clade across all nuclear or mitochondrial topologies. Conflicting tips between data types (consensus nuclear vs. consensus mitochondrial) are indicated with red lines between the topologies. Biogeographic regions are color coded, as are subgeneric classifications based on morphotypes, as defined by Findley (1972). Species with more than one sample are designated with a superscript that is referenced in Table 1. Specimens derived from whole genome alignments are designated with a superscript “G.” [08:41 13/7/2018 Sysbio-OP-SYSB170072.tex] Page: 239 236–250 Downloaded from https://academic.oup.com/sysbio/article/67/2/236/4095671 by DeepDyve user on 13 July 2022 240 SYSTEMATIC BIOLOGY VOL. 67 TABLE 2. Matrix composition analysis Matrix ID (% taxa Min. num. Number of loci Parsimony- Variable Alignment Optimum with data from locus) of taxa characters informative uninformative length partitions 15 5 3648 38,718 62,588 13,77,262 31 25 9 3379 37,894 57,672 12,60,248 37 35 12 3232 37,259 56,453 12,27,093 34 45 16 3064 36,539 54,782 11,87,492 31 55 20 2890 35,284 53,200 11,44,471 27 65 24 2668 34,373 51,148 10,91,620 27 75 27 2481 33,031 48,732 10,41,099 20 85 31 2034 29,179 41,711 903,903 16 95 35 1193 18,288 24,189 575,321 15 100 37 212 3480 4778 112,125 6 General alignment information. For a subset of analyses a series of alignments were generated based on the number of taxa represented by each locus. Thirty-seven taxa were examined so an alignment with all 37 taxa was considered 100% complete. Parsimony-informative characters make up a small portion of the total alignment. The optimum partitioning scheme was calculated with PartitionFinder. correlated with average locus length (Supplementary concatenation methods only recovered the same tree in Material Fig. 1 available on Dryad at http://dx.doi.org/ 1 of 100 attempts. 10.5061/dryad.5g205). On average, only 2.6% of char- Finally, we used weighted and unweighted statist- acters in each bin were parsimony-informative (Sup- ical binning to combine individual gene trees into plementary Material Fig. 1 available on Dryad at supergenes, estimate the supergene phylogeny, and then http://dx.doi.org/10.5061/dryad.5g205). Each of the infer the species tree from the supergene trees. The ten length-based alignments recovered slightly different 3648 loci were combined into 528 binned loci with 480 topologies. Terminal relationships were generally stable bins containing 7 loci each and 48 bins containing 6 across analyses with the majority of differences between loci each. Binning methods increases the normalized topologies found in the early bifurcations of Myotis. quartet scores from an average of 0.547 with gene trees Generally, longer alignments produced well resolved to 0.672 for the binned-unweighted and 0.673 for the and similar/identical topologies with significant nodal binned-weighted supergene trees. Given the relatively support regardless of the phylogenetic method or parti- even distribution of loci into bins the negligible differ- tioning scheme used. In contrast, smaller data sets were ence in quartet/species tree discordance between the more likely to yield unique topologies. unweighted and weighted supergene trees is expected. Given the goal of fully exploring nuclear tree space Binning methods have been shown to recover incorrect and generating as many reasonable nuclear UCE based species trees under coalescent methods by altering gene topologies for comparison with the mitochondrial tree, tree frequencies. In addition, it is possible that support we decided to randomly sample a limited number for incorrect trees can actually increase with the number of nuclear UCE loci to create small alignments that of input gene trees (Liu and Edwards 2015). With these would be more likely to result in unique topologies. We concerns in mind, both binning methods recovered the therefore randomly subsampled UCE loci to create 100 same topology which was the most common nuclear unique data sets of 365 loci. Loci were concatenated in topology observed across all analyses (Fig. 1c). each replicate data set and analyzed using maximum likelihood in RAxML. Of the 100 alignments analyzed, 80 unique topologies were generated (mean Robinson– Foulds distance = 4.3). Topological Comparisons In addition to concatenation methods we explored phylogenetic tree space of the UCE data set using All mitochondrial and nuclear trees are available species tree methods. Normalized quartet scores from in Supplementary Material File 1 available on Dryad ASTRAL-II (Mirarab and Warnow 2015) analyses were (http://dx.doi.org/10.5061/dryad.5g205). When visu- consistent among all analyses, with scores ranging from alizing all topologies in tree space, nuclear trees co- 0.540 to 0.553, and the number of induced quartet localized and were distinct from mitochondrial topolo- gene trees ranging from 7,745,739 (100% complete 212 gies (Fig. 2a). Comparisons of Robinson–Foulds sym- gene trees) and 63,042,410 (15% 3648 loci). SVDquartets metrical differences show 26 symmetrical differences (Chifman and Kubatko 2014) sampled all 66,045 quartets. between the most similar nuclear and mitochondrial On average, the total weight of incompatible quartets trees (Fig. 2b). For comparison only 9 of 43,070 and 0 of was 2.84%. Similar to the concatenated analysis, we 378 pairwise nuclear versus nuclear and mitochondrial inferred coalescent-based species from the same 100 versus mitochondrial tree comparisons, Figure 2c and subsamples of 365 loci described above. Despite being d, respectively, exhibited more than 26 symmetrical generated from the same underlying data, summary, and differences. [08:41 13/7/2018 Sysbio-OP-SYSB170072.tex] Page: 240 236–250 Downloaded from https://academic.oup.com/sysbio/article/67/2/236/4095671 by DeepDyve user on 13 July 2022 2018 PLATT ET AL.—PHYLOGENOMICS OF MYOTIS BATS 241 FIGURE 2. Comparison of mitochondrial and nuclear topologies (Fig. 1a). All trees recovered from mitochondrial and nuclear analyses were visualized in tree space using multidimensional scaling of Robinson–Foulds distances. Nuclear trees (green outlines) were clustered separately from mitochondrial trees (blue outlines). Robinson–Foulds distances between nuclear and mitochondrial trees (Fig. 1b), nuclear versus nuclear trees, (Fig. 1c) and mitochondrial versus mitochondrial trees (Fig. 1d) indicate the overall distinctions between the nuclear and mitochondrial topologies. Of the 45 analyses that recovered the most frequently parameters to produce a range of reasonable nuclear and observed topology (Fig. 1b), 38 were Bayesian and mitochondrial topologies that may be more useful than RAxML searches that varied by matrix completeness overreliance on a tree resulting from one or two analyses. and partitioning scheme. The fact that these analyses Rather than rejecting concordance between the two data recovered the same topology is expected given that they types from a single analysis we took steps to analyze are not independent. For example, a RAxML analysis both data sets to generate many viable topologies to of the 25% complete data set uses 1.26 Mb of the potentially reveal hidden overlap between marker types 1.38 Mb of data present in the 15% complete matrix, masked by phylogenetic error. sharing 91% of the data. Analyses that directly varied We recovered 175 and 14 unique nuclear and mito- the alignments and/or sampled less data (e.g. randomly chondrial topologies from 294 and 28 analyses using sampling loci) were more likely to generate unique multiple methodologies, sampling strategies, and phylo- topologies than the nested analyses described above. genetic parameters. By varying our phylogenetic meth- Of the 200 analyses that randomly sampled UCE loci, ods we show that recovering a consistent topology is 164 unique topologies were observed. This implies that difficult even when using the same marker type, thus when analyses of large data sets produce well-resolved making it difficult to determine the “best” topology. trees with significant nodal support, sampling smaller Rather than choosing single trees to represent the portions of the data, may provide a mechanism for nuclear and mitochondrial data, we use consensus trees creating phylogenetic uncertainty not represented by to represent the ambiguity that is present through typical tree scoring metrics. all analyses, in addition to phylogenies conforming to the most common topology (Fig. 1). Despite our efforts to identify as many plausible, distinct nuclear and mitochondrial trees as possible we were unable to DISCUSSION recover overlapping topologies between marker types We analyzed 3648 UCE loci and mitochondrial gen- suggesting that nuclear UCE loci and the mitochondrial omes of 37 Myotis bats and outgroups. The mitochon- genomes of Myotis have distinct evolutionary histories. drial and nuclear UCE phylogenies recovered distinct A comparison of topologies in tree space illustrates the topologies (Fig. 1), whether comparing the most com- differenced in trees generated for each data set (Fig. 2a). monly recovered or consensus topologies. Previous Intramarker discordance, as measured by differences in work with UCE loci demonstrated that support for nuclear versus nuclear or mitochondrial versus mito- deep divergences varied based on the number of loci chondrial topologies is low. Pairwise comparisons of examined (McCormack et al. 2013). Further, bootstrap nuclear topologies contain an average of 8.3 symmetrical replicates and clade probability values can be inaccurate differences, mitochondrial comparisons contain an aver- metrics of nodal support (Douady et al. 2003, Hedtke age of 10.8 differences. By contrast, pairwise differences et al. 2006). We varied the input data and phylogenetic between nuclear and mitochondrial topologies contain [08:41 13/7/2018 Sysbio-OP-SYSB170072.tex] Page: 241 236–250 Downloaded from https://academic.oup.com/sysbio/article/67/2/236/4095671 by DeepDyve user on 13 July 2022 242 SYSTEMATIC BIOLOGY VOL. 67 an average of 33.3 symmetrical differences. The two most agrees with Stadelmann et al. (2007)exceptinthe similar nuclear and mitochondrial topologies contain placement of M. septentrionalis and M. auriculus. Our 26 symmetrical differences. Of the 43,072 pairwise nuclear UCE topology places M. septentrionalis and M. comparisons of nuclear topologies, only 27 contain more auriculus sister to a clade containing M. californicus, M. than 26 symmetrical differences. Likewise of the 378 leibii, and M. melanorhinus. Even though this placement pairwise comparisons of mitochondrial topologies only of M. septentrionalis and M. auriculus disagrees with one contains more than 26 symmetrical differences. Stadelmann et al. (2007) it is the same relationship These observations indicate that discordance between recovered by Haynie et al. (2016). Thus, even though marker types (Fig. 2b) is much greater than within there were some concerns and miscalled nucleotides marker type (Fig. 2c,d). All but the most divergent in the MITObim assemblies, these inaccuracies have nuclear topologies are more similar to each other than not caused the mitochondrial tree to deviate from any nuclear versus mitochondrial comparison. the expectations derived from previous work. Given Conflict between mitochondrial and nuclear data may the consistency of our mitochondrial phylogeny with be driven by error in phylogenetic estimation or may previously published and independently recovered reflect genuine discordance between the two marker topologies, we are confident that the mitochondrial types (Degnan and Rosenberg 2009; Huang et al. 2010). phylogeny we recovered here, and by others, reflects We relied on multiple tree-inference methods (e.g. the true mitochondrial tree. However, the mitochondrial summary vs. concatenation), manipulated phylogenetic topology may not adequately reflect the species history, parameters (e.g. partitioning strategy), and sampling particularly when considering the factors that cause criteria (e.g. loci sampled in all taxa) to minimize the incongruence between nuclear and mitochondrial gene impacts of phylogenetic error on the data sets. In most trees. Causes of conflicting gene trees include horizontal cases, varying parameters or methodologies generated transfer, gene duplication, introgressive hybridization, unique topologies with minor differences. In most cases and incomplete lineage sorting. The rapid radiation of unique topologies resulted from rearrangements of a this clade as well as other biological factors suggests few terminal taxa. For example, placement of M. volans that some of these phenomena are more likely to have and M. brandtii was often either sister to the remaining influenced the Myotis radiation than others. New World Myotis or as an early bifurcation between Horizontal transfer of genes is thought to be rare in the Nearctic and Neotropical clades. Myotis vivesi was at eukaryotes, but, vespertilionids in general (Thomas et times found as sister to the clade containing M. lucifugus, al. 2011; Platt et al. 2014; Platt et al. 2016), and Myotis M. occultus, and M. fortidens or as sister to the clade in particular (Pritham and Feschotte 2007; Ray et al. containing the neotropical Myotis. 2007; Ray et al. 2008; Pagan et al. 2010), have experi- Summary methods failed to recover the most com- enced horizontal transfer of DNA transposons. These mon nuclear topology presented in Figure 1b except events would not be reflected in our phylogeny because when loci were binned prior to gene tree estimation. repetitive sequences were removed prior to phylogenetic Summary methods also tended to recover more unique analyses. More generally, gene duplications could create topologies than concatenation methods when analyzing conflicting signal among individual UCE markers (e.g. data from the same gene(s). For example, the random comparing nonorthologous UCE loci), but the number sample analyses recovered 80 unique topologies using of gene duplication events would have to be very high concatenation (RAxML) and 89 unique topologies with to impact enough of the 3648 UCE loci to confound the summary methods (ASTRAL-II). This likely has to do mitochondrial and nuclear phylogenies. Further ruling with the limited number of informative characters per out gene duplication events as the dominant cause of locus and by extension limited phylogenetic signal per conflicting phylogenetic signal is the fact that such events gene tree (Supplementary Material Fig. 1 available on are likely depressed in Myotis as evidenced by their Dryad at http://dx.doi.org/10.5061/dryad.5g205). In smaller genome size (~2.2 Gb) and trend towards DNA these instances, limited phylogenetic signal per gene loss (Kapusta et al. 2017) combined with low rates of would likely lead to increased opportunity for phylogen- paralogy in UCEs in general (Derti et al. 2006). etic error in gene tree estimation. Further supporting this Introgressive hybridization and reticulation could idea, binning of compatible UCE loci may have indirectly significantly influence the phylogenies of Myotis in a increased phylogenetic signal resulting in the same topo- way that leads to conflicting signal between the nuclear logy that many of the concatenation analyses recovered. and mitochondrial genomes (Sota 2002; Good et al. No other summary/coalescent method recovered this 2015). Hybridization in bats may be relatively common topology. given their propensity to swarm at cave entrances for Previous studies that used primarily mitochondrial breeding purposes. In European Myotis, swarming has data recovered effectively the same relationships among allowed for high degrees of hybridization between M. Myotis as our mitochondrial analyses (Ruedi and Mayer brandtii, M. mystacinus, and M. alcathoe (Bogdanowicz et 2001; Stadelmann et al. 2007; Roehrs et al. 2010; Larsen al. 2012). Further, M. evotis, M. thysanodes, and M. keeni all et al. 2012; Ruedi et al. 2013; Haynie et al. 2016). experienced historical gene flow during their divergence Differences between our tree and any single previ- (Carstens and Dewey 2010; Morales et al. 2016). We could ously published trees are minor. For example, our tree also explain the differences between the mitochondrial [08:41 13/7/2018 Sysbio-OP-SYSB170072.tex] Page: 242 236–250 Downloaded from https://academic.oup.com/sysbio/article/67/2/236/4095671 by DeepDyve user on 13 July 2022 2018 PLATT ET AL.—PHYLOGENOMICS OF MYOTIS BATS 243 and nuclear UCE phylogenies if Myotis experienced relationships link species with different morphotypes. extensive incomplete lineage sorting during their radi- Based on these results, it appears that the three major ation. Two factors influence the rate of lineage sorting, morphotypes in Myotis are indeed a result of convergent the fixation rate and the speciation rate (Hudson et al. evolution, as suggested by previous work (Ruedi and 2002). Increasing the time to fixation and/or decreasing Mayer 2001; Stadelmann et al. 2007) despite the conflict- the amount of time between cladogenic events will ing mitochondrial and nuclear phylogenies recovered. increase the likelihood of incomplete lineage sorting. Among the more dramatic differences between the Myotis are generally long-lived species (Dzeverin 2008) nuclear and mitochondrial topologies is the placement of and underwent a rapid radiation between 5 and 10 M. volans and M. brandtii as sister to all New World taxa MYA (Lack et al. 2010), suggesting that Myotis species by the nuclear data. Our mitochondrial analyses place are likely to experience higher levels of lineage sorting. M. volans within a Nearctic clade and M. brandtii directly The importance of these events—introgressive hybrid- inbetween the Nearctic and Neotropical bifurcations ization and incomplete lineage sorting—in driving the (Fig. 1a,c) as has been previously reported (Stadelmann differences between the mitochondrial and nuclear et al. 2007). In contrast the consensus UCE tree (Fig. 1d) phylogenies should be further explored. places M. brandtii, a species distributed throughout the Old World, within a polytomy at the base of the New World Myotis radiation and M. volans sister to all New Evolutionary History of Myotis World Myotis (including M. brandtii). The placement of M. brandtii in the nuclear UCE consensus tree does Our understanding of relationships within Myotis is not necessarily contradict the mitochondrial phylogeny heavily biased by mitochondrial data because nuclear since it is an unresolved polytomy. However, it is markers were harder to collect and produced fewer worth noting that the most commonly recovered nuclear informative sites (Ruedi and Mayer 2001; Stadelmann et al. 2007; Lack et al. 2010; Roehrs et al. 2010; Larsen topology (Fig. 1b) places M. brandtii sister to all New et al. 2012; Ruedi et al. 2013; Haynie et al. 2016)or World Myotis (excluding M. volans). This relationship nuclear phylogenies with large numbers of makers would more closely affiliate M. brandtii with other Old contained limited numbers of taxa (Platt et al. 2015). World taxa. The placement of M. volans as sister to all Our UCE-based results indicate that nuclear trees vary New World taxa (including M. brandtii) in the most substantially from the mitochondrial tree. Given that the common nuclear tree is a significant departure from nuclear and mitochondrial trees are different, we find previous work and, at first glance, does not make as it necessary to re-evaluate Myotis in the context of the much sense in a biogeographic framework. Myotis volans nuclear data. is distributed across western and northwestern North Paraphyly of M. nigricans and M. albescens was inferred America as far as far north as Alaska. Myotis brandtii is from previous mitochondrial phylogenies. Larsen et al. distributed across much of Northern Europe and into (2012) identified a minimum of 4 and potentially 12 the extreme eastern regions of Siberia. The key may lie lineages in M. albescens and M. nigricans. Our sampling in understanding the biogeography and phylogenetics included three M. albescens and two M. nigricans, of a third species, M. gracilis, a species that we were compared with Larsen’s 17 and 29 samples. Despite unfortunately unable to include. different mitochondrial and nuclear topologies overall, Myotis gracilis, along with M. brandtii,are theonlytwo our mitochondrial and nuclear phylogeny recovered the Myotis geographically distributed in the Old World, but same paraphyletic clade of three M. albescens samples phylogenetically affiliated with the New World Myotis and M. levis. Close relationships between these taxa (Stadelmann et al. 2007). If future work verifies the sister were found in previous work and are expected. More relationship between M. brandtii and M. gracilis, then we can envision a scenario where M. gracilis, M. brandtii, importantly we did not find that M. albescens was and M. volans are the result of cladegenic events that paraphyletic across much of neotropical Myotis. We also occurred during the transition of Myotis from the Old found that M. nigricans is monphlyletic in the nuclear World to the New World. It is important to remember that tree, but paraphyletic in the mitochondrial tree. These results from M. nigricans and M. albescens are interesting this interpretation relies on a fairly dramatic departure but further inference is limited due to low sample sizes from the currently accepted mitochondrial relationships of these taxa. of M. volans (represented here by a single sample) to The original subgeneric taxonomy of Myotis was based other Myotis species and abandons the nuclear UCE on three morphotypes that were later shown to be the consensus tree for the most frequently recovered UCE result of convergent evolution (Ruedi and Mayer 2001). topology. In addition, our taxonomic sampling excludes If lineage-sorting affected the mitochondrial phylogeny, Myotis species from the Old World, Ethiopian clade, the it is possible that the morphotypes truly are mono- sister clade to New World Myotis (Ruedi and Mayer 2001; phyletic. However, superimposing the previous sub- Stadelmann et al. 2007; Lack et al. 2010; Ruedi et al. 2013). generic/morphological classification onto the species Taxa from the Ethiopian clade are needed to properly tree shows an interspersed distribution of morphotypes root the New World Myotis clade and understand its throughout the most common and consensus nuclear biogeographic origins. With these caveats in mind, the topologies (Fig. 1a,b). Many strongly supported terminal hypothesis presented here should be viewed as highly [08:41 13/7/2018 Sysbio-OP-SYSB170072.tex] Page: 243 236–250 Downloaded from https://academic.oup.com/sysbio/article/67/2/236/4095671 by DeepDyve user on 13 July 2022 244 SYSTEMATIC BIOLOGY VOL. 67 speculative. Increasing the number of sampled Myotis to reflect the true species tree as any single-UCE locus lineages will shed additional light on this hypothesis. chosen at random. Phenomena such as lineage sorting, Other taxa with conflicting positions between marker reticulation, and introgression have likely influenced types include the M. lucifugus + M. occultus clade, M. the genomes of Myotis and should be accounted for in fortidens, and M. vivesi. In general, these relationships subsequent work. It is possible that the Myotis radiation are characterized by very short branches (Fig. 1d) and is more accurately reflected as a hard polytomy or a are the most likely to be affected by incomplete lineage phylogenetic network rather than a strictly bifurcating sorting or limited phylogenetic information. This could phylogeny. explain the strong support with the mitochondrial tree compared with the nuclear species tree, while allowing MATERIALS AND METHODS for a number of nuclear loci to disagree with the species tree, as well. Taxon Selection There are a number of monophyletic groups iden- Taxa were selected to span the major phylogenetic tified with nuclear data (Fig 1a) that share general break points with emphasis on the Nearctic and Neo- morphologies. For example, all of the long-eared bats tropical bifurcation as recovered in previous mitochon- (M. septentrionalis, M. auriculus, M. evotis, M. thysanodes, drial phylogenies (Stadelmann et al. 2007; Ruedi et al. and M. keenii) represent a monophyletic group of higher 2013). In addition, multiple individuals morphologically elevation, forest-dwelling species that glean insects off identified as M. nigricans and M. albescens were included of surfaces (Fitch and Shump 1979; O’Farrell and Studier to test for paraphyly as suggested by Larsen et al. 1980; Warner 1982; Manning and Jones 1989; Caceres and (2012). Three Old World species of Myotis and the Barclay. 2000). The group represented by M. fortidens, outgroup, Eptesicus fuscus, were included to root phylo- M. lucifugus, and M. occultus represent a relatively genetic analyses. All field identifications were confirmed long-haired form of Myotis. While having a distinct from museum-voucher specimens. Information for all dental formula, fortidens was historically described as a specimens examined is available in Table 1. subspecies of M. lucifugus (Miller and Allen 1928), and M. occultus has alternately represented its own species or been considered a subspecies of lucifugus (Hollister Library preparation, sequencing, and processing 1909; Valdez et al. 1999; Piaggio et al. 2002). The clade Genomic DNA was extracted from 33 samples using consisting of M. keaysi, M. oxyotus, M. ruber, M. riparius, either a Qiagen DNEasy extraction kit or a phenol- and M. diminutus represents a neotropical group of chloroform/ethanol precipitation. DNA was fragmen- primarily woolly haired bats (LaVal 1973). None of these ted using the Bioruptor UCD-300 sonication device relationships are monophyletic in the consensus or most (Diagenode, Denville, NJ, USA). Libraries were prepared common mitochondrial topologies. If the mitochondrial using the Kapa Library Preparation Kit KR0453-v2.13 genome has been subjected to phenomena that obscure (Kapa Biosystems, Wilmington, MA, USA) following the true species tree then these species groups, along the manufacturer’s instructions with five minor modi- with their synapomorphic morphological features, can fications. First, we used half volume reactions. Second, be reevaluated. subsequent to end repair, we added Sera-Mag Speed- beads (Thermo-Scientific, Waltham, MA, USA; prepared Conclusions according to Glenn et al. 2016) at a ratio of 2.86:1 for Relationships within Myotis, which until now have end repair cleanup. Third, we ligated universal iTru relied heavily on mitochondrial data, served as the y-yoke adapters (Glenn et al. 2016) onto the genomic basis for species identification (Puechmaille et al. 2012), DNA. Fourth, following adapter ligation, we performed evolutionary hypotheses (Simões et al. 2007), and even one postligation cleanup followed by Dual-SPRI size conservation recommendations (Boyles and Storm 2007). selection using 55 μL of speedbead buffer (22.5mM Previous studies using nuclear data have largely been PEG, 1M NaCl) and 25 μL of Speedbeads. Finally, we uninformative or utilized too few samples to draw performed a PCR at 95°C for 45 s, then 14 cycles of 9 °C definitive conclusions. Trees estimated from ~3650 for 30 s, 60°C for 30 sec, 72°C for 30 s, then 72°C for a nuclear loci and 295 phylogenetic analyses recovered 5 min final extension and a 12°C hold using iTru5 and 175 topologies, none of which are congruent with the iTru7 primers to produce Illumina TruSeqHT compatible mitochondrial phylogeny of Myotis. Conflict between libraries (Glenn et al. 2016). the mitochondrial and nuclear trees as well as among Libraries were quantified on a Qubit 2.0 (Life Techno- individual nuclear loci suggest that the Myotis radiation logies) and 83 ng from each library was added to create may have been accompanied by high levels of incomplete 5 pools of 6 or 7 libraries each. We then split the pools in lineage sorting and possible hybridization. Rather than two. One subsample was enriched for UCE loci, the other placing emphasis on the mitochondrial tree, it may be was not. UCE loci in the enriched library pools were cap- more appropriate to consider it for what it really is: a tured using Tetrapods 5K version 1 baits from MYcroar- single gene on par with a single-UCE locus, albeit one ray (Ann Arbor, MI, USA) following their MYbaits with many more phylogenetically informative charac- protocol v. 2.3.1 with overnight incubations (Faircloth ters. If true, then the mitochondrial genome is as likely et al. 2012). Enriched libraries were quantified with a [08:41 13/7/2018 Sysbio-OP-SYSB170072.tex] Page: 244 236–250 Downloaded from https://academic.oup.com/sysbio/article/67/2/236/4095671 by DeepDyve user on 13 July 2022 2018 PLATT ET AL.—PHYLOGENOMICS OF MYOTIS BATS 245 Qubit and pooled with other unrelated samples prior partition for gene B, and a partition for the overlapping to sequencing on an Illumina HiSeq 3000 to produce nucleotides of gene A and B. paired-end reads of 151 bases. The unenriched samples were sequenced on a separate run using a single lane of Assembly and Phylogenetic Analysis of Nuclear UCEs Illumina HiSeq 2500. All samples were demultiplexed Quality filtered raw sequence reads from the UCE- with Illumina’s fastq2bcl software. Reads were quality enriched libraries were assembled into contigs using the filtered by removing any potential adapter sequence and Trinity assembler (Grabherr et al. 2011) with a minimum trimming read ends once the average Phred quality over kmer coverage of 2. We used Phyluce to identify a four-base window score dropped below 20 using the those assembled contigs that were UCE loci. We also Fastx toolkit (Gordon and Hannon 2010). harvested UCE loci from E. fuscus (GCA_000308155.1), Myotis brandtii (GCA_000412655.1), M. davidii (GCA_000327345.1), and M. lucifugus (GCF_000147115.1) Assembly, Annotation, and Phylogenetic Analysis of the genome assemblies using the Phyluce package (Faircloth Mitochondrial Genome 2016). Once extracted from Trinity and genome Raw reads from the unenriched libraries were used assemblies, we aligned all UCE loci using MAFFT (Katoh to generate mitochondrial genomes via MitoBim (Hahn et al. 2002) and trimmed the aligned data with gBlocks et al. 2013). This program used MIRA (Chevreux et al. (Castresana 2000). Repetitive sequences (i.e. transpos- 1999) to map reads to a M. brandtii reference mitochon- able elements) in each alignment were identified with drial genome (Genbank accession number KT210199.1). RepeatMasker and trimmed where found. Alternative methods of mitochondrial genome assembly Each alignment was analyzed using three different were used when MitoBim assembly failed. These taxa partitioning schemes. Unpartitioned alignments were include M. albescens (TK61766), M. albescens (TK 101723), simply concatenated UCE loci treated as a single unit. M. albescens (RDS 7889), M. fortidens,M. keeni, M. melan- These alignments are referred to herein as “unparti- orhinus, M. nigricans (QCAZ 9601), M. septentrionalis, tioned.” Fully partitioned alignments were concatenated M. simus, M. velifer, and M. volans. For these samples, alignments mitochondrial genes that were partitioned by we first identified reads that were mitochondrial in locus. These alignments are referred to herein as “locus origin using BLAST searches against the M. brandtii partitioned.” Finally, PartitionFinder v1.1.1 (Lanfear et mitochondrial genome (KT210199.1). Those reads were al. 2012) was used to combine individual loci into an assembled using Trinity v2.2.0 with the “–single” option optimal partitioning scheme. The combined partition- to treat reads as unpaired. For taxa where we used ing schemes for each alignment were identified with NCBI genome assemblies to recover UCE loci in silico PartitionFinder v1.1.1 (Lanfear et al. 2012). Rather than mitochondrial genomes from genbank were used to searching for best-fit substitution models for each locus in the mitochondrial analyses GenBank as follows: M. or partition, the GTR+ model of sequence evolution brandtii (KT210199.1), E. fuscus (KF111725.1), M. lucifugus was assigned to all loci (Darriba and Posada 2015) (KP273591.1), and M. davidii (KM233172.1). except in the case of amino acid alignments where the Once assembled, each mitogenome was annotated MtMam model was used. Initial trees for PartitionFinder via MITOS (Bernt et al. 2013). Annotated genes were were generated using RAxML v7.4.1 (Stamatakis 2006) manually validated via BLAST to confirm sequence with linked branch lengths. Partitioning schemes were identity and length. Protein coding genes were checked heuristically searched using the hcluster algorithm. for stop codons using EMBOSS’s transeq program (Rice Partitioning schemes were chosen using the Bayeisan et al. 2000). When a stop codon was found, we used information criterion. the raw reads to verify the sequence. We used BWA Finally, each alignment and partitioning scheme v0.7.12 (Li and Durbin 2009) to align the reads to the was analyzed using Bayesian inference and maximum Mitobim assembled mitogenome to verify base calls likelihood phylogenetic methods. Bayesian phylogenies from Mitobim. The protein coding rRNA and tRNA were generated with the MPI version of ExaBayes genes from each assembly were aligned using MUSCLE (v1.4.1) using 4 independent runs of 4 chains each. and concatenated into three different alignments con- ExaBayes runs were terminated after 10 million gen- taining only protein coding genes, only rRNA and tRNA erations only if the average standard deviation of genes, or all protein coding and RNA genes. A fourth split frequencies was less than 0.01. The first 25% alignment of all protein coding genes was translated to of samples were discarded after which every 1000th amino acids and concatenated into a single alignment. generation was sampled. Proper sampling, post burn-in Several different partitioning schemes were examined was inspected via Tracer v1.6. (Rambaut and Drum- for each of the mitochondrial alignments. Alignments mond 2014) and effective sample sizes greater than 200 were either partitioned by gene, by codon, or by gene were considered acceptable. Posterior probability values type (rRNA and tRNA vs. protein coding). Genes were greater than 95% were considered to be significant. partitioned individually except in the instances where RaxML (v8.1.3) was used to estimate and score the two genes overlapped. These regions were partitioned maximum likelihood phylogeny with the rapid boot- separately from the individual genes resulting in three strapping option and 10,000 bootstrap replicates. We partitions for the two genes: a partition for gene A, a define strongly supported bipartitions as those present [08:41 13/7/2018 Sysbio-OP-SYSB170072.tex] Page: 245 236–250 Downloaded from https://academic.oup.com/sysbio/article/67/2/236/4095671 by DeepDyve user on 13 July 2022 246 SYSTEMATIC BIOLOGY VOL. 67 in 95–100% of bootstrap replicates and moderately sup- Partitioning strategies.—Alignments were analyzed ported bipartitions are present in 85–95% of bipartitions using three different partitioning schemes—single, (Wiens et al. 2008). locus, and combined—similar to the mitochondrial Alignment types.—Different combinations of partitioning schemes described above. Unpartitioned UCE loci were used to create unique matrices alignments were simply concatenated UCE loci treated based on 1) locus distribution, 2) locus length, as a single-genetic unit. Rather than searching for best-fit or 3) random locus sampling. The first set of substitution models for each UCE locus or partition, the alignments (locus distribution) was determined by GTR+ model of sequence evolution was assigned to all the number of taxa represented in the UCE alignments loci (Darriba et al. 2015). Initial trees for PartitionFinder (phyluce_align_get_only_loci_with_min_taxa; Faircloth were generated using RAxML v7.4.1 (Stamatakis 2006) 2016), or degree of completeness. Matrices were with linked branch lengths. Partitioning schemes were constructed using loci which were present in 100% heuristically searched using the hcluster algorithm (number of specimens, n = 37), 95% (n = 35), 85% and the best scheme was chosen using the Bayesian (n = 31), 75% (n = 27), 65% (n = 24), 55% (n = 20), information criterion. 45% (n = 16), 35% (n = 12), 25% (n = 9), and 15% Inference methods.—Phylogenetic trees were generated (n = 5) of specimens examined. These 10 groups were with three different phylogenetic inference methods nonexclusive, so a locus that was assembled in all across five different inference implementations includ- specimens (100% complete) would also be included ing concatenation and summary tree methods. with loci present in only 55% of specimens. On the Maximum likelihood trees were inferred for each other hand, a locus found in only 55% of specimens alignment and partitioning combination using RAxML would not be included in the 100% complete data v8.1.3 (Aberer et al. 2014). The best scoring (lowest - set. Each set of UCE alignments was concatenated lnL) tree from each data set was identified from 100 using phyluce_align_format_nexus_files_for_raxml random starting trees and bootstrapped 100 times using and a nexus character block was created using the the GTR+ in both cases. The autoMRE function in phyluce_align_format_nexus_files_for_raxml—charsets RAxML v8.1.3 was used to determine the need for option. These data sets then served as the basis for additional bootstrap replicates beyond the initial 100 downstream phylogenetic analyses. For example, when (Pattengale et al. 2009). A stopping criterion was set a a partitioning methodology (discussed below) was priori if the weighted Robinson–Foulds distance was less tested, it was performed for each of the 100%, 95%, 85%, than 5% in 95% of random permutations of computed etc. alignments. In addition to partitioning schemes, the bootstrap replicates (Pattengale et al. 2009). If necessary, effect of missing data was examined using Bayesian and an additional 100 bootstrap replicates were computed maximum likelihood methods. until the convergence stopping criteria were met. Finally, The second alignment criterion combined UCE loci bipartition frequencies of bootstrap replicates were of similar sizes. Pervious coalescent analyses of UCE drawn onto the best scoring tree from the initial data showed that subsampling the most informative RAxML searches for each of the respective data sets. loci can result in different topologies (Meiklejohn et Alignments based on UCE length or randomly sampled al. 2016). Rather than using coalescent based analyses, were analyzed using slightly different methods. In both we used concatenation of UCE loci to identify different cases UCE loci were partitioned individually and nodal topologies based on length. Under these assumptions, support was calculated using 100 bootstrap replicates UCE loci were sorted into 10 groups based on their using the RAxML fast-bootstrapping option. For all max- length and the predicted correlation between length imum likelihood analyses, we define strongly supported and number of informative characters was confirmed bipartitions as those present in 95–100% of bootstrap (Supplementary Material Fig. 1 available on Dryad at replicates and moderately supported bipartitions are http://dx.doi.org/10.5061/dryad.5g205). UCE loci in present in 85–95% of bipartitions (Wiens et al. 2008). the same size cohort were combined into a single Bayesian analyses were conducted using ExaBayes alignment, partitioned by locus, and analyzed with v1.4.1 (Aberer et al. 2014). For all Bayesian analyses RAxML using the methods described below. four independent runs of four chains each were run The final set of alignments was generated by random in parallel for a minimum of one million generations sampling of UCE loci. In large phylogenetic analyses, sampling every thousandth generation and applying systematic error can result in highly supported, but aGTR+ substitution model for each partition. After incorrect topologies due to compounding of nonphylo- one million generations, analyses continued until the genetic signal (Rodríguez-Ezpeleta et al. 2007). By standard deviation of the split frequency between chains randomly reducing the data set and replicating the max- was less than 0.01. The “-M 3” option was used to reduce imum likelihood analyses, we can reduce the potential the memory footprint of all ExaBayes runs. Proper effects of compounding error. Roughly 10% of the data sampling, post burn-in, was inspected via Tracer v1.6. set, 365 loci, were randomly sampled, concatenated, and (Rambaut et al. 2014). Effective sample sizes greater than partitioned by locus to create 100 new alignments, which 200 were considered acceptable. Posterior probability were then analyzed with RAxML using the methods values greater than 95% were considered significant. described below. An extended majority rule consensus tree was created [08:41 13/7/2018 Sysbio-OP-SYSB170072.tex] Page: 246 236–250 Downloaded from https://academic.oup.com/sysbio/article/67/2/236/4095671 by DeepDyve user on 13 July 2022 2018 PLATT ET AL.—PHYLOGENOMICS OF MYOTIS BATS 247 from all trees after the first 25% of trees were discarded analysis. For example, ASTRAL-II branch lengths are using TreeAnnotator v1.7.0 (Rambaut et al. 2013) and representative of coalescent units and ASTRID does not parameter estimates across all runs were calculated with calculate branch lengths. For accurate tree comparisons, Tracer v1.6 (Rambaut et al. 2014). branch lengths were stripped from all trees. Pairwise Species trees were calculated from gene trees for unweighted Robinson–Foulds distances were calcu- individual UCE loci recovered in five or more taxa using lated among all trees. Robinson–Foulds distances were the GTR+ substitution model and fast bootstrapping transformed into two dimensions using the stochastic (-f a) option in RAxML and 1000 bootstrap replicates. In CCA algorithm for nonlinear dimension reduction in general, gene trees were classified based on the degree TreeScaper v1.09 (Huang et al. 2016). Coordinates were of completeness (i.e. number of taxa represented) similar then visualized in R using hexagonal binning in the to the way we treated individuals as described above. hexbin library v1.27.1 (Lewin-Koh 2011). Nuclear and Species trees were estimated and bootstrapped using mitochondrial 50% majority rule consensus trees were three different programs. ASTRAL-II v4.10 (Mirarab et generated with PAUP v4.0a150 (Swofford 2003). al. 2015) was used to build a summary tree. Support values for bipartitions in the tree were generated from SUPPLEMENTARY MATERIAL 100 bootstrap replicates using site as well as site and locus resampling (Seo 2008). Species trees were estimated from Data available from the Dryad Digital Repository: ASTRID v1.4 (Vachaspati et al. 2015) using bionj and http://dx.doi.org/10.5061/dryad.5g205. bootstrapped for 100 replicates. SVDquartets (Chifman et al. 2014), as implemented in PAUP v4.0a150 (Swofford FUNDING 2003), was used to estimate a species trees from a random subset of 200,000 quartets and 1000 bootstrap replicates. This work was supported by the National Science UCE loci are relatively short markers with few Foundation [DEB-1355176]. Additional support was informative characters from which to build gene trees. provided by College of Arts and Sciences at Texas Tech Errors in gene tree estimation may reduce the accur- University. acy of summary methods and phylogenetic inference (Liu et al. 2009, Leaché Rannala 2011, DeGiorgio and ACKNOWLEDGMENTS Degnan 2014, Mirarab et al. 2016). We used weighted and unweighted statistical binning to combine genes We would like to thank the following museums, into compatible supergenes, increasing the number of collection managers, and collaborators for the tis- informative characters per “locus” and reducing the sue loans necessary to complete this work: Joseph phylogenetic error (Mirarab et al. 2014, Bayzid et al. Cook (Museum of Southwestern Biology), Museum of 2015). The gene trees used for the summary tree methods Vertebrate Zoology, Manuel Ruedi (Natural History described above were used rather than re-estimating Museum of Geneva), Santiago Burneo (Pontificia Uni- trees. Bifurcations supported by more than 50% of the versidad Catolica del Ecuador Museo de Zoologia), bootstrap replicates were retained for each gene tree. Heath Garner, Robert Bradley and Caleb Phillips (Texas Alignments from compatible trees were concatenated Tech University Natural Science Research Laboratory), into a single-supergene alignment. Trees for supergenes Link Olson (University of Alaska Museum of the were estimated using RAxML. The best trees for each North), and Cody Thompson and Priscilla Tucker supergene, as defined by log likelihood score, were (University of Michigan Museum of Zoology). In addi- retained from 500 searches. Bipartition support was tion, we would like to thank the Texas Tech HPCC estimated from 500 bootstrap replicates. For all analyses, (http://www.depts.ttu.edu/hpcc/) for providing the the GTR+ model of substitution was used and each computational resources necessary to complete this gene in the supergene alignment was partitioned separ- project. Finally we would like to thank the reviewers. ately. The resulting supertrees were then used for species Their efforts greatly improved this manuscript. tree estimation using ASTRAL-II. For the unweighted analysis, all supertrees were included in the pool of trees. REFERENCES For the weighted analysis, supertrees were weighted according to the number of genes combined in the Aberer A.J., Kobert K., Stamatakis A. 2014. ExaBayes: massively parallel supergene alignment. For example, if a supergene was Bayesian tree inference for the whole-genome era. Mol. Biol. Evol. a composite of six genes, the supertree was present six 31:2553–2556. Bayzid M.S., Mirarab S., Boussau B., Warnow T. 2015. Weighted times compared with a composite of five genes which statistical binning: enabling statistically consistent genome-scale would be represented only five times. Support for the phylogenetic analyses. PLoS One. 10:e0129183. weighted and unweighted species trees was estimated Bernt M., Donath A., Jühling F., Externbrink F., Florentz C., Fritzsch G., by site and locus resampling (Seo 2008) for 100 bootstrap Pütz J., Middendorf M., Stadler P.F. 2013. MITOS: improved de novo metazoan mitochondrial genome annotation. Mol. Phylogenet. replicates in ASTRAL-II. Evol. 69:313–319. Topological comparisons.—Trees recovered from all Bogdanowicz W., Piksa K., Tereba A. 2012. Hybridization hotspots at analyses were compared with each other in order to bat swarming sites. PLoS One 7:e53334. quantify the differences between topologies. Branch Boyles J.G., Storm J.J. 2007. The perils of picky eating: dietary breadth lengths have different meanings based on the type of is related to extinction risk in insectivorous bats. PLoS One 2:e672. [08:41 13/7/2018 Sysbio-OP-SYSB170072.tex] Page: 247 236–250 Downloaded from https://academic.oup.com/sysbio/article/67/2/236/4095671 by DeepDyve user on 13 July 2022 248 SYSTEMATIC BIOLOGY VOL. 67 Caceres M.C., Barclay R.M.R. 2000. Myotis septentrionalis. Mamm. Gresham C.R., Peterson D.G., Schmitz J., Pollock D.D., Haussler D., Species 634:1–4. Triplett E.W., Zhang G., Irie N., Jarvis E.D., Brochu C.A., Schmidt Carstens B.C., Dewey T.A. 2010. Species delimitation using a combined C.J., McCarthy F.M., Faircloth B.C., Hoffmann F.G., Glenn T.C., coalescent and information-theoretic approach: an example from Gabaldón T., Paten B., Ray D.A. 2014. Three crocodilian genomes North American Myotis bats. Syst Biol. 59:400–414. reveal ancestral patterns of evolution among archosaurs. Science Castresana J. 2000. Selection of conserved blocks from multiple 346:1254449. alignments for their use in phylogenetic analysis. Mol. Biol. Evol. Hahn C., Bachmann L., Chevreux B. 2013. Reconstructing mitochon- 17:540–552. drial genomes directly from genomic next-generation sequencing Chevreux B., Wetter T., Suhai S. 1999. Genome sequence assembly using reads—a baiting and iterative mapping approach. Nucleic Acids trace signals and additional sequence information. Proceedings of Res. 41:e129. the German Conference on Bioinformatics 99:45–56. Haynie M.L., Tsuchiya M.T.N., Ospina-Garcés S.M., Arroyo-Cabrales Chifman J., Kubatko L. 2014. Quartet inference from SNP data under J., Medellín R.A., Polaco O.J., Maldonado J.E. 2016. Placement of the the coalescent model. Bioinformatics 30:3317–3324. rediscovered Myotis planiceps (Chiroptera: Vespertilionidae) within Darriba D., Posada D. 2015. The impact of partitioning on the Myotis phylogeny. J. Mammal. 97:701–712. phylogenomic accuracy. bioRxiv. Available from: URL Hedtke S.M., Townsend T.M., Hillis D.M. 2006. Resolution of phylo- https://doi.org/10.1101/023978. genetic conflict in large data sets by increased taxon sampling. Syst. DeGiorgio M., Degnan J.H. 2014. Robustness to divergence time Biol. 55:522–529. underestimation when inferring species trees from estimated gene Hollister N. 1909. Two new bats from the southwestern United States. trees. Syst. Biol. 63:66–82. Proc. Biol. Soc. Wash. 22:43–44. Degnan J.H., Rosenberg N.A. 2009. Gene tree discordance, phylogen- Hosner P.A., Faircloth B.C., Glenn T.C., Braun E.L., Kimball R.T. etic inference and the multispecies coalescent. Trends Ecol. Evol. 2016. Avoiding missing data biases in phylogenomic inference: An 24:332–340. empirical study in the landfowl (Aves: Galliformes). Mol. Biol. Evol. Derti A., Roth F.P., Church G.M., Wu C.T. 2006. Mammalian 33:1110–1125. ultraconserved elements are strongly depleted among segmental Huang H., He Q., Kubatko L.S., Knowles L.L. 2010. Sources of duplications and copy number variants. Nat. Genet. 38:1216–1220. error inherent in species-tree estimation: impact of mutational and Douady C.J., Delsuc F., Boucher Y., Doolittle W.F., Douzery E.J.P. coalescent effects on accuracy and implications for choosing among 2003. Comparison of Bayesian and maximum likelihood bootstrap different methods. Syst. Biol. 59:573–583. measures of phylogenetic reliability. Mol. Biol. Evol. 20:248–254. Huang W., Zhou G., Marchand M., Ash J.R., Morris D., Van Dooren Dzeverin I. 2008. The stasis and possible patterns of selection in P., Brown J.M., Gallivan K. A., Wilgenbusch J.C. 2016. TreeScaper: evolution of a group of related species from the bat genus Myotis visualizing and extracting phylogenetic signal from sets of trees. (Chiroptera, Vespertilionidae). J. Mamm. Evol. 15:123–142. Mol. Biol. Evol. 33:3314–3316. Edwards S., Bensch S. 2009. Looking forwards or looking backwards Hudson R.R., Coyne J.A., Huelsenbeck J. 2002. Mathematical con- in avian phylogeography? A comment on Zink and Carrowclough sequences of the genealogical species concept. Evolution 56:1557– 2008. Mol. Ecol. 18:2930–2933. 1565. Faircloth B.C. 2016. PHYLUCE is a software package for the analysis of Kapusta A., Suh A., Feschotte C. 2017. Dynamics of genome size conserved genomic loci. Bioinformatics 32:786–788. evolution in birds and mammals. Proc. Natl. Acad. Sci. U. S. A. Faircloth B.C., Branstetter M.G., White N.D., Brady S.G. 2015. Target 114:E1460–E1469. enrichment of ultraconserved elements from arthropods provides Katoh K., Misawa K., Kuma K.I., Miyata T. 2002. MAFFT: a novel a genomic perspective on relationships among Hymenoptera. Mol. method for rapid multiple sequence alignment based on fast Fourier Ecol. Resour. 15:489–501. transform. Nucleic Acids Res. 30:3059–3066. Faircloth B.C., McCormack J.E., Crawford N.G., Harvey M.G., Lack J.B., Roehrs Z.P., Stanley C.E., Ruedi M., Van Den Bussche Brumfield R.T., Glenn T.C. 2012. Ultraconserved elements anchor R.A. 2010. Molecular phylogenetics of Myotis indicate familial- thousands of genetic markers spanning multiple evolutionary level divergence for the genus Cistugo (Chiroptera). J. Mammal. timescales. Syst. Biol. 61:717–726. 91:976–992. Findley J.S. 1972. Phenetic relationships among bats of the genus Lanfear R., Calcott B., Ho S.Y.W., Guindon S. 2012. Partition- Myotis. Syst. Biol. 21:31–52. Finder: combined selection of partitioning schemes and substi- Fitch J.H., Shump K.A. 1979. Myotis keenii. Mamm. Species Archive tution models for phylogenetic analyses. Mol. Biol. Evol. 29: 121:1–3. 1695–1701. Glenn T.C., Nilsen R., Kieran T.J., Finger J.W., Pierson T.W., Bentley Larsen R.J., Knapp M.C., Genoways H.H., Khan F.A.A., Larsen P.A., K.E., Hoffberg S., Louha S., Garcia-De-Leon F.J., Angel del Rio Wilson D.E., Baker R.J. 2012. Genetic diversity of neotropical Myotis Portilla M., Reed K., Anderson J.L., Meece J.K., Aggery S., Rekaya (Chiroptera: Vespertilionidae) with an emphasis on South American R., Alabady M., Belanger M., Winker K., Faircloth B.C. 2016. species. PLoS One 7:e46578. Adapterama I: universal stubs and primers for thousands of dual- LaVal R.K. 1973. A revision of the Neotropical bats of the genus Myotis indexed Illumina libraries (iTru & iNext). bioRxiv. Available from: Los Angeles County: Natural History Museum. URL https://doi.org/10.1101/049114. Leaché A.D., Rannala B. 2011. The accuracy of species tree estima- Good J.M., Vanderpool D., Keeble S., Bi K. 2015. Negligible nuclear tion under simulation: a comparison of methods. Syst. Biol. 60: introgression despite complete mitochondrial capture between two 126–137. species of chipmunks. Evolution 69:1961–1972. Leavitt D.H., Marion A.B., Hollingsworth B.D., Reeder T.W. 2017. Gordon A., Hannon G. 2010. Fastx-toolkit. Available from: URL Multilocus phylogeny of alligator lizards (Elgaria, Anguidae): http://hannonlab.cshl.edu/fastx_toolkit. testing mtDNA introgression as the source of discordant Grabherr M.G., Haas B.J., Yassour M., Levin J.Z., Thompson D.A., molecular phylogenetic hypotheses. Mol. Phylogenet. Evol. 110: Amit I., Adiconis X., Fan L., Raychowdhury R., Zeng Q., Chen Z., 104–121. Mauceli E., Hacohen N., Gnirke A., Rhind N., di Palma F., Birren Lewin-Koh N. 2011. Hexagon binning: an overview. Available from: B.W., Nusbaum C., Lindblad-Toh K., Friedman N., Regev A. 2011. URL ftp://ftp.naist.jp/pub/lang/R/CRAN/web/packages/ Full-length transcriptome assembly from RNA-Seq data without a hexbin/vignettes/hexagon_binning.pdf. reference genome. Nat. Biotech. 29:644–652. Li G., Davis B.W., Eizirik E., Murphy W.J. 2016. Phylogenomic evidence Green R.E., Braun E.L., Armstrong J., Earl D., Nguyen N., Hickey G., for ancient hybridization in the genomes of living cats (Felidae). Vandewege M.W., St. John J.A., Capella-Gutiérrez S., Castoe T.A., Genome Res. 26:1–11. Kern C., Fujita M.K., Opazo J.C., Jurka J., Kojima K.K., Caballero Li H., Durbin R. 2009. Fast and accurate short read alignment with J., Hubley R.M., Smit A.F., Platt R.N., Lavoie C.A., Ramakodi M.P., Burrows–Wheeler transform. Bioinformatics 25:1754–1760. Finger J.W., Suh A., Isberg S.R., Miles L., Chong A.Y., Jaratlerdsiri Liu L., Edwards S.V. 2015. Comment on “Statistical binning enables W., Gongora J., Moran C., Iriarte A., McCormack J., Burgess an accurate coalescent-based estimation of the avian tree”. Science S.C., Edwards S.V., Lyons E., Williams C., Breen M., Howard J.T., 350:171. [08:41 13/7/2018 Sysbio-OP-SYSB170072.tex] Page: 248 236–250 Downloaded from https://academic.oup.com/sysbio/article/67/2/236/4095671 by DeepDyve user on 13 July 2022 2018 PLATT ET AL.—PHYLOGENOMICS OF MYOTIS BATS 249 Liu L., Yu L., Kubatko L., Pearl D.K., Edwards S.V. 2009. Coalescent Rambaut A., Drummond A. 2013. TreeAnnotator v1. 7.0. Available methods for estimating phylogenetic trees. Mol. Phylogenet. Evol. from: URL http://beast.bio.ed.ac.uk. 53:320–328. Rambaut A., Suchard M., Xie D., Drummond A. 2014. Tracer v1.6. Manning R.W., Jones J.K. 1989. Myotis evotis. Mamm. Species Archive Available from: URL http://beast.bio.ed.ac.uk. 329:1–5. Ray D.A., Feschotte C., Pagan H.J.T., Smith J.D., Pritham E.J., McCormack J.E., Harvey M.G., Faircloth B.C., Crawford N.G., Glenn Arensburger P., Atkinson P.W., Craig N.L. 2008. Multiple waves of T.C., Brumfield R.T. 2013. A Phylogeny of birds based on over recent DNA transposon activity in the bat, Myotis lucifugus. Genome 1,500 loci collected by target enrichment and high-throughput Res. 18:717–728. sequencing. PLoS One 8:e54848. Ray D.A., Pagan H.J.T., Thompson M.L., Stevens R.D. 2007. Bats with McGee M.D., Faircloth B.C., Borstein S.R., Zheng J., Darrin Hulsey C., hATs: evidence for recent DNA transposon activity in genus Myotis. Wainwright P.C., Alfaro M.E. 2016. Replicated divergence in cichlid Mol. Biol. Evol. 24:632–639. radiations mirrors a major vertebrate innovation. Proc. R. Soc. Lond. Rice P., Longden I., Bleasby A. 2000. EMBOSS: the European molecular B Biol. Sci. 283:20151413. biology open software suite. Curr. Trends 16:276–277. Meiklejohn K.A., Faircloth B.C., Glenn T.C., Kimball R.T., Braun Rodríguez-Ezpeleta N., Brinkmann H., Roure B., Lartillot N., Lang B.F., E.L. 2016. Analysis of a rapid evolutionary radiation using Philippe H. 2007. Detecting and overcoming systematic errors in Ultraconserved Elements: Evidence for a bias in some multispecies genome-scale phylogenies. Syst. Biol. 56:389–399. coalescent methods. Syst. Biol. 65:612–627. Roehrs Z.P., Lack J.B., Van Den Bussche R.A. 2010. Tribal phylogenetic Miller G.S. Jr., Allen G. 1928. The American bats of the genus Myotis relationships within Vespertilioninae (Chiroptera: Vespertilionidae) and Pixonyx. Bull. US Nat. Mus. 144:1–218. based on mitochondrial and nuclear sequence data. J. Mammal. Mirarab S., Bayzid M.S., Boussau B., Warnow T. 2014. Statistical binning 91:1073–1092. enables an accurate coalescent-based estimation of the avian tree. Ruedi M., Mayer F. 2001. Molecular systematics of bats of the genus Science 346:1250463. Myotis (Vespertilionidae) suggests deterministic ecomorphological Mirarab S., Bayzid M.S., Warnow T. 2016. Evaluating summary convergences. Mol. Phylogenet. Evol. 21:436–448. methods for multilocus species tree estimation in the presence of Ruedi M., Stadelmann B., Gager Y., Douzery E.J.P., Francis C.M., Lin incomplete lineage sorting. Syst. Biol. 65:366–380. L.-K., Guillén-Servent A., Cibois A. 2013. Molecular phylogenetic Mirarab S., Warnow T. 2015. ASTRAL-II: coalescent-based species tree reconstructions identify East Asia as the cradle for the evolution estimation with many hundreds of taxa and thousands of genes. of the cosmopolitan genus Myotis (Mammalia, Chiroptera). Mol. Bioinformatics 31:i44–i52. Phylogenet. Evol. 69:437–449. Morales A.E., Jackson N.D., Dewey T.A., O’Meara B.C., Carstens B.C. Seo T.-K. 2008. Calculating bootstrap probabilities of phylogeny using 2016. Speciation with gene flow in North American Myotis bats. Syst. multilocus sequence data. Mol. Biol. Evol. 25:960–971. Biol. 66:440–452. Simões B.F., Rebelo H., Lopes R.J., Alves P.C., Harris D.J. 2007. Patterns O’Farrell M.J., Studier E.H. 1980. Myotis thysanodes. Mamm. Species of genetic diversity within and between Myotis d. daubentonii and 137:1–5. M. d. nathalinae derived from cytochrome b mtDNA sequence data. Pagan H.J.T., Smith J.D., Hubley R.M., Ray D.A. 2010. PiggyBac-ing on Acta Chiropt. 9:379–389. a Primate genome: novel elements, recent activity and horizontal Sota T. 2002. Radiation and reticulation: extensive introgressive transfer. Genome Biol. Evol. 2:293–303. hybridization in the carabid beetles Ohomopterus inferred Pattengale N.D., Alipour M., Bininda-Emonds O.R.P., Moret B.M.E., from mitochondrial gene genealogy. Population Ecol. 44: Stamatakis A. 2009. How many bootstrap replicates are necessary? 0145–0156. In: Batzoglou S, editor. Research in Computational Molecular Stadelmann B., Lin L.K., Kunz T.H., Ruedi M. 2007. Molecular phylo- Biology: 13th Annual International Conference, RECOMB 2009, geny of New World Myotis (Chiroptera, Vespertilionidae) inferred Tucson, AZ, USA, May 18-21, 2009. Proceedings. Berlin, Heidelberg: from mitochondrial and nuclear DNA genes. Mol. Phylogenet. Evol. Springer Berlin Heidelberg, p. 184–200. 43:32–48. Piaggio A.J., Valdez E.W., Bogan M.A., Spicer G.S. 2002. Systematics Stamatakis A. 2006. RAxML-VI-HPC: maximum likelihood-based of Myotis occultus (Chiroptera: Vespertilionidae) inferred from phylogenetic analyses with thousands of taxa and mixed models. sequences of two mitochondrial genes. J. Mammal. 83:386–395. Bioinformatics 22:2688–2690. Platt R.N., Mangum S.F., Ray D.A. 2016. Pinpointing the vesper bat Swofford D.L. 2003. PAUP*. Phylogenetic analysis using parsimony transposon revolution using the Miniopterus natalensis genome. (* and other methods). Version 4. http://paup.sc.fsu.edu/ Mob. DNA 7:12. Thomas J., Sorourian M., Ray D., Baker R.J., Pritham E.J. 2011. The Platt R.N., Vandewege M.W., Kern C., Schmidt C.J., Hoffmann F.G., limited distribution of Helitrons to vesper bats supports horizontal Ray D.A. 2014. Large numbers of novel miRNAs originate from DNA transfer. Gene 474:52–58. transposons and are coincident with a large species radiation in bats. Vachaspati P., Warnow T. 2015. ASTRID: Accurate Species TRees from Mol. Biol. Evol. 31:1536–1545. Internode Distances. BMC Genomics 16:S3. Platt R.N., Zhang Y., Witherspoon D.J., Xing J., Suh A., Keith M.S., Valdez E.W., Choate J.R., Bogan M.A., Yates T.L. 1999. Taxonomic status Jorde L.B., Stevens R.D., Ray D.A. 2015. Targeted capture of of Myotis occultus. J. Mammal. 80:545–552. phylogenetically informative ves SINE insertions in genus myotis. Warner R.M. 1982. Myotis auriculus. Mamm. Species 191:1–3. Genome Biol. Evol. 7:1664–1675. Wiens J.J., Kuczynski C.A., Smith S.A., Mulcahy D.G., Sites J.J.W., Pritham E.J., Feschotte C. 2007. Massive amplification of rolling-circle transposons in the lineage of the bat Myotis lucifugus. Proc. Natl. Townsend T.M., Reeder T.W. 2008. Branch lengths, support, and Acad. Sci. U. S. A. 104:1895–1900. congruence: testing the phylogenomic approach with 20 nuclear Puechmaille S.J., Allegrini B., Boston E.S.M., Dubourg-Savage M.-J., loci in snakes. Syst. Biol. 57:420–431. Evin A., Knochel A., Le Bris Y., Lecoq V., Lemaire M., Rist D., Teeling Willis S.C., Farias I.P., Ortí G. 2014. Testing mitochondrial capture and E.C. 2012. Genetic analyses reveal further cryptic lineages within the deep coalescence in Amazonian cichlid fishes (Cichlidae: Cichla). Myotis nattereri species complex. Mammal. Biol. 77:224–228. Evolution 68:256–268. [08:41 13/7/2018 Sysbio-OP-SYSB170072.tex] Page: 249 236–250

Journal

Systematic BiologyOxford University Press

Published: Mar 1, 2018

There are no references for this article.