Mimicry and mitonuclear discordance in nudibranchs: New insights from exon capture phylogenomics

Abstract Phylogenetic inference and species delimitation can be challenging in taxonomic groups that have recently radiated and where introgression produces conflicting gene trees, especially when species delimitation has traditionally relied on mitochondrial data and color pattern. Chromodoris, a genus of colorful and toxic nudibranch in the Indo‐Pacific, has been shown to have extraordinary cryptic diversity and mimicry, and has recently radiated, ultimately complicating species delimitation. In these cases, additional genome‐wide data can help improve phylogenetic resolution and provide important insights about evolutionary history. Here, we employ a transcriptome‐based exon capture approach to resolve Chromodoris phylogeny with data from 2,925 exons and 1,630 genes, derived from 15 nudibranch transcriptomes. We show that some previously identified mimics instead show mitonuclear discordance, likely deriving from introgression or mitochondrial capture, but we confirm one “pure” mimic in Western Australia. Sister–species relationships and species‐level entities were recovered with high support in both concatenated maximum likelihood (ML) and summary coalescent phylogenies, but the ML topologies were highly variable while the coalescent topologies were consistent across datasets. Our work also demonstrates the broad phylogenetic utility of 149 genes that were previously identified from eupulmonate gastropods. This study is one of the first to (a) demonstrate the efficacy of exon capture for recovering relationships among recently radiated invertebrate taxa, (b) employ genome‐wide nuclear markers to test mimicry hypotheses in nudibranchs and (c) provide evidence for introgression and mitochondrial capture in nudibranchs.

problematic in taxonomic groups where species delimitation has traditionally relied on mitochondrial data and color pattern because introgression is undetectable without information from the nuclear genome and it can facilitate the evolution of mimicry (Enciso-Romero et al., 2017). This is particularly true for tropical nudibranchs, where systematic relationships remain problematic due to poor resolution, high levels of diversity, and the recent discovery of mimicry in some groups (Layton, Gosliner, & Wilson, 2018;Padula et al., 2016). Although mimicry has been extensively studied in many terrestrial taxa (e.g., butterflies, Flanagan et al., 2004;Heliconius Genome Consortium, 2012;millipedes, Marek & Bond, 2009;and frogs, Chouteau, Summers, Morales, & Angers, 2011), less is known about mimicry in marine invertebrates.
Chromodorid nudibranchs have evolved to selectively sequester toxic chemicals from their sponge prey and many display aposematic signals (e.g., Cimino & Ghiselin, 1999;Rudman, 1991). Previous work has identified a Müllerian mimicry ring in eastern Australia, where multiple chemically defended Chromodoris species appear very similar (Cheney et al., 2016), and a second quasi-Batesian system in chromodorids where some species are unpalatable but nontoxic (Winters, White, et al., 2018). The type genus, Chromodoris, is brightly colored and exhibits blue, black, yellow, and white color patterns, similar to the warning colors found in other mimetic taxa (e.g., Heliconius, Mallet, 2010). Müllerian mimicry can evolve by mutation or through introgression among species, the latter of which has been well documented in butterflies (Enciso-Romero et al., 2017;Pardo-Diaz et al.,2012), but identifying the origin of mimicry is difficult for groups lacking genomic resources. Previous studies have relied onmitochondrial data for phylogenetic reconstruction in Chromodoris, which prevents the detection of introgression, and the resulting phylogenies were poorly resolved with short branch lengths, indicating that a recent radiation has occurred (Johnson & Gosliner, 2012;Layton et al., 2018;Turner & Wilson, 2008;Wilson & Lee, 2005). Layton et al. (2018) also identified up to four distinct color patterns associated with a single mitochondrial clade that matched other congenerics, indicating the presence of mimicry in this genus. Given the lack of resolution in the existing Chromodoris phylogeny, and the evidence for mimicry and a recent radiation, new methods are needed to clarify evolutionary relationships and patterns in this complex genus.
Next-generation sequencing techniques produce large amounts of genome-wide data that can be used to improve phylogenetic resolution (e.g., Hackett et al.,2008;McCormack et al.,2013) and reduced representation techniques that target-specific sections of the genome are becoming increasingly popular for this work. Exon capture, which utilizes transcriptomes to establish a universal set of loci that can be targeted with probes, has been employed for phylogenetic reconstruction in several taxonomic groups (e.g., Bi et al.,2012;Bragg, Potter, Bi, & Moritz, 2016;Hugall, O'Hara, Hunjan, Nilsen, & Mousalli, 2016;Quattrini et al.,2018;Teasdale, Kohler, Murray, O'Hara, & Moussalli, 2016), but its utility for resolving phylogenetic relationships at shallow scales of divergence and in cases where species delimitation is problematic has rarely been tested.
Our objectives in this study were to (a) assess the efficacy of exon capture for resolving phylogenetic relationships in a group of recently radiated nudibranchs, (b) test hypotheses of mimicry in the genus, and (c) assess topological variability among different gene sets and phylogenetic methods.

