Gudgeon fish with and without genetically determined countershading coexist in heterogeneous littoral environments of an ancient lake

Abstract Countershading, characterized by a darker dorsal surface and lighter ventral surface, is common among many animals. This dorsoventral pigment polarity is often thought to be adaptive coloration for camouflage. By contrast, noncountershaded (melanistic) morphs often occur within a species due to genetic color polymorphism in terrestrial animals. However, the polymorphism with either countershaded or melanistic morphs is poorly known in wild aquatic animals. This study explored the genetic nature of diverged color morphs of a lineage of gudgeon fish (genus Sarcocheilichthys) in the ancient Lake Biwa and propose this system as a novel model for testing hypotheses of functional aspects of countershading and its loss in aquatic environments. This system harbors two color morphs that have been treated taxonomically as separate species; Sarcocheilichthys variegatus microoculus which occurs throughout the littoral zone and Sarcocheilichthys biwaensis which occurs in and around rocky areas. First, we confirmed that the divergence of dorsoventral color patterns between the two morphs is under strict genetic control at the levels of chromatophore distribution and melanin‐related gene expression under common garden rearing. The former morph displayed sharp countershading coloration, whereas the latter morph exhibited a strong tendency toward its loss. The crossing results indicated that this divergence was likely controlled by a single locus in a two‐allele Mendelian inheritance pattern. Furthermore, our population genomic and genome‐wide association study analyses detected no genome‐wide divergence between the two morphs, except for one region near a locus that may be associated with the color divergence. Thus, these morphs are either in a state of intraspecific color polymorphism or two incipient species. Evolutionary forces underlying this polymorphism appear to be associated with heterogeneous littoral environments in this lake. Future ecological genomic research will provide insight into adaptive functions of this widespread coloration, including the eco‐evolutionary drivers of its loss, in the aquatic world.


