Pattern blending enriches the diversity of animal colorations

Comprehensive analyses of fish colorations illuminate a simple mechanism creating intricate patterns on animals.


INTRODUCTION
Animal color patterns play critical roles in animal survival, adaptation, intra-and interspecific communication, and speciation (1)(2)(3)(4)(5)(6). Traditionally, and reasonably, biologists have used animal coloration as an important key to identification of a species (7,8). We often see apparent discreteness or qualitative differences in color patterns even between closely related animals (9)(10)(11)(12)(13) and recognize these differences as strong evidence for novel/distinct species. At the same time, many color pattern motifs-from simple spot or stripe patterns to complex and camouflaged labyrinthine patternsappear repeatedly across a wide range of taxonomic groups (14)(15)(16). As a result, we see many cases in which a certain animal can more closely resemble distantly related species than its own close relatives (Fig. 1A). This extraordinary level of homoplasy-ubiquitous, long-range similarity over local divergence-is characteristic of animal coloration, which is rarely seen in other morphological traits, such as size and shape of body parts.
The evolutionary and developmental mechanisms underlying this peculiar, cross-taxonomic diversity of animal color patterns have remained largely unexplored until recently. In the past few years, however, progress has been made by identifying genes that are involved in similar skin patterns across multiple animal species (14)(15)(16). Some simple motifs, such as longitudinal stripes and wing spots, have been shown to be formed by using genetically encoded prepatterns originating from spatial landmarks or structural components of other body parts (12,17). Accordingly, with respect to these simple pattern motifs, intuitive explanations for their repeated occurrence across clades have been provided: Within each lineage, the painting of prepatterned regions is just switched on and off by particular genes (12,15). By contrast, more complex colorations such as camouflaged labyrinthine patterns have apparently no obvious relationship with any structures or body parts and cannot be formed simply from prepatterned positional information. Hence, it still remains a mystery how these intricate patterns have evolved and why unexpectedly similar, elaborate pattern motifs repeatedly and ubiquitously emerged across distantly related animal taxa.
A theoretical approach (18) has suggested an intriguing possibility to unravel this mystery. Using a specific class of the reactiondiffusion (RD) system (19) first proposed by Alan Turing (20), it has been predicted that complex and camouflaged markings can be formed by the morphological "blending" of simple patterns by hybridization (18). For example, crossing between virtual organisms with inverted spot patterns-one with light spots on a dark background and the other with dark spots on a light background-will necessarily result in intricate labyrinthine patterns as morphologically "blended" intermediate phenotypes ( Fig. 1B and fig. S1).
If this prediction is widely applicable to real animals, one can expect that organisms having complex and camouflaged motifs can abruptly emerge from events of hybridization between species that have simple patterns. This "pattern blending" hypothesis would account for the enigmatic repeated occurrence of similar, intricate pattern motifs across distantly related taxa, while simultaneously offering a plausible explanation for the discrete, qualitative pattern differences between closely related species. Although the model prediction has been partially confirmed by empirical evidence such as artificial hybridization of salmonid species (18), natural instances supporting the pattern blending hypothesis so far are very limited, and it is unclear whether, or to what extent, pattern blending contributes to the actual diversity of animal color patterns.
In this study, I test the pattern blending hypothesis by exploring the natural diversity of color patterns in fish lineages. First, I analyze the skin patterns of more than 18,000 species of marine and freshwater fishes and examine phylogenetic collocations of pattern motifs to infer the mechanistic associations among them. I show that apparent phenotypic similarity does not necessarily correspond to the mechanistic association between pattern motifs. On the contrary, strong positive associations between phenotypically dissimilar pattern motifs are present, which are consistent with the pattern blending hypothesis. Next, I provide demonstrative examples of pattern blending by using color pattern quantification and comparative genomic analyses targeting representative fish species. Last, by applying phylogenetic comparative methods, I show that pattern blending has contributed to the diversity of color patterns in multiple major fish lineages.