| Sample collection in RNAlater
Two Chromodoris (C. magnifica, C. westraliensis), ten chromodorids from other genera, and one Actinocyclus were collected by hand using SCUBA at depths of 7-17 m from sites in San Diego

| RNA extraction, transcriptome sequencing, and de novo assembly
Total RNA was extracted from RNAlater samples using either the QiagenRNeasy kit following manufacturer's protocols or by homogenizing in TRIzol reagent for 1min with zirconium oxide beads and subsequently following manufacturer's protocols for the Zymo Direct-zol purification kits including the optional DNAse step. The quantity of RNA was determined with Nanodrop, Qubit and gel visualization. RNA purity was quantified using Nanodrop 260/280 nm absorption ratios with ratios close to 2.0 considered pure RNA and a volume of 40μl was used for multiplexed sequencing. Transcriptomes were sequenced on an Illumina HiSeq 2000 platform, pooling 5-6 libraries per lane, using 100bp paired end reads by the Australian Genome Research Facility (Melbourne). In order to enhance computational efficiency, a random subset of 20,000,000 paired reads were selected from the four transcriptomes with more than 20M reads (Actinocyclus verrucosus, Chromodoris magnifica, Felimida macfarlandi, Verconia norba) prior to assembly. Randomly subsampling reads did not affect downstream analysis because rare transcripts were not assessed in this study. Two additional nudibranch transcriptomes (Prodoris clavigera, Doris kerguelenensis) were downloaded from GenBank. Both the subsetted transcriptomes (N = 4) and the full transcriptomes (N = 11) were assembled de novo using the transcriptome pipeline in Agalma (Dunn, Howison, & Zapata, 2013)  were removed after assembly to remove weakly expressed isoforms of transcripts or incorrectly identified isomers due to sequencing errors.

| Identifying single copy homologous genes
A maximum-likelihood phylogeny of 1,126 genes from all transcriptomes (see "all nudibranchs" below) was constructed using a single partition and a GTR20 + F substitution model in IQtree v1.6.8 (Nguyen, Schmidt, von Haeseler, & Minh, 2015). This transcriptome phylogeny was used to construct a framework for interpreting gene and taxon selection in downstream analyses. Multiple methods of target gene identification were used to identify, compare, and maximize the number of phylogenetically informative genes for Chromodoris.
A workflow for target gene search and identifying exon sequences is available in Figure 1a. These were identified from three comparisons across different taxonomic levels; (a) between Chromodoris and closely related genera (five transcriptomes; "Fast 5"), (b) across the family Chromodorididae and other nudibranch outgroups (fifteen transcriptomes; "all nudibranchs"), and (c) across eupulmonate gastropods (Teasdale et al., 2016). Single copy homologous genes were identified from two different searches using the phylogeny pipeline in Agalma. The first set of homologous genes was recovered with a nucleotide search of five closely related species (C. magnifica, C. westraliensis, Doriprismatica atromarginata, Goniobranchus coi, G. fidelis).
The second set of homologous genes was then recovered with an amino acid search of all fifteen nudibranch transcriptomes, and data were extracted from the Agalma output using the matrix2genes function. The Agalma phylogeny pipeline selects genes that are phylogenetically informative and produce gene trees with no more than one terminal per taxa (i.e., paralogy pruning). The motivation for conducting these two searches was to target genes that would be informative for resolving both deep and shallow-scale evolutionary relationships. We calculated pairwise distances in Geneious v9.0 (Kearse et al., 2012) between the two most complete transcriptomes (C. westraliensis, G. coi) in order to evaluate the distribution of evolutionary rates between the results from the two Agalma searches since we predicted that these gene sets would resolve relationships at different evolutionary scales. Finally, a third set of target genes was taken from Teasdale et al. (2016) who identified a set of 500 genes from eupulmonate gastropods that were predicted to have broad phylogenetic utility.