| INTRODUC TI ON
The adaptive functions of animal coloration are a widespread and fascinating topic in ecology and evolutionary biology (Cuthill et al., 2017).
Countershading, characterized by a light ventral surface and darker dorsal surface, is among the most common color patterns in terrestrial and aquatic animals including mammals, birds, reptiles, fishes, and insects (Kiltie, 1988;Thayer, 1896). This dorsoventral pigment polarity is ecologically explained as an adaptation for camouflage through the self-shadow concealment and/or background matching, although it may have other benefits such as enhancing thermoregulation or protection against ultraviolet light (reviewed by Rowland, 2009;Ruxton et al., 2004). By contrast, noncountershaded (uniformly dark-colored) phenotypes, which are associated with increased melanin pigmentation on the ventral surface, known as melanism, are also widespread among terrestrial animals (Allen et al., 2012;Brown et al., 2017;Caro, 2005). Both genetically determined countershading and melanistic morphs can occur within a single population, and morph frequency variation has been associated with environmental heterogeneity (e.g., Nachman, 2005;Potash et al., 2020). Such color polymorphisms often have a simpler genetic basis than other ecologically relevant traits (McKinnon & Pierotti, 2010;Orteu & Jiggins, 2020). Therefore, they provide excellent opportunities for understanding the adaptive significance of animal coloration and elucidating the range of evolutionary processes that contribute to the maintaining genetic polymorphisms within and between populations (Gray & McKinnon, 2007;White & Kemp, 2016).
Many aquatic organisms are also countershaded, and the concealing function of countershading may be more beneficial underwater than in terrestrial habitats (Kiltie, 1988;Rowland, 2009;Ruxton et al., 2004).
Thus, countershading is thought to be a more effective antipredator adaptation among aquatic prey species, where predators can be attacked from all directions, because dark dorsal pigmentation provides a good match to the dark background when the prey is viewed from above, whereas a light ventral surface would match the downwelling light from overhead when it is viewed from below. Partly for this reason, the loss of countershading is far less common in aquatic habitats than in terrestrial systems (Venables et al., 2019). However, the hypotheses about the ecological and evolutionary mechanisms of using countershading to attain crypsis remain poorly tested in wild aquatic organisms (Kelley & Merilaita, 2015;Kelley et al., 2017). Although some examples of noncountershaded coloration have been reported, the ecological drivers and evolutionary genetic basis underlying dark pigmentation remain unknown (Caro et al., 2011;Lindgren et al., 2014;Venables et al., 2019).
A productive way to address these questions including the functional aspects of such color patterns would be to examine aquatic animals that exhibit stable color polymorphism with either countershaded or melanistic morphs and are amenable to ecological and genetic analyses.
As a potential candidate for this, we focused on a lineage of the gudgeon fish (genus Sarcocheilichthys, belonging to the subfamily Gobioninae, the family Cyprinidae), which inhabit the littoral zones of Lake Biwa, an ancient temperate lake in central Japan. Lake Biwa (surface area 670 km 2 ; maximum depth 104 m) is the oldest lake in East Asia dating to about 4 million years ago (Kawabe, 1994;Yokoyama, 1984), and its littoral zone comprises complicated bottom environments including sandy, pebbly, and rocky areas (see Komiya et al., 2011). The lake harbors 16 endemic fish species/subspecies (hereafter species) including two Sarcocheilichthys species (reviewed by Okuda et al., 2014). These two species have divergent body colors which are not sex-linked and distributional patterns (Hosoya, 1982;Nakamura, 1969). Sarcocheilichthys variegatus microoculus has a yellow-grayish body and occurs throughout the littoral zone and in rivers flowing into Lake Biwa, whereas Sarcocheilichthys biwaensis has a brownish body with yellowbrownish and brownish-black coloration, and occurs strictly in and around rocky areas. Thus, the two species co-occur in rocky areas, showing spatial variation in their frequency, and S. biwaensis could be found up to about 30% in these areas (Nakamura, 1969). At least visually, the former species display typical countershading and the latter species seems to exhibit dark pigmentation over its entire body including the ventral region and lack countershading pattern (see Figure 1). Although these fishes have been treated taxonomically as separate species, previous phylogeographic and population genetic analyses using neutral DNA markers have suggested no genetic differentiation between these species and the panmictic status of Sarcocheilichthys in Lake Biwa (Komiya et al., 2011(Komiya et al., , 2014Tabata et al., 2016). Therefore, we propose this Sarcocheilichthys system as a model of stable color polymorphism, with countershaded and melanistic morphs. However, to date, no study has examined the color pattern divergence or genetic aspects of color polymorphism in this system in detail, except for Nakamura (1969)'s brief description (data not shown).
Here we assessed the following questions in order to establish the Sarcocheilichthys system as a novel model for studying the evolutionary mechanisms of countershading color polymorphism in the aquatic world. First, we examined whether the divergence of dorsoventral color patterns in this system is under strict genetic control at the levels of chromatophore distribution and melanin-related gene expression. Second, we further examined whether this phenomenon is controlled by a few genes exhibiting Mendelian segregation, as has been shown in other taxa. Last, we applied restriction siteassociated DNA sequencing (RAD-seq) to explore the nature of genome-wide pattern of genetic divergence between these color morphs (i.e., taxonomic species) in an area where they coexist, to obtain evidence of color polymorphism within a single interbreeding population for this system.

K E Y W O R D S
genetic color polymorphism, incipient species, Lake Biwa, melanism, spatially heterogeneous environments 2 | MATERIAL S AND ME THODS

| Crossing experiment and common garden rearing
As mentioned above, regardless of their taxonomic status, the Sarcocheilichthys system in Lake Biwa appears to represent color polymorphism within a single population or a pair of incipient species with weak state of premating isolation (see also Discussion).
Therefore, for convenience, we will refer to S. variegatus microoculus as the S V morph and S. biwaensis as the S B morph. Intra-and intermorph crossing of the Sarcocheilichthys system can be achieved through artificial insemination using mature wild individuals (Nakamura, 1969). Therefore, to create experimental families for common garden rearing, backcrossing, and F 2 hybrid populations, we collected parental fish (three males and three females) of the S V morph from a tributary of the Ane River, which flows into Lake Biwa, and parental fish (three males and two females) of the S B morph around off Onoe in Lake Biwa. These sampling locations are indicated in Komiya et al. (2014) as locality codes 9 and 7, respectively.
The parental individuals were transferred to 180-L tanks, where they were maintained until artificial insemination.
We created two full-sibling families of the S V morph (F V1 and F V2 ), two full-sibling families of the S B morph (F B1 and F B2 ), and a hybrid family (F H1 ) between the two morphs. We also created a mixed family (F BM ) from fertilized eggs of the S B morph from captive broodstock maintained at the Lake Biwa Museum (Japan). Artificial insemination was performed from May to June 2015 using eggs from each ovulated female and sperm from each male, which were obtained by pushing the belly. An F 1 intercross family (F H1 ) was created using a male S V morph that was not used for creation of full-sibling families and a female S B morph that was also used to create F B2 . The fertil- To determine the inheritance pattern of divergent coloration, we generated another F 1 intercross family (F H2 ) between a male S B morph and a female S V morph that were not used in previous crosses, in May 2014. After hatching, the progeny were reared until maturity in a 180-L tank. From May to June 2016, we created a backcross (BC) population from artificially fertilized eggs using sperm from the father individual and eggs from 18 mature hybrid females and an F 2 hybrid population by crossing one mature F 1 hybrid male and 18 F 1 hybrid females. The progeny were reared in 60-L tanks (about 50 individuals per tank) under the light, temperature, and feeding F I G U R E 1 Divergent color phenotypes with or without countershading in Sarcocheilichthys inhabiting Lake Biwa, Japan. Photographs of a countershaded morph taxonomically described as S. variegatus microoculus taken from the side (a) and below (b). Photographs of a melanistic morph taxonomically described as S. biwaensis taken from the side (c) and below (d). A variety of body color patterns, especially in terms of yellow coloration, were observed for each morph. Photographs in (a) and (c) are provided courtesy of the Lake Biwa Museum

| Chromatophore quantification and melaninrelated gene expression
We counted the individual chromatophores on the scales of each morph and compared their dorsoventral distribution throughout the body. Like the zebrafish (Danio rerio), which is a model cyprinid fish in a wide variety of biological fields (Hirata et al., 2005), Sarcocheilichthys chromatophores include melanophores, xanthophores, and iridophores ( Figure 2c). We easily detected a dense, uniform sheet of iridophores; however, we did not analyze iridophores in this study because they were difficult to count and less relevant to the objective of this study. To observe melanized melanophores and xanthophores on the scales, six individuals of each experimental family (F V1 , F V2 , F B1 , F B2, F BM , and F H1 ) were anesthetized in a 0.05% 2-phenoxyethanol solution, and scales were collected from To confirm divergence in the expression patterns of melaninrelated genes between the two morphs, we selected two important genes regulating melanin production, tyrosinase (TYR) gene and dopachrome tautomerase (DCT). DCT is conserved as a singlecopy gene across teleost genomes (Braasch et al., 2007). Although TYR gene is also present as a single copy in the zebrafish genome, many teleosts have two copies of this gene (Braasch et al., 2007). Therefore, at this time, the possibility that Sarcocheilichthys fishes have two TYR copies cannot be ruled out completely. The elongation factor 1 alpha (EF1α) gene was used as an internal control. EF1α gene has stable expression level during development and across tissue types, and is recommended as housekeeping gene in zebrafish (McCurley & Callard, 2008).
We used six individuals from each experimental family (excluding F BM ) for the reverse transcription-quantitative PCR (RT-qPCR) analysis. Total RNA was extracted from 10 scales on the dorsal (A) and ventral (E) areas of each individual using the RNeasy Micro Kit (Qiagen). Genomic DNA was digested using the RNase-Free DNase Set (Qiagen). cDNA was synthesized from extracted RNA (100 ng) using the High-Capacity cDNA Reverse Transcription Kit (Life Technologies). To design primers for RT-qPCR, the partial F I G U R E 2 Color phenotypes of a pure cross of the S. variegatus microoculus (S V ) morph, a pure cross of the S. biwaensis (S B ) morph, an F 1 hybrid cross, and an F 2 hybrid population between the two morphs under common garden rearing (a). The frequencies of individuals expressing either S V or S B color phenotypes in a backcrossed (BC) population and an F 2 intercross population, respectively (b). Chromatophores (Me: melanophore; Xa: xanthophore; Ir: iridophore) on the scales of Sarcocheilichthys fish (c)

Numbers of cells / mm 2
All statistical analyses were performed using R version 4.0.2 (R Core Team, 2020). Significance was determined at a level of p < 0.05.

| Population genomic and genome-wide association study (GWAS) analyses
To assess the degree of genome-wide divergence between the two morphs, we performed population genomic analysis of coexisting population using double-digest RAD-seq (ddRAD-seq; Peterson et al., 2012). Sixteen individuals of each morph were collected around off Onoe in Lake Biwa (see above). This site has been known as the co-occurring area of both morphs (Hosoya, 1982;Komiya et al., 2011). Most of these individuals were used in previous studies (Komiya et al., 2011(Komiya et al., , 2014. Genomic DNA was extracted from the pectoral or caudal fins using DNeasy Blood & Tissue Kit (Qiagen). A ddRAD library was prepared according to Peterson's protocol with slight modifications (Sakaguchi et al., 2015). In this study, genomic DNA was digested using MseI and BglII (New England Biolabs) and the library size was selected into 350-to 550-bp fragments by running on a 2% agarose gel. The samples were sequenced on the Illumina HiSeq X Ten with 150 bp paired-end reads. Finally, 41,804 SNPs were obtained and used in subsequent population genetic analyses.
Clustering analysis was performed by ADMIXTURE v.1.3.0 (Alexander et al., 2009). ADMIXTURE was run with K assumption ranging from 1 to 4. The results were summarized and graphically displayed using CLUMPAK (Kopelman et al., 2015). Cross-validation (CV) scores were used to determine the best K value. Next, the principal component analysis (PCA) was performed using the R package adegenet v.2.1.3 (Jombart, 2008;Jombart & Ahmed, 2011). We also calculated F ST value between the two morphs by GenoDive v.3.04 (Meirmans, 2020). Significance tests were conducted by 10,000 rounds of permutations.
In addition, a GWAS was conducted to define the genetic determinant of the two color phenotypes using the ddRAD SNP dataset.
We used the method genome-wide efficient mixed-model associa-

| Genetic crosses
All individuals within S V families (n = 26 for F V1 ; n = 30 for F V2 ) exhibited countershading, whereas those of S B families (n = 21 for F B1 ; n = 30 for F B2 ; n = 28 for F BM ) did not, as evaluated by appearance alone (Figure 2a). In each of these families, the number of melanophores tends to decrease from the dorsal to ventral regions (one-way repeated-measures ANOVA, F 4,20 = 5.1-219.7, p < 0.006; Figure 3b).
Although it did not differ significantly for A to C regions between the S V and S B morphs (nested ANOVA where family was nested within morphs, F 1,3 = 0.001-1.15, p > 0.36), the latter morph had a greatly higher number for D and E regions than the former morph, which nearly lack melanophores in ventral regions (F 1,3 = 20.8, p < 0.02 for D, F 1,3 = 213.2, p < 0.0001 for E). Although xanthophore numbers differed significantly among regions within each family (one-way repeated-measures ANOVA, F 4,20 = 3.9-28.8, p < 0.018), we observed neither a clear dorsoventral pattern nor a difference between morphs in the xanthophore numbers ( Figure 3b). In addition, there was no significant difference between the two morphs except for D region (nested ANOVA where family was nested within morphs, The offspring color ratios of hybrid crosses indicated that color divergence between these morphs was likely controlled by a single locus, that is, a two-allele Mendelian inheritance pattern. All progeny of the two F 1 hybrid families (n = 29 for F H1 ; n = 32 for F H2 ) exhibited countershading (see Figure 2a). Patterns in the overall distribution of melanophores and melanin-related gene expression were almost identical between the F H1 family and S V families (F V1 and F V2 ), suggesting that the allele for countershading is dominant (Figures 3, 4). In a backcrossed population, 102 and 107 individuals showed countershaded and melanistic color phenotypes, respectively, which was not significantly different from the expected 1:1 ratio (two-tailed exact binomial test, p > 0.78; Figure 2b). In an F 2 intercross population, 364 and 114 individuals showed countershaded and melanistic color phenotypes, respectively, which was not significantly different from the expected 3:1 ratio (two-tailed exact binomial test, p > 0.59; Figure 2b).

| Genetic divergence in coexistence populations
All of our population genetic analysis results indicated a lack of genome-wide divergence and population structure between the two morphs in a region where both morphs coexist. In the clustering analysis with ADMIXTURE using 41,804 SNP loci from ddRAD-seq, K = 1 was estimated to be best fit because its CV score was the lowest ( Figure S1). The pairwise estimate of genetic differentiation between the two morphs was very low (F ST = 0.001) and not significant (p = 0.079 after 10,000 permutations). The PCA results were consistent with these findings (Figure 5a). The levels of expected heterozygosity (H E ) and nucleotide diversity (π) were almost similar between the two morphs, respectively (H E and π: 0.286 and 0.296 for S V morph; 0.284 and 0.295 for S B morph).
The GWAS result revealed that each one SNP on two RAD loci 29,926 and 14,635 (p = 1.04 × 10 -6 and p = 0.039, respectively, following the Bonferroni correction) was significantly associated with color phenotype, showing an extremely highly significant result for the former locus (Figure 5b, see also Figures S2 and S3). We also detected statistically significant genotypic linkage disequilibrium between the two loci (p = 0.0009 by GENEPOP online version 4.7; Raymond & Rousset, 1995). Figure 3a) were selected for RT-qPCR.

| D ISCUSS I ON
One objective of this study was to determine whether the diverged color phenotypes of the Sarcocheilichthys system in Lake Biwa have a genetic basis. Under common garden rearing, color divergence was clearly shown to be under strict genetic control and not environmentally driven plasticity (e.g., Hochkirch et al., 2008;Roulin, 2014).
Although both morphs tended to decrease from the dorsal to ventral regions both at the levels of melanophore distribution and gene expression which regulates melanin production, this dorsoventral gradation of melanin pigmentation on the body differed significantly and showed discontinuity between the two morphs. The S V morph exhibited a near-total absence of melanized melanophores in the ventral region; conversely, the S B morph is characterized by melanophores occurrence and melanin production throughout the body.
Thus, our findings clarified that the S V morph exhibits sharp countershading, whereas the S B morph exhibits a strong (but not complete) tendency toward to lose its dorsoventral countershading. By contrast, xanthophore showed a complex distribution pattern throughout the body and their density differed widely among families of the same morph. For example, individuals of the family F V1 had many more xanthophores than the other families. Within this system, we observed considerable, continuous variation in the amount of yellow coloration in each morph (see Figure 1).
Both morphs examined in this study, which are treated as taxonomically separated species that are endemic to Lake Biwa, are presumed to have evolved from the ancestral riverine (sub)species Sarcocheilichthys variegatus variegatus, which is widely distributed in western Japan, except in Lake Biwa (Hosoya, 1982). Phylogeographic analyses using neutral DNA markers have suggested that the Lake Biwa population was created by multiple colonizations of this lake by riverine lineages, although a reconsideration of the taxonomic status of these Sarcocheilichthys fishes is necessary (see Komiya et al., 2014). Because the ancestral riverine lineages exhibit countershading, at least visually, the melanistic phenotype of the S B morph must have evolved from the countershading phenotype of the ancestral lineages.
The results of our F 1 hybrid, BC, and F 2 hybrid crosses indicated that the body color divergence was likely controlled by a single-locus, two-allele Mendelian inheritance pattern, and the melanistic allele of the S B morph is recessive, which causes a melanistic phenotype when homozygous. In several mammals and birds, it is known that melanism results from allelic variation in one of two genes, melanocortin receptor 1 (MC1R) and an antagonist for MC1R, agouti-signaling protein (ASIP), which are involved in the melanocortin system, a complex neuroendocrine signaling mechanism involved in numerous physiological processes including vertebrate pigmentation (e.g., Janssen & Mundy, 2017;Kingsley et al., 2009;McRobie et al., 2019;Nachman et al., 2003;Vidal et al., 2010). In mammals, ASIP binds to MC1R, which is the receptor for α-melanocyte-stimulating hormone (α-MSH), and binding of ASIP to MC1R lowers the ratio of eumelanin (dark/brown pigment) to pheomelanin (red/yellow pigment) produced in the melanocytes (mammalian melanophores) (reviewed by d 'Ischia et al., 2015). Thus, spatially regulated ASIP expression generates the countershaded pattern by driving a switch between the production of chemically distinct melanins in melanocytes in dorsal and ventral regions. Although fish melanophores synthesize only eumelanin, and their cellular and biochemical mechanisms are different from those of mammalian species, dorsoventral gradients of melanin pigmentation in many fishes also appear to depend on a dorsoventral expression gradient of ASIP1, the fish ortholog of mammalian ASIP (reviewed by Cal et al., 2017;Irion & Nüsslein-Volhard, 2019). A recent study reported that ASIP1 knockout zebrafish created using CRISPR-Cas9 technology exhibited a noncountershaded phenotype similar to melanism (Cal et al., 2019). By contrast, the gene expression of two melanocortin receptors (MC1R and MC5R), but not ASIP1, is higher in darker, more pigmented areas of the Japanese flounder (Matsuda et al., 2018). Each of these genes may be a putative genetic factor involved in color divergence within the Sarcocheilichthys system.
Population genetic analyses using 41,804 RAD-SNP markers indicated no genome-wide genetic differentiation between the two Sarcocheilichthys color morphs. This result strongly supported the findings of a previous study that was based on a restricted dataset (i.e., partial nucleotide sequences of mitochondrial cytochrome b gene and 14 microsatellite loci; Komiya et al., 2011). The S B morph has been taxonomically described as an independent species (S. biwaensis), primarily based on its body (and fin) color, body shape, and head morphology (Hosoya, 1982). Indeed, the S B morph has a longer head and deeper body than its congener. However, the ecomorphological traits of the S V morph (S. variegatus microoculus), which utilizes several types of bottom environments including sandy, pebbly, and rocky areas, show a remarkable continuous variation in relation to its habitat types, and its head and body shapes in rocky areas largely overlap with those of the S B morph (Komiya et al., 2011). These traits of both morphs inhabiting rocky areas are considered to be advantageous for capturing cryptic and/or attached prey in rocky areas with complex structure (Komiya et al., 2011). Thus, sympatric S B and S V morphs exhibit similar morphologies, except for body coloration, and cannot be discerned by genome-wide DNA markers, causing the difficult issue of their status.
Generally, the speciation process is continuous, with sympatric divergence occurring along a speciation continuum that ranges from a relatively homogeneous population containing multiple phenotypes (intraspecific polymorphism) to reproductively isolated sister species with clear genetic differentiation (Hendry et al., 2009;Nosil, 2012). Therefore, the panmictic status of the Sarcocheilichthys system suggested by genome-wide analysis lies somewhere along this speciation continuum, between genetic color polymorphism within a single population and development of some reproductive isolation mechanism between color morphs indicating that the isolation is either relatively recent and has not yet accumulated genetic differences (i.e., early stage of speciation or incipient speciation).
Although degree of reproductive isolation is typically used as a measure of how far the process of speciation has proceeded, premating isolation between these color morphs, that is, whether or not assortative mating exists between divergent color phenotypes under natural conditions, has not been well examined in the wild. However, we have a limited evidence of weak or lacking premating isolation between the color morphs. In a preliminary artificial crossing experiment using mixed sperm from multiple males and mixed eggs from multiple females of the S V morph collected in a rocky area, some offspring presented the S B color phenotype (Kokita et al., unpublished data). This phenomenon strongly suggests that some S V individuals have heterozygous genotypes (a dominant S V allele and a recessive S B allele) for the causal gene underlying color divergence, considering the genetic inheritance pattern of body color as mentioned above. Although the scope of this study did not include discerning intraspecific genetic polymorphism or incipient speciation, we assert that the maintenance of genetically determined color morphs within a homogenous population is likely for this system.
In any case, it is possible that the only gene specifically associated with body color differs between the two color morphs. Our GWAS actually detected that a single SNP on RAD locus 29,926 was strongly associated with the phenotype of each morph. Because RAD locus 29,926 did not show complete genotype-phenotype association based on a single-locus, two-allele Mendelian inheritance pattern, with the S B morph having the recessive allele (see Figure 5b), this result does not indicate that this SNP lies within a region in complete linkage disequilibrium of the causal locus underlying phenotypic variation between these morphs. However, there is a strong possibility that this SNP will exist near a locus potentially associated with the color divergence. Because of the genome-wide synteny between cyprinid fish species belonging to the subfamily Gobioninae and zebrafish (Kakioka et al., 2013), we performed a BLAST search of the two RAD consensus sequences (29,926 and 14,635; see Figure S3) containing significant SNPs against the zebrafish reference genome (GRCz11). However, the highest E-value of the BLAST hit sequences was 0.13; therefore, we were unable to find any homologous regions or syntenic regions. We are currently conducting de novo whole-genome sequencing and assembly, quantitative trait locus (QTL) mapping using the BC cross created in this study, comparative transcriptome analysis using RNA-Seq from scale skin tissues, and GWAS analysis using whole-genome resequencing data of wild-caught individuals. These genetic and genomic studies are likely to identify the actual causal gene and mutation for the loss of dorsoventral countershading in the S B morph in the near future.
Further population genetic analyses using the causal gene or its mutation will provide insight into the state along the speciation continuum of this Sarcocheilichthys system. Whether the Sarcocheilichthys system examined in this study represents intraspecific polymorphism or two incipient species associated with some degree with premating isolation, this system can be a potentially good model for investigating the functional significance of countershading and the eco-evolutionary drivers that cause its loss in the aquatic world, where these such studies are rare. As mentioned above, the melanistic morph is a rock-dwelling specialist that inhabits darker substrates; conversely, the countershaded morph utilizes a wide variety of bottom environments including sandy and pebbly zones with lighter substrates (Komiya et al., 2011). Therefore, the previous study speculated that the uniform dark coloration of the melanistic morph acts as camouflage, albeit briefly, in poorly lit rocky environments (Hosoya, 1982;Komiya et al., 2011). Certainly, since variation in countershading is dependent on the light environment, a dark average tone has often been associated with dark, closed environments (e.g., Allen et al., 2012;Caro et al., 2011;Cuthill et al., 2017). Most such studies have focused on camouflage as an antipredator adaptation (Cuthill, 2019;Cuthill et al., 2017), and the melanistic morph of Sarcocheilichthys may also have evolved to prevent its detection and/or recognition by predators. The littoral areas of Lake Biwa contain native piscivorous fishes such as the rockdwelling catfish Silurus lithophilus and piscivorous chub Opsariichthys uncirostris (Okuda et al., 2014), as well as the carnivorous soft shell turtle Pelodiscus sinensis (Matsui, 2012;Takigawa et al., 2020); these aquatic animals are potential predators of Sarcocheilichthys. Because Sarcocheilichthys fish usually swim biased to the subbenthic column (Nakamura, 1969), they face predation threats from the side and above. The complicated bottom environments of the littoral regions of Lake Biwa, including sandy, pebbly, and rocky area, were established about 0.4 million years ago (Meyers et al., 1993). Thus, uniformly dark coloration may offer an evolutionary advantage through matching the visual background in poorly lit, closed rocky environments. Conversely, according to the primary explanation for the evolution of countershading, the countershaded morph may attain concealment by reducing shadows (or background matching) in welllit environments like sandy and pebbly areas with lighter substrate. Therefore, the melanistic morph appears to have a cryptic disadvantage in these areas.
However, this scenario cannot explain the distributional patterns of the countershading morph that utilizes a wide variety of bottom environments, especially the rocky environments. The two morphs co-occur in rocky areas, showing spatial variation in their frequencies, with the melanistic morph representing 10-30% of the total population in these areas. (Nakamura, 1969). Our qualitative observations using SCUBA equipment in 2015 and 2019 confirmed that these color morphs certainly occur together in multiple rocky areas near two localities (codes 7 and 13; Komiya et al., 2014) (Kokita et al., unpublished data). However, although these color morphs co-occur, they may utilize different microhabitats within rocky areas. Rocky areas in the littoral zone in Lake Biwa exhibit great structural and environmental complexity, and there are many vacant, sandy spaces among the rocks, creating patchy, well-lit microhabitats. These contrasting microhabitats within the rocky zone may influence the relative fitness of each morph. Gene flow can also play an important role in maintaining polymorphisms between populations, regardless of selective processes acting within population (e.g., Bittner & King, 2003;Saccheri et al., 2008). Komiya et al. (2014) suggested the extensive gene flow among local populations of the countershaded morph inhabiting different bottom environments in Lake Biwa.
If this phenomenon represents intraspecific genetic color polymorphism within a homogenous population, several different and complex mechanisms may contribute to maintenance of this stable polymorphism (see Gray & McKinnon, 2007). At this time, we forgo a detailed discussion because it is just speculation.
In conclusion, the results of this study confirmed the coexistence of genetically determined color morphs with and without countershading in the Sarcocheilichthys system of Lake Biwa. This color divergence is likely controlled by a single locus, with melanistic color as the recessive alleles. In general, hypotheses explaining the evolution of color polymorphism are often based on the assumption that color morphs are genetically determined (Gray & McKinnon, 2007;White & Kemp, 2016). Therefore, this novel system provides an ideal opportunity for testing hypotheses about the functional aspects of the countershading and its loss in aquatic environments, and for understanding the mechanisms underlying the evolutionary maintenance of color polymorphism in the spatially heterogeneous littoral environments of this lake. Further studies of this system, integrating ecological, genetic, and genomic analyses are needed.

ACK N OWLED G M ENTS
We thank Y. Hayasaki and N. Kinoshita for their help in conducting the rearing and genetic experiments. We are also grateful to T.
Komiya and fishermen of Lake Biwa for their assistance with the field collections. Computations were partially performed on the NIG supercomputer at ROIS National Institute of Genetics. The present study was partly supported by JSPS KAKENHI Grant numbers 26291079, 17H03720, and 20H03009.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflicts of interest.