Diversity of color patterns in marine and freshwater fish
To examine the diversity of color patterns in fish lineages, I performed pattern annotation on images of 18,114 fish species referenced from public databases (FishBase: https://www.fishbase.org/ and FishPix: http://fishpix.kahaku.go.jp/fishimage-e/). I used 11 wellrecognized classes of pattern motifs for annotation (7,8): 4 types of stripe patterns (St-H, horizontal stripes; St-D, diagonal stripes/bands; St-V, vertical stripes/bars; Maze, labyrinthine stripes), 3 types of spot patterns (Sp-D, small dark spots on a light background; Sp-L, small light spots on a dark background; Eyes, eyespots/spots larger than eyes), and 4 other patterns (Area, area fill pattern; Sddl, saddle-like pattern; Bltc, blotch/speckle pattern; Mono, monotone/patternless). Each image was binary labeled for each class of pattern motifs using a web-based in-house annotation system, and the presence/absence of each pattern motif for each species and taxonomic group was recorded ( Fig. 2A, fig. S2, and data file S1). The compiled data revealed that all of the 11 color pattern motifs were widespread across the lineages, illustrating the ubiquitous long-range similarity over local divergence, although the range and abundance varied among the pattern motifs ( Fig. 2A).
Using this dataset, I examined the mechanistic associations among color pattern motifs by evaluating their phylogenetic co-occurrences. It is reasonable to assume that mechanistically related motifs, wherein pattern transition can readily occur, are likely to co-occur within each small taxonomic group, such as genus. Then, the mechanistic association between pattern motifs can be inferred based on this "phylogenetic collocation" (i.e., within-genus co-occurrence of pattern motifs). This association should not necessarily or simply be a measure of superficial, phenotypic similarity between pattern motifs; rather, it should reflect the underlying mechanistic connections related to evolutionary, developmental, or ecological constraints and/or structures. I used three different association measures to evaluate the within-genus co-occurrence between each pairwise combination of the 11 classes of color pattern motifs (Fig. 2, B and C). The Jaccard index (JI) has been classically used for evaluating similarities between two sets of binary data in many fields, including community ecology (21). The log-likelihood (LL) and z score (Z) have been frequently used for quantifying collocational strength in the field of corpus linguistics (22). According to these indices, strong associations were present among spot-like motifs (Sp-D, Sp-L, and Bltc; JI = 0.31 to 0.40, LL = 144 to 203, Z = 6.9 to 9.7), which are concordant with their phenotypic similarity. In contrast, despite the obvious morphological similarity among stripe motifs (St-H, St-D, St-V, and Maze), their associations were relatively weak, if any (JI = 0.13 to 0.24, LL = 9.8 to 88, Z = 2.4 to 7.7). This implies that developmental and/or ecological constraints may exist within each class of stripe motifs, so that the transition to another stripe class is restricted. For example, the development of horizontal stripes in fish may depend heavily on prepatterns related to structural body parts (9), as is the case in mammals and birds (14,17).
The most notable associations were found in distinctly different motif pairs: labyrinthine and spot motifs (Maze and Sp-D: JI = 0.29, LL = 269, Z = 11.8; Maze and Sp-L: JI = 0.28, LL = 207, Z = 12.1) (Fig. 2, B and C). Considering the apparent phenotypic differences between labyrinths and spots, their strong association seems counterintuitive in terms of morphology. On the other hand, it is consistent with the pattern blending hypothesis and in silico hybridization experiments ( Fig. 1B and fig. S1). This consistency is more evident in triple co-occurrence analyses of pattern motifs: When focusing on the genera in which dark (Sp-D) and light (Sp-L) spots co-occur, the pattern motif most strongly and substantially associated with these two spot motifs was indicated to be the labyrinthine motif (Maze) by all three association measures (Fig. 2D). To confirm these results, I used three other measures [Sørensen-Dice coefficient (SDC), Simpson similarity index (SSI), and T score (T)] and obtained qualitatively similar results ( fig. S3). By contrast, for motif pairs other than dark and light spots, the strongest association was detected with motifs other than the labyrinthine patterns ( fig. S4).

Pattern quantification and genomic analyses of pufferfish
To investigate the possibility of pattern blending in more detail, I performed an in-depth analysis of the skin pattern diversity in a representative small taxonomic group: pufferfish species of the genus Arothron, which are known for the compact size of their genome (less than 400 Mb, smallest among vertebrates), as well as their remarkable and diverse color patterns (Fig. 3) (23,24). Pattern quantification targeting a total of 120 individuals from nine species revealed that pattern variations observed among Arothron fishes were in good agreement with numerical simulations based on the RD model (Fig. 4A, figs. S5 and S6, and data file S1). Above all, the labyrinthine patterns of Arothron species, Arothron multilineatus Matsuura, 2016, Arothron carduus (Cantor, 1849), and Arothron mappa (Lesson, 1831) (Fig. 3, D to F) were shown to have "intermediate" color tones between two types of spot patterns (a darker background tone for light spot patterns and a lighter background tone for dark spot patterns) (Fig. 4B) and a "transgressive" property as to the pattern complexity (labyrinthine patterns being more complex than both of spot patterns) (Fig. 4C).
This twofold character of the labyrinthine patterns of Arothron species shows a remarkable correspondence with that of the blended patterns emerging in the in silico hybrids between light spots and dark spots ( Fig. 3A and fig. S1). According to the pattern blending  To test this prediction, I performed comparative genomic analyses of Arothron species. For 12 species of the genus Arothron, including 3 labyrinthine species, whole-genome libraries were constructed and sequenced using freshly collected samples and/or museum specimens. I generated a de novo genome assembly of one Arothron species, A. firmamentum (aroFir_1.0: the draft genome size is 335 Mb, N50 contig and scaffold sizes are 15.9 and 135 kb, respectively), and mapped the reads of 20 individuals from 12 Arothron species to this reference genome (mapping rates of all samples ranged from 96 to 98%, except for museum samples). Twenty-eight million single-nucleotide polymorphisms (SNPs) were identified across the Arothron species in the form of genotype likelihoods, and complete or nearly complete mitochondrial DNA (mtDNA) genome sequences were also reconstructed from the reads for each species.
According to the phylogenetic tree based on mtDNA genomes (Fig. 5A), the labyrinthine species, A. carduus, is very closely related to the white-spotted species, A. reticularis, while another labyrinthine species, A. multilineatus, appears almost identical to the black-spotted species, A. stellatus. On the other hand, whole-genome analyses indicated that these labyrinthine species have mixed ancestry: Principal components analysis (PCA) based on whole-genome SNPs showed that they were both located near the midpoint of the clusters of the white-spotted A. reticularis and the black-spotted A. stellatus in the plot of the first two components (Fig. 5B). Furthermore, from the admixture analysis based on the same SNP dataset, the ancestry proportion of A. multilineatus was estimated to be approximately half of both A. stellatus and A. reticularis (Fig. 5C). A. carduus also showed mixed ancestry between A. stellatus and A. reticularis, with a small contribution from A. firmamentum (Fig. 5C). I also estimated interclass heterozygosity (25) based on 5 million diagnostic sites that segregate A. stellatus and A. reticularis. Most of the callable diagnostic sites were heterozygous in the labyrinthine individuals (A. multilineatus: 97 and 99%, A. carduus: 87%), indicating that these species are in an early generation of admixture.
These data strongly support the prediction from the pattern blending hypothesis: Two labyrinthine Arothron species, A. multilineatus and A. carduus, are both derived from the interspecific hybridization of the white-spotted A. reticularis and the black-spotted A. stellatus, although the parental combinations and color patterns were inverted (Fig. 3, D and E; white lines on a black background and black lines on a white background, respectively). According to the original description (24,26), no morphological traits that can distinguish these two species from other Arothron species have been recognized besides the peculiar skin patterns. Their reticulated colorations may, thus, indicate their reticulated origins, rather than their species identities.
As for the other labyrinthine species, A. mappa, it seems unlikely that it was in the early generations of admixture (Fig. 5, A to C). However, the ABBA-BABA test (27) [Patterson's D statistic and its derivative, D min (28)] detected traces of past hybridization or introgression between multiple combinations of this and other Arothron species, including black-spotted and white-spotted relatives [78% of all possible trios (28 of 36) have a significantly positive D min score (Holm-Bonferroni familywise error rate, <0.01); table S1], suggesting the probable occurrence of pattern blending events during the evolution of this Arothron species.

Phylogenetic comparative analysis of color pattern motifs in fish lineages
If pattern blending events have contributed substantially to the evolution and diversity of animal skin patterns, their traces could be detected as signals of correlated evolution between distinct pattern motifs. For example, under the condition of frequent pattern blending, a labyrinthine motif is more likely to emerge in taxonomic groups that contain both dark and light spotted species than in groups without one or both of these spot motifs.
To examine the possibility of such correlated evolution between labyrinthine and spot motifs in fish lineages, I applied phylogenetic comparative methods (29) to the genus-level trait data for each of the top 10 major fish orders: Cypriniformes, Siluriformes, Perciformes, Cichliformes, Characiformes, Gobiiformes, Cyprinodontiformes, Blenniiformes, Pleuronectiformes, and Tetraodontiformes (ranked based on the number of genera within each order). I tested independent and dependent models of pattern motif evolution ( fig. S7): Under the independent model, all three motifs (Sp-D, Sp-L, and Maze) evolve independently without affecting one another, while under the dependent model, the evolution of one motif can be affected by the state of the other two motifs. Bayesian analyses (30) found substantial support for the dependent model (correlated evolution of labyrinthine and spot motifs) in 8 of 10 orders (log Bayes factor = 4.18 to 43.1; table S2). Furthermore, the estimated rates of transition from "without Maze" to "with Maze" (Maze gain rate, a to d in Fig. 6A) were higher when both dark and light spot motifs were included in the pattern portfolio (rate d), than when one or both of these spot motifs happened to be absent (rates a to c) (Fig. 6B). By contrast, transition rates from "with Maze" to "without Maze" (Maze loss rate, e to h in Fig. 6A) tended to be lower when both spot motifs coexisted (Fig. 6C). This means that the coexistence of dark and light spot motifs increases the probability of the emergence of labyrinthine motif and reduces the probability of its loss. These tendencies were also confirmed in further analyses, including comparisons with other motif combinations performed on the entire fish lineage (figs. S8 and S9). These findings suggest the frequent occurrence of pattern blending events in multiple major fish lineages and its contribution to the enrichment of color pattern diversity.

DISCUSSION
In this study, I presented empirical evidence supporting the pattern blending hypothesis. Comprehensive analysis of pattern diversity among fish skins revealed strong mechanistic links for the qualitatively distinct pattern motifs, which are consistent with the hypothesis. Furthermore, the pufferfish examples clearly demonstrated that pat-tern blending actually occurs in the wild: Hybridization between species with simple spot patterns can create camouflaged, labyrinthine organisms that appear like distinct or novel species, just as the model predicted.
In nature, there are many species that have labyrinthine motifs. I found more than 900 fish species with labyrinthine motifs from >160 families, accounting for 5% of all fish species and 30% of all fish families examined here ( Fig. 2A). It is reasonable to expect that at least a part of these labyrinthine animals originated from interspecific hybridization between spotted species, as is the case for Arothron pufferfishes. I anticipate that some of them may be just hybrids and have been deceiving taxonomists with their camouflaged coloration to be given unworthy taxonomic positions as novel/distinct species.
Another evolutionarily more intriguing possibility could be hybrid speciation (31). In the case of Arothron pufferfish, many traces of interspecific hybridization/introgression were detected between labyrinthine species and their spotted relatives. Moreover, in multiple major fish lineages, I found strong support for the correlated evolution of labyrinthine and spot motifs: Coexistence of dark and light spot species substantially increased the probability of the emergence of a labyrinthine motif, suggesting the frequent occurrence of pattern blending by hybridization and consequential hybrid speciation.
Accumulating evidence now indicates that hybridization is a potential driver of biological diversity (13,28,(31)(32)(33)(34)(35). Recently, it was shown that interspecific hybridization can form a new species very rapidly, even within a few generations (36). Considering that animal markings are frequently used as visual cues for the recognition of a species and mate choice (5, 6), qualitative and drastic changes in color patterns from parent species may easily result in reproductive isolation under conditions of assortative mating (37,38). In addition, a conspicuous yet disguising appearance can have multiple advantages in certain circumstances: camouflage from predators/prey by means of background matching or disruptive coloration, and aposematic signals if present with strong defensive tools, such as poisons or spines (1-3). Thus, both of the two major requirements for hybrid speciation, namely, reproductive isolation from the parent species and exploitation of new niches (31), can be fulfilled with regard to the camouflaged labyrinthine color patterns.
Although I do not rule out any other evolutionary, developmental, and ecological factors that may have caused camouflaged colorations to arise, it could be evolutionary feasible that a group of hybrid individuals with "intermediate and transgressive" phenotypes, such as labyrinthine patterns, evolved as a distinct lineage of hybrid origin to form a new species [speciation by fusion (38)]. I expect that such reticulate evolution of color patterns may account for a key part of the general pattern diversity-ubiquitous, long-range similarity over local divergence-observed among animals. Further exploration of diverse biological groups, including rare species in museum collections, should shed light on the rich structure of the phenotypic space for this seemingly two-dimensional trait, and possibly uncover a secret path through the evolutionary labyrinth of animal color patterns.

Mathematical modeling and simulation
Mathematical modeling and numerical simulation were performed as previously described (18). Briefly, I used a class of RD system described as  ∂ v where u and v are the concentrations of hypothetical factors, f and g are the reaction kinetics, and D u and D v are the hypothetical diffusion coefficients (or their mathematical equivalents) for u and v, respectively. The reaction rate, R, was introduced for the convenience of parameter adjustment. The reaction kinetics and parameters are as follows Python library were used for rendering. Note that similar pattern transitions to those observed in these simulations can be reproduced using a variety of other models (39,40).

Annotation of fish color patterns
Fish images with identified species names were referenced from public databases (FishBase: https://www.fishbase.org/ and FishPix: http://fishpix.kahaku.go.jp/fishimage-e/). A total of 19,755 images of 18,114 species from 4069 genera and 559 families were subjected to pattern annotation using a web-based in-house annotation system. Based on the descriptions of pattern motifs commonly used in original papers, pictorial books, and other publications on this topic (7,8), I formulated the following 11 classes of pattern motifs: 4 types of stripe patterns (St-H, horizontal or longitudinal stripes running parallel to the body axis; St-D, diagonal stripes/bands running at an angle of roughly 20° to 70° to the body axis; St-V, vertical stripes/ bars oriented perpendicular to the body axis; and Maze, labyrinthine stripes with less or no directional anisotropy), 3 types of spot patterns (Sp-D, small dark spots that are evenly distributed on a light background; Sp-L, small light spots on a dark background; and Eyes, ocellated spots or one or a few spots larger than the actual eyes), and 4 other patterns (Area, area fill pattern in which large areas are painted in a single color; Sddl, saddle-like markings that line up on the back of the body; Bltc, blotch/speckle pattern with irregularly shaped and not evenly spaced markings; and Mono, monotone/ patternless). Pattern annotation was carried out by multilabel classification, that is, each image was binary labeled for each class of pattern motifs (1, containing the pattern motif; 0, not containing the pattern motif).

Co-occurrence analysis of color pattern variations
For each genus and each class of pattern motifs, the occurrence of species having that pattern motif was counted and binary labeled (1, containing species with the pattern motif; 0, not containing species with the pattern motif) (fig. S2A). Next, within-genus pattern co-occurrence between each pairwise combination of pattern classes was examined and counted ( fig. S2B). The significance and strength of the pattern motif co-occurrences were evaluated using six association measures (21,22): JI, SDC, SSI, LL, T score (T), and z score (Z). The association measures for the motif A and B pair are defined as where O ij is the observed frequencies, i and j denote occurrence (1) or nonoccurrence (0) of motifs A and B, respectively (e.g., O 11 is the number of genera where motifs A and B co-occur, whereas O 10 is the number of genera where motif A occurs but motif B does not), A i and B j are the total numbers of occurrence (i, j = 1) or nonoccurrence (i, j = 0) of motifs A and B, respectively, and E ij = A i B j / ij O ij is the expected frequencies. A total of 2384 genera containing two or more species was used for the calculation. For the calculation of association measures for triple co-occurrence, a total of 1651 genera containing three or more species was used, and a combination of two motifs [e.g., dark (Sp-D) and light (Sp-L) spots] was treated as motif A in the above definition (i.e., the occurrence of motif A means the co-occurrence of Sp-D and Sp-L; fig. S2C).

Quantitative analysis of color patterns
Various methods and metrics have been so far developed and proposed for the quantification of color patterns (18,41,42). For a detailed quantitative analysis of the body patterns of Arothron species, I focused on two indicators: pattern complexity (based on elementwise circularity/elongation) and overall color tone (pattern lightness measured as the proportion of unpigmented area in a binarized image) (18). The pattern complexity score (PCS) is defined based on the area-weighted mean isoperimetric quotient of the contours extracted from each image where Q i = 4S i /L i 2 is the isoperimetric quotient (or circularity) of each contour, w i = S i / i S i is the area weight, and S i and L i are the area and perimeter of each contour, respectively. Here, pattern complexity was calculated on the basis of the area and contour length for each pattern element that makes up the overall pattern. The area of a shape enclosed by a contour is minimized when the shape takes its simplest form, a circle. As the shape becomes more complex, the contour becomes longer. Therefore, as a measure of the complexity of a shape, we can use the "circularity" of the shape, which is calculated by comparing the area of the shape enclosed by a contour of a certain length with the area of a circle with a circumference of the same length. The circularity takes a maximum value of 1 for a circle (the simplest shape) and decreases as the shape deviates from the circle. The complexity of the entire pattern can be calculated from this circularity obtained for each element. However, using a simple average of the element-wise circularity may be inappropriate because the overall impression of the complexity of a pattern depends on the pattern elements that occupy a large area of the entire pattern. For example, when a region of interest contains one complex, dominantly large element and a large number of simple, small elements, the pattern complexity calculated from the simple average of the element-wise circularity may underestimate the overall complexity of that region. Therefore, in this study, I used an area-weighted average of circularity and calculated the PCS by subtracting this value from 1. These quantification analyses were performed on binarized images using the OpenCV library and in-house Python script.
The PCS used in this study is considered appropriate as a quantitative index that can reliably capture the characteristics of patterns such as spots and maze and the transitions between them (Fig. 4); however, it is difficult to measure the differences in complexity between other combinations of patterns such as horizontal or vertical stripes and maze. For the latter, it may be useful to use other indices such as the aspect ratio and orientation of each pattern element (41) and/or their variances.
Sampling, DNA extraction, library preparation, and sequencing A total of 21 individual/tissue specimens from 12 Arothron species were obtained from various sources, including museum/aquarium collections and local aquarium shops. Identification of the species and nomenclature followed Matsuura (23,24). Total DNA was extracted from tissue samples of pectoral fin and/or skeletal muscle using the DNeasy Blood & Tissue Kit (Qiagen, 69504). Doublestranded genomic DNA libraries for high-throughput sequencing were prepared either with the TruSeq DNA PCR-Free/Nano DNA Library Prep Kit (Illumina, FC-121-3001/4001) or the NEBNext Ultra II DNA Library Prep Kit (New England Biolabs, E7645S).
Because the two labyrinthine species, A. carduus and A. multilineatus, are extremely rare, I used formalin-fixed specimens deposited in museum collections (HUMZ 35438 and KAUM-I. 13606, collected in 1973 and 2009, respectively). For the extraction of very short, degraded DNA fragments [mode length, 30 to 40 base pairs (bp)] from these museum specimens, the DNAs-ici!-F (Rizo, DS-0005) was used. In addition to the double-stranded DNA libraries, single-stranded DNA libraries were constructed for museum samples, based on the methods developed for ancient DNA analysis (43) with the following modifications: (i) in the extension step (steps 12 to 14), the extension primer CL9 was replaced with a newly designed extension primer for P5 (EPP5: 5′-ACACTCTTTCCCTACACGACGCTCTTCCGATCT-3′), so that widely used universal primers could be used for easy outsourcing of DNA sequencing, (ii) the concentration of 2′-deoxyadenosine 5′-triphosphate (dATP) was increased by fivefold to promote the generation of 3′ adenine overhangs, so that (iii) the blunt-end repair steps (steps 16 to 19) could be omitted, and the second adaptor ligation steps (steps 20 to 22) could be performed using the Ligation Module in NEBNext Ultra II DNA Library Prep Kit (New England Biolabs) and the NEBNext Adaptors. To avoid potential contamination, DNA extraction and library preparation for formalin-fixed museum specimens were conducted in an isolated room with a clean bench.
Whole-genome sequencing was performed at Macrogen using the Illumina HiSeq X platform. All animal experiments were conducted in accordance with the institutional guidelines of Osaka University.
Inferring population structure, admixture proportions, and introgression PCA was performed on the output file from ANGSD in Beagle format using PCAngsd v0.95 (57), which can handle genotype likelihoods as input data for improved accuracy in samples with low and variable sequencing depth. An in-house Python script was used for plotting. Genome-wide admixture proportions for each individual were estimated using NGSadmix v32 (58), with K = 4 to 6 ancestral populations based on the same dataset.
Interclass heterozygosity (the proportion of loci with alleles from both ancestral populations) of possible hybrid individuals were estimated using ANGSD and an in-house script based on a total of 5,198,905 diagnostic biallelic sites that segregate five individuals each of A. reticularis and A. stellatus.
ABBA-BABA analysis was performed using the D. nigroviridis genome (tetNig2) as an outgroup and reference; reads from Arothron species were remapped to the tetNig2 reference genome using BWA-MEM or BWA-backtrack (options are the same as those used for the aroFir_1.0 reference), and then processed with ANGSD (angsd -doAbbababa 1 -doCounts 1 -blockSize 5000000 -anc tetNig2.fa). Patterson's D statistic (27) for each combination of Arothron species was calculated, and the significance was evaluated with z scores obtained by the block-jackknifing using 5-Mb blocks. The D min statistic (28)

Analyses of correlated evolution of pattern motifs
For the phylogenetic comparative analyses, I used the time-calibrated phylogeny of ray-finned fishes (59) obtained from https://fishtreeoflife. org by using the R package "fishtree" (60). I pruned the tree by randomly selecting a single representative species per each genus for which genus-level pattern annotation data were available and prepared genus-level tree samples for each order by repeating this process 500 times per order.
Analyses of correlated evolution among pattern motifs were performed using BayesTraits V3.0.1 (30). I defined the portfolio of pattern motifs in a genus as the trait of that genus. The combinations of the presence/absence of each of the three pattern motifs (Sp-D, Sp-L, and Maze) were coded into eight different states (i.e., 000, 100, …, 111), and the Multistate module with reversible-jump Markov chain Monte Carlo (MCMC) in BayesTraits was used to estimate the evolutionary transition rates among these states for each order. I tested independent (IND) and dependent (DEP) models of pattern motif evolution ( fig. S7): Under the independent model, all three motifs evolve independently without affecting one another, while under the dependent model, the evolution of one motif can be affected by the state of the other two motifs. I used a uniform hyperprior ranging from 0 to 10 to seed the mean of the exponential prior. Each MCMC analysis was run three times for 600 million iterations with a burn-in of 100 million and sampling every 50,000 iterations. The marginal likelihood was estimated by the stepping stone sampler with 100 stones and 10,000 iterations per stone. Log Bayes factors were calculated as twice the difference between the log marginal likelihood of the dependent and independent models. I took log Bayes factors >2 as positive evidence, >5 as strong evidence, and >10 as very strong evidence (30).
Further analyses were performed on the entire fish lineage (2837 genera for which trait data were available) and various sets of three motif combinations. I tested partially independent models (indX, indY, and indZ), in which two motifs are correlated (evolve dependently), while the evolution of the other motif is independent (figs. S8A and S9A), as well as the independent and dependent models described above, for the following sets of three motif combinations: (i) dark spots, light spots, and one of the other nine motifs  S9, B and C). For these analyses, I used the same hyperprior as above and performed three MCMC runs each for 6 million iterations with a burn-in of 1 million and sampling every 5000 iterations. The marginal likelihood was estimated by the stepping stone sampler with 100 stones and 1000 iterations per stone. Log Bayes factors were calculated as twice the difference between the log marginal likelihood of each model.