| Identifying exon/intron boundaries
The Lottia gigantea gene IDs listed in Teasdale et al. (2016) were downloaded from Genbank and combined with the two sets of target genes identified with the Agalma v1.0 pipeline to generate a single set of genes. These targets were then blasted against the Aplysia F I G U R E 1 Gene search. (a) Workflow outlining steps for target gene search and identifying Chromodoris exon sequences. (b) Homologous genes retrieved from two Agalma searches in this study (Fast 5, All nudibranchs) and a published dataset in Teasdale et al. (2016). (c) Pairwise distance of gene alignments of Chromodoris westraliensis and Goniobranchus coi from two Agalma searches (All nudibranchs amino acid, Fast 5 nucleotide) and the Teasdale et al. (2016) dataset californica CDS genome to identify corresponding genes in Aplysia.
The A. californica genome was selected as a reference because it was the best-assembled genome in Heterobranchia, the subclass that Nudibranchia belongs in. We then used Exonerate v2.2.0 (Slater & Birney, 2005) to find introns by specifying the Aplysia genes as the query and the Aplysia genome (AplCal2.0 assembly Broad Institute) as the target. The output from Exonerate, a text file with concatenated exons per gene, was parsed into individual exons and exons >200 bp were saved using a custom python script. To find the corresponding exons in the Chromodoris transcriptomes, Exonerate was run again using these Aplysia exons as the query and the C. magnifica and C. westraliensis transcriptomes as targets. Exons were discarded from the Exonerate output if they were less than 115 bp in length, less than 65% similar to the Aplysia query exons, or longer than the original query exons. The single best sequence for each exon from either C. magnifica or C. westraliensis was selected based on highest score and all resulting exons were combined into a single file. These "best" exons were reciprocally blasted against each Chromodoris transcriptome to select genes that were variable (and thus informative) between the two Chromodoris species. A custom python script was then used to filter the results of the reciprocal blast based on percent identity. Only BLAST hits that were 92%-99% similar were retained, as hits that were 100% similar would not be informative within Chromodoris and hits <92% similar may have been paralogs.
If an exon blasted to any other exon in the target set, both were removed from the target set to prevent baits from hybridizing to multiple targets. Sequences of target exons from C. westraliensis were sent to Arbor Biosciences (formerly MYcroarray) in Ann Arbor, Michigan for bait design and synthesis.

| DNA extraction and targeted sequencing
The majority of DNA extractions for exon capture were taken from previous work by Layton et al. (2018) using samples of foot tissue from 65 Chromodoris specimens preserved in 96% ethanol, and an additional specimen preserved in 4% glutaraldehyde. Four additional DNA extractions were performed in this study for outgroup taxa (Ardeadoris egretta, D. atromarginata, G. coi, G. fidelis) using Qiagen DNeasy blood and tissue kits with an elution volume of 100-200 μl following manufacturer's protocols. Extractions were quantified on a Qubit, and in cases where concentrations were low, several extractions were prepared from the same animal and combined and concentrated using a Zymo DNA Clean and Concentrator kit following manufacturer's protocols. A total of 69 samples were sent to Arbor Biosciences for library preparation and target capture using the MYbaits-1-12 targeted sequencing kit protocols. DNA samples were sheared with sonication using a Qsonica Q800R instrument to an average insert length of 250bp, and eight libraries were pooled per reaction for a total of nine multiplex capture reactions. Sequencing was performed on a half lane of the Illumina HiSeq 2,500 platform with 100bp paired end reads also at Arbor Biosciences.
2.6 | Bioinformatic processing, phylogenetic analysis, and SNP calling Trimmomatic v0.36 (Bolger, Lohse, & Usadel, 2014) was used to remove adapter sequences, exon capture reads with a quality score below 15 in a 4-bp sliding window, and reads shorter than 26 bp. HybPiper v1.3.1 (Johnson et al., 2016) was employed to assemble the cleaned reads into contigs of the targeted regions of the genes. Briefly, the reads were mapped to a reference file of concatenated bait sequences using Burrows-Wheeler Aligner (BWA) (Li & Durbin, 2009). Mapped reads were then assembled de novo by gene using SPAdes v3.13.0 (Bankevich et al., 2012), and the resulting contigs were trimmed to include only the targeted exons with Exonerate v2.2.0. HybPiper produced an unaligned fasta file for each gene, containing a DNA sequence for each sample, and a series of summary statistics. For example, Hybpiper uses BWA to map the reads to contigs to present a value for percent reads on target.
Genes that did not enrich or enriched poorly (genes whose contigs were <50% of the reference) were removed (n = 125). The resulting gene files were aligned using MAFFT v7 (Katoh & Standley, 2013).
All gene alignments were trimmed with the "strict" method in tri- Phylogenetic analyses were conducted separately on two different datasets; (a) the full dataset containing all loci identified in this study, including loci identified by Teasdale et al. (2016), and (b) only those loci identified by Teasdale et al. (2016). We chose to assess the informativeness of just those loci identified by Teasdale et al. (2016) as these were predicted to have broad phylogenetic utility across Gastropoda. A mitochondrial (COI, 16S) dataset derived from Layton et al. (2018) was also used for comparing tree landscapes in the Treespace package in R (Jombart, Kendall, Almagro-Garcia, & Colijn, 2017). Individual gene alignments were concatenated for model testing (Kalyaanamoorthy, Minh, Wong, von Haeseler, & Jermiin, 2017) and ML analysis in IQtree v1.6.8 (Nguyen et al., 2015) with 1,000 ultrafast bootstrap replicates (Hoang, Chernomor, von Haeseler, Minh, & Vinh, 2018). Nodes with <50% bootstrap support were collapsed in the ML phylogeny. Additionally, individual gene ML trees were constructed with the same methods and then used as an input for summary coalescence analysis in ASTRAL II (Mirarab & Warnow, 2015) with 100 bootstrap replicates. The cophylo function in the Phytools package in R (Revell, 2012) was used to examine congruence in multiple sets of trees. Bootstrap trees and consensus trees for each exon dataset (i.e., full gene, Teasdale gene) and each tree-building method (i.e., concatenated ML, ASTRAL) were used to explore topological variability in the Treespace package in R using the Kendall Colijn distance metric to produce the PCA (Jombart et al., 2017). For the mtDNA dataset, bootstrap and consensus trees from an ML analysis in Layton et al. (2018) were explored in Treespace. The diversity of phylogenetic tools and genetic datasets employed in this study provided a unique opportunity for comparing topological variability among methodologies.
HybPiper produced a binary alignment map (BAM) file for each sample, which was used to call variants with the Genome Analysis Tool Kit v4 (GATK) (McKenna et al., 2010) following the GATK best practices for variant calling. PCR duplicates were removed from the BAM files which were then sorted and used to call SNP variants using GATK haplotype caller in gVCF mode. Variants were selected for analysis after a hard filtering step (QD < 2.0, MQ < 40.0, FS > 60.0, SOR > 3.0, MQRankSum < −12.5, ReadPosRankSum < −8.0) using the GATK SelectVariants tool. The resulting VCF file was pruned for linkage disequilibrium in PLINK 2.0 (Chang et al., 2015) by removing SNPs that were correlated (r 2 > 0.2) to any other SNP in a 50bp sliding window. The pruned VCF file was used with a custom python script to find alleles that were fixed differently in putative parental lineages for individuals where hybridization was suspected. We report this as the percentage of heterozygous fixed alleles. In F1 or F2 hybrids, we would expect the heterozygosity rate of alleles that are fixed differently between the two parental lineages to be very high.
To determine whether a signal of introgression was present among closely related species, we ran hybridization detection using phylogenetic invariants (HyDe) (Blischak, Chifman, Wolfe, & Kubtako, 2018) with the concatenated exon alignment and all individuals. We ran the full hybridization detection analysis using the run_hyde.py script to test all possible triplet combinations (Blischak et al., 2018).

| Searching for homologous loci in nudibranchs
Post hoc comparisons showed that RNA extractions significantly differed in concentration depending on extraction protocol (onetailed t-test, p = .01), with an average of 177.9 ng/μl using the Trizol protocol (N = 7) (range = 22.7-311.1) and 50.6 ng/μl using the Qiagen RNeasy protocol (N = 6) (range = 15.2-81.8). The total number of raw sequence reads ranged from 32.8 to 47.5 million per sample, with an average of 39.9 million, and the total number of transcripts ranged from 88,185 to 180,070, with an average of 137,528 (Table 1). The mean length of assembled transcripts ranged from 670 to 954 bp, with an average of 787 bp. The length of assembled transcripts significantly differed between samples that were extracted with the Trizol and Qiagen protocols (one-tailed t test, p = .04), with an average of 823 bp in the former and 757 bp in the latter. There was no significant correlation between number of reads and mean transcript length (p = .09). A maximum-likelihood (ML) analysis of all fifteen nudibranch transcriptomes is presented in Figure 2.  Figure 1b. After genes were filtered for phylogenetic informativeness, the number of genes and exons (>120 bp) targeted for sequence capture totaled 1,774 and 2,925, respectively. The targeted genes ranged in length from 120 to 4,163 bp. The average number of targeted exons per gene was 1.78. Pairwise distances of gene alignments from C. westraliensis and G. coi and each Agalma search are presented in Figure 1c.
The total volume of DNA for exon capture ranged from 10 to 180 µl, and the concentration ranged from 6.84 to 500 ng/µl. A total of three samples originated from pooled extractions; the original concentrations of the extractions ranged from 0.27 to 7.88 ng/µl prior to pooling. Capture efficiency varied between species, with the number of genes with contigs ranging from 1,061 (60%) in the out- Chromodoris aspersa (WAMS67676) had the lowest number of genes with contigs (1,456, 82.1%) within Chromodoris. The number of reads and genes mapped per sample is available in Table 2. Percentage of reads on target ranged from 41.3% in the outgroup G. fidelis to 78.3% in C. kuiteri (WAMS103139), with an average of 60% for all Chromodoris. A total of 144 genes were removed that did not enrich, or enriched poorly, likely due to intron sites that differed between The majority of genes sequenced from baits (N = 1,085) were derived from loci in the Fast 5 dataset.

| Exploring mitonuclear discordance in Chromodoris phylogeny
The full-gene ASTRAL phylogeny recovered most species-level entities and sister-species relationships with high support (Figure 3).
The number of species recovered in both the ASTRAL and concatenated ML analyses was largely congruent with a recent mtDNA analysis (Layton et al., 2018), with the exception of four individuals that showed mitonuclear discordance (Figure 3). In all four cases, the nuclear signal and color pattern were concordant, but the mitochondrial signal was that of another species, including two individuals that were previously identified as mimics in Layton et al. (2018).
This includes the "burni" morph of C. colemani (C. colemani × burni), which had the mitogenome of C. colemani but with the nuclear signal and color pattern of C. burni, and the "magnifica" morph of C. joshi (C. joshi × magnifica), which had the mitogenome of C. joshi but with the nuclear signal and color pattern of C. magnifica. Additionally, a specimen that was previously identified as C. magnifica (C. magnifica × meso) had the mitogenome of C. magnifica but with the nuclear signal and color pattern of a new mesophotic reef species (C. sp. meso), and a specimen that was previously identified as C. elisabethina had the mitogenome of C. elisabethina but with the nuclear

| Assessing introgression among closely related species
After pruning for linkage disequilibrium, 3,811 SNPs were retained from a total 66,247. For individuals with mitonuclear discordance, the heterozygosity of fixed alleles ranged from a minimum of 6.03%

| An informative bait set for Chromodoris nudibranchs
This study has established an informative bait set to improve phylo-  (Teasdale et al., 2016). It is unlikely that a nucleotide search returned more genes due to incorrect identification of paralogs as homologs because (a) a paralogy-pruning step in the Agalma Phylogeny pipeline ensures that mostly orthologous genes are included in subsequent matrix construction, and (b) the Hybpiper analysis did not produce any paralog warnings, meaning that there were no cases of multiple contigs assembled de novo from the capture data mapping to >85% of the same gene target.
We also included a set of genes from Teasdale et al. (2016) that were expected to have broad phylogenetic utility, and we confirmed that a subset of these genes (N = 149) were informative for resolving sister-species relationships in Chromodoris. Although we originally predicted that the different gene sets would be informative at different scales of divergence, there was no difference in the distribution of evolutionary rates between the sets of genes.

| Mitonuclear discordance and mimicry
This study employed phylogenomic data to test hypotheses of mimicry that derive from Layton et al. (2018). We uncovered mitonuclear discordance in two of these purported mimics, while in a third mimic, the mitochondrial and nuclear signals were concordant.
Concordance among mitochondrial and nuclear signals renders this individual a true mimic, while discordance among molecular signals suggests that the others derive from introgression or mitochondrial capture. Only a few cases of mitonuclear discordance have been reported in marine molluscs to date (i.e., mussels, Quesada, Wenne, & Skibinski, 1999;Rawson & Hilbish, 1998;octopus, Amor et al., 2019), and our study adds the first known cases in nudbranchs where the mitochondrial signal is incongruent with both the nuclear signal and color pattern. Mitonuclear discordance can arise from incomplete lineage sorting (ILS), but it is likely that discordance in this study has arisen through introgression or mitochondrial capture, although population-level sampling would help differentiate these patterns (Twyford & Ennos, 2012). First, Layton et al. (2018) demonstrate substantial mtDNA divergence among species that is largely supported with nuclear data in this study, suggesting that there is coalescence in these markers (Hinojosa et al., 2019). Secondly, these individuals exhibit morphologies that are strikingly similar to one of the putative parental lineages, which may support a scenario of asymmetric introgression (e.g., Wallace et al., 2011). Lastly, hybridization detection analysis indicates that two individuals with mitonuclear discordance are hybrids, although the parent species identified in this analysis do not match the expected parental lineages which might reflect ancestral hybridization or ILS. In any case, a lack of nuclear introgression in two other cases of discordance points toward a pattern of mitochondrial capture (Good, Vanderpool, Keeble, Bi, 2015). Introgression and mitochondrial capture are likely scenarios for this group given that they often occur sympatrically, are closely related, and have only recently radiated (e.g., Foltz, 1997;Mallet, Besansky, & Hahn, 2015). Although hybridization has been considered rare in marine invertebrates, previous work has shown evidence for introgression in several broadcast spawning invertebrates, including bivalves (e.g., Gardner, Skibinski, & Bajdik, 1993), ascidians (Nydam & Harrison, 2010;Nydam et al., 2017), and corals (Vollmer & Palumbi, 2002;Willis, van Oppen, Miller, Vollmer, & Ayre, 2006). Evidence for introgressive hybridization in Chromodoris is particularly intriguing given that they do not broadcast spawn but rather engage in copulation as simultaneous hermaphrodites. Note: Museum registration, species name, collection locality, DNA concentration, mean quality score, number of reads mapped, proportion of reads on target, and number of genes mapped. Data generated by the Hybpiper pipeline are highlighted in blue. Mimics appear in bold and individuals with mitonuclear discordance appear underlined, with names in brackets representing discordant molecular signals. Species names derived from Layton et al. (2018).

TA B L E 2 (Continued)
Rare instances of interspecific matings in Nudibranchia have been observed by divers, but because no study has uncovered hybrid individuals, it was generally assumed that postmating-prezygotic reproductive isolation would prevent fertilization. This study is the first to uncover evidence for introgression and mitochondrial capture in nudibranchs, demonstrating that species delimitation relying solely on mitochondrial data is problematic for this group.
Introgression has been crucial in the evolution of mimicry in other taxa (e.g., Pardo-Diaz et al., 2012;Heliconius Genome Consortium, 2012), and thus future work will investigate the role of hybridization and mimicry in facilitating speciation in nudibranchs.
Color patterns have been employed for species delimitation in many groups of nudibranchs, but mimicry complicates these patterns.
The nearly identical phenotypic match between mimetic C. colemani and the model C. westraliensis (differing only by a white line separating the black dorsum and orange margin in the former) in Western Australia is likely due to the lack of congenerics in the region, given that a mosaic color pattern would be expected in the presence of multiple co-existing species (Akcali, Pérez-Mendoza, Kikuchi, & Pfennig, 2019).
The close match between mimic and model also confers a fitness advantage because predators learn this pattern more quickly (Mallet & Singer, 1987), although it remains unknown whether this mimicry is Müllerian or Batesian in nature. Although mimicry has been documented in tropical fish, with color patterns driving speciation in some groups (e.g., Puebla, Bermingham, Guichard, & Whiteman, 2007), less is known about mimicry in nudibranchs and in molluscs in general (but see Cheney et al., 2016;Winters et al., 2017;Winters, White, et al., 2018;. A study by Padula et al. (2016) uncovered mimicry in another chromodorid genus (Felimida), but that study was based on mitochondrial data and further work is needed to determine whether this constitutes a "pure" case of mimicry or whether this pattern has arisen through introgression. Two mimics in this study were recovered as putative hybrids in a hybridization detection analysis, but additional genomic resources and population-level sampling would advance our understanding of the evolutionary processes governing mimicry in this genus.

F I G U R E 3
Mimicry and mitonuclear discordance in nudibranchs. ASTRAL summary coalescent tree of 1,630 genes with 100 bootstrap replicates. Bootstrap support is indicated at each node, with an asterisk marking nodes with 100% support. Triangles represent collapsed clades. Individuals with mitonuclear discordance are marked with a square, and mimics are marked with a circle (see inset). Names derived from Layton et al.(2018), and "txt" denotes a transcriptome sample. Inset: cases of mimicry and mitonuclear discordance from this study. Colored squares around each image match the legend provided, where identification derives from either mitochondrial signal, nuclear signal or color pattern. C. colemani (Photographer: Gary Cobb) and Cwestraliensis (Photographer: Bruce Potter) images derived from the Sea Slug Forum

| Topological incongruence and implications for Chromodoris phylogeny
This study reports incongruence between phylogenies generated with concatenated ML and ASTRAL summary coalescent methods, where the former produced phylogenies that were discordant between different datasets (i.e., sets of genes) and the latter produced highly concordant phylogenies consistent across different datasets.
The highly variable phylogenies produced by the concatenated ML analyses suggest that this method underperforms for taxonomic groups that have recently diverged and where gene trees may be in conflict. In fact, Kubatko and Degnan (2007) demonstrated that concatenating loci when gene trees are in conflict can produce incorrect phylogenetic estimates, and therefore the choice of phylogenetic method is important for resolving relationships in recently radiated taxa. Bragg et al. (2018) found that concatenated ML and ASTRAL phylogenies of a rapid radiation of Australian skinks were largely congruent, contradicting results in this study. Nonetheless, incorporating data from a large number of loci means that species trees will often converge on a single topology (e.g., Blom, Bragg, Potter, & Moritz, 2016).
Some sister-species relationships differ between the full-gene 2014). The evolution of spots could also be traced on a phylogeny of all chromodorid nudibranchs in order to determine the evolutionary status of this character in a broader context.
The exon capture approach employed in this study has significantly improved resolution of Chromodoris phylogeny and demonstrates, for the first time, mitonuclear discordance in nudibranchs likely originating from introgression or mitochondrial capture. This study also confirms a "pure" mimic in Western Australia, prompting additional investigation into the genomic basis of mimicry and warning colouration in these enigmatic animals.

ACK N OWLED G EM ENTS
We are grateful to our collaborators who have contributed speci-

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.