Oncogenic allelic interaction in Xiphophorus highlights hybrid incompatibility.

Mixing genomes of different species by hybridization can disrupt species-specific genetic interactions that were adapted and fixed within each species population. Such disruption can predispose the hybrids to abnormalities and disease that decrease the overall fitness of the hybrids and is therefore named as hybrid incompatibility. Interspecies hybridization between southern platyfish and green swordtails leads to lethal melanocyte tumorigenesis. This occurs in hybrids with tumor incidence following progeny ratio that is consistent with two-locus interaction, suggesting melanoma development is a result of negative epistasis. Such observations make Xiphophorus one of the only two vertebrate hybrid incompatibility examples in which interacting genes have been identified. One of the two interacting loci has been characterized as a mutant epidermal growth factor receptor. However, the other locus has not been identified despite over five decades of active research. Here we report the localization of the melanoma regulatory locus to a single gene, rab3d, which shows all expected features of the long-sought oncogene interacting locus. Our findings provide insights into the role of egfr regulation in regard to cancer etiology. Finally, they provide a molecular explainable example of hybrid incompatibility.

In the late 1920s, three investigators, Myron Gordon, Georg Haeussler, and Kurt Kosswig, independently found that hybrids between two distant Xiphophorus fish species, Xiphophorus maculatus (southern platyfish), and Xiphophorus hellerii (green swordtail), develop spontaneous and lethal pigment cell tumors that were later determined to be melanoma (1–3). Since its establishment, this model system has been intensively studied to assess the underlying genetic contributions to tumor etiology. The development of hybridization-induced tumor has been viewed as a representation of the genome incompatibility hypothesis known as the Bateson–Dobzhansky–Muller (BDM) model (4–8). The BDM model states that negative epistatic interactions in hybrids serve as the molecular genetic mechanisms underlying genome incompatibility and is associated with problems in hybrid fitness. Although BDM incompatibility was identified in a few model organisms (9), Xiphophorus and mice represent the only vertebrate systems in which the incompatible loci have been identified (10, 11).

This Xiphophorus interspecies hybrid system, now termed the “Gordon–Kosswig–Anders (GKA) model,” was first described by Gordon and Kosswig in 1920s. This model employs crossing of X. maculatus and X. hellerii to produce F1 interspecies hybrids. X. maculatus exhibits a nevus-like pigmentation pattern in its dorsal fin (spotted dorsal, Sd), while X. hellerii does not exhibit this trait. In the F1 hybrid, the Sd pigmentation pattern becomes expanded, with melanin pigmentation covering the entire dorsal fin due to melanocyte hyperplasia (12, 13). Backcrossing the F1 hybrid to the X. hellerii parent leads to three distinct phenotypes among the backcross (BC) progeny that follow Mendelian distributions: 25% of hybrids exhibit hyperplasia of pigmentation pattern as observed in the F1 hybrid, 25% exhibit lethal and invasive nodular exophytic melanoma, and the remaining 50% of progeny do not display a black pigmentation pattern (14) (Fig. 1). The hybridization-induced disease observed in Xiphophorus interspecies hybrid represents a type of genetic incompatibility. In 1950s, Anders argued this spontaneous tumorigenesis is due to segregation of two loci from X. maculatus; one was named Tu for “tumor” and another locus was named R for “repressor” or Diff for “differentiation” [hereafter referred to as R(Diff)]. These concepts led to what we now know as oncogenes and tumor suppressors (15, 16). In the late 1980s, it was shown the Tu gene encodes a mutant duplicate copy of egfr, and this gene was named as Xiphophorus melanoma regulatory kinase (xmrk) (17, 18). The xmrk oncogene is tightly linked to or part of the Sd locus and controls melanocyte proliferation.

Fig. 1.

GKA model crossing scheme. The crossing scheme shows the Xiphophorus species used to produce F1 and BC1 interspecies hybrids. X. maculatus Jp163 A and X. hellerii are used to produce F1 hybrids artificially. The F1 hybrids are subsequently backcrossed to X. hellerii to produce BC hybrid progeny. BC hybrids exhibiting melanocyte hyperplasia and heterozygous R(Diff) are used as parents for next-generation BC.

In addition, the GKA model offers us a natural two-hit melanoma model wherein the oncogenic effect of xmrk can be fully eliminated by a regulatory locus that must have coevolved with xmrk (19). The EGFR gene is one leading oncogene of many human cancers (20). It is preproliferative and is an upstream activator of BRAF and NRAS signaling, which are the driver oncogenes in over 50% of all human melanomas (21). This evidence promotes the characterization of R(Diff) as having significant implications in cancer etiology. Therefore, identifying the R(Diff) gene will highlight the genetic interactions underlying the tumor induction due to hybrid incompatibility and in addition may forward novel molecular target(s) in regulation of EGFR function for human disease control.

However, the R(Diff) gene has not been identified. Previous effort to define the R(Diff) locus forwarded a chromosomal region that includes cdkn2ab, an ortholog of the human tumor suppressor genes CDKN2A and B, which is mutant in 10% of high-melanoma-risk families (22–27). In this study, we found cdkn2ab is not the R(Diff) gene but is tightly linked to it. More importantly, we have identified the long-hypothesized R(Diff) locus to a strong candidate gene, rab3d.

Results

A Region on Chromosome 5 Determines Malignancy of Melanocytic Lesions.

Hybrids between X. maculatus and X. hellerii exhibit enhanced pigmentation in their dorsal fin (i.e., melanocyte hyperplasia) due to hemizygosity for the xmrk oncogene and heterozygosity of the R(Diff) locus [i.e., xmrk ; R(Diff) ]. This phenotype and genotype are also present in ∼25% of BC hybrid progeny (Fig. 1). Due to meiotic recombination, successive backcrossing of such animals to X. hellerii should stepwise reduce heterozygosity in the advanced BC progeny and finally result in a xmrk ; R(Diff) isogenic line (i.e., introgression). In this manner, advanced BC (BCn) fish were produced in order to validate the candidate gene cdkn2ab (chromosome 5, 15.8 Mbp) as a locus carrying R(Diff) function. The maternal parent (i.e., interspecies hybrid) for each successive BC generation was genotyped for inheritance of cdkn2ab and only heterozygous hybrids (i.e., xmrk , cdkn2ab ) exhibiting the enhanced dorsal fin pigmentation were selected for a next round of backcrossing. χ2 tests were performed on genotyping data of each variant site collected from 90 BCn hybrids that exhibited melanocyte hyperplasia to locate genes that displayed ancestral allele linkage disequilibrium and predominately showed a heterozygous inheritance pattern (). Two genomic regions were found to correlate to the hyperplasia phenotype (Fig. 2): one on chromosome 21 which encompasses xmrk, the melanoma driver oncogene that induces melanocyte proliferation, and a second region on chromosome 5 that corresponds to a previously mapped R(Diff) region of 5.8-Mbp region harboring the candidate cdkn2ab gene (Fig. 2). An average loss of 50% heterozygous loci per BC generation is expected. Therefore, one expects to see a heterozygous region only accounting for an average of 5.5 and 2.7 Mbp of the 700-Mbp genome for BC7 and BC8 individuals, respectively. However, haploid maps produced from these advanced BC animals show the heterozygous content to be much higher than this expectation. In addition, heterozygous loci are predominantly surrounding the cdk2ab region (), suggesting selection of individuals that exhibited melanocyte hyperplasia and a genotype of cdkn2ab for further backcrossing coselected an adjacent locus on chromosome 5, rendering higher-than-expected heterozygosity in BCn hybrids. These observations indicate that cdkn2ab itself is less likely to be R(Diff), while the coselected locus with cdkn2ab is the true R(Diff).

Fig. 2.

Genetic mapping of heterozygous loci in advanced BC hybrids. A total 90 BC hybrids (BC2 to BC8) that exhibited dorsal fin melanocyte benign hyperplasia and were produced by crossing cdkn2ab genotyped hybrid (cdkn2ab ) with X. hellerii (cdkn2ab ). (A) Manhattan plot showing –log10P value (χ2 test) across the genome. The y axis represents –log10P value and the x axis represents amplicon chromosomal coordinates, which are labeled as red or blue. Only the −log10 P value of loci that exhibited higher X. maculatus allele frequency is plotted due to introgression. The light blue dashed line represents a P value of 0.01 that is suggestive of statistical significance. χ2 test P values were corrected using Bonferroni method across the genome-wide data. The red dashed line represents adjusted P value of 0.05 corresponding to 3.6 × 10−5. (B) A zoom-in view of the chromosome 5 10,029,746- to 15,895,617-bp region. This region is highly correlated to the pigmentation phenotype observed in BC hybrids.

Genetic Mapping of a Mutant EGFR Regulator Locus.

To locate the R(Diff) candidate gene(s) independent of artificial selection of genetic marker, we produced BC interspecies hybrid progeny of the GKA model, that is, X. hellerii (Rio Sarabia) × [X. maculatus Jp163 A × X. hellerii (Rio Sarabia)] and performed targeted genomic sequencing on BC progeny that developed two distinct melanocyte phenotypes: benign hyperplasia and melanoma tumor. Association analyses between parental allele inheritance and pigment cell phenotypes were performed on a total of 66 BC1 progeny (19 hyperplasia and 47 melanoma). Both melanocyte phenotypes are known to be xmrk-dependent; therefore, no association between xmrk and melanocyte phenotypes was observed, as expected (Fig. 3). Genotypes of four linked polymorphic sites (10,582,852 10,582,855 10,582,868, and 10,582,870) on chromosome 5 are the most significantly correlated to melanocyte phenotypes, where individuals exhibiting benign pigment cell lesions inherited both parental alleles, while melanoma-bearing individuals only inherited the X. hellerii alleles (Fig. 3).

Fig. 3.

Genetic mapping of the R(Diff) locus. BC1 hybrids were produced by crossing F1 hybrid to X. hellerii. Pigmented hybrids were classified into two categories independent of any molecular marker. Sixty-six hybrids that include 19 exhibiting dorsal fin melanocyte hyperplasia and 47 displaying melanoma were genotyped. The number of animals of each genotype does not reflect the BC hybrid phenotypical distribution because the Xiphophorus Genetic Stock Center preferentially collects tumor-bearing fish for research purposes. (A) Manhattan plot showing –log10 P value of χ2 test across the genome. The y axis represents –log10 P value and the x axis represents amplicon order on each chromosome, which is labeled as red or blue. the light blue dashed line represents P value of 0.01 that is suggestive of statistical significance. χ2 test P values were corrected using Bonferroni method across the genome-wide data. the red dashed line represents adjusted P value of 0.05 corresponding to 3.6 × 10−6. Both melanocyte hyperplasia and melanoma are driven by xmrk, and therefore chromosome 21 is not related to separation of the two melanocyte phenotypes. A strong peak corresponding to linkage disequilibrium is found on chromosome 5, with polymorphisms located at 10,582,852, 10,582,855, 10,582,868, and 10,582,870 bp exhibited the top –log10 P value. These polymorphisms are adjacent to three gene models: dand5, tspan, and rab3d. (B) Chromosome 5 haploid maps of all BC hybrids. Blue and red bars represent genotypes in term of parental allele inheritance (red: heterozygous for both parental alleles; blue: homozygous for X. hellerii allele). Their locations on the bar graph correspond to amplicon physical location. Colored lines underneath the haploid maps represent phenotypes, with the red line corresponding to pigment cell lesion and the blue line corresponding to melanoma. A locus located between 10,515,844 and 10,617,563 bp is free of recombination in both groups of hybrids. Hybrids exhibiting melanocyte hyperplasia inherited both ancestral alleles in this locus, and tumor-bearing hybrids inherited only a recurrent parental (i.e., X. hellerii) allele.

Haploid maps produced from all BC individuals supported the result of association analyses and forwarded a region (i.e., 10,515,844 to 10,617,563) that is free of chromosomal cross-overs in both cohorts (Fig. 3). The R(Diff) locus is assigned to this 101.7-kbp region (Fig. 3). As expected, this locus maps to the vicinity (5.2 Mbp) upstream of cdkn2ab.

The rab3d Gene Is the Functional Carrier of the R(Diff) Function.

The 101.7-kbp R(Diff) locus encodes three gene models on the reverse strand: differential screening-selected gene aberrative in neuroblastoma (DAN) domain family member 5 (dnad5), tetraspanin-1 (tspan1), and ras-related protein Rab-3D (rab3d). Expression of a gene is the prerequisite for display of genetic function. To determine the gene expression pattern of these three candidate genes and to assess gene expression changes in parental nevus-like dorsal fin melanocyte spots, melanocyte hyperplasia, and melanoma tumor, transcription profiling was performed on these tissues. The rab3d gene is the only gene expressed among the three candidates in melanocyte spots of parental X. maculatus (three pools of dorsal fin spots), dorsal fin melanocyte hyperplasia (n = 22), and melanoma tumor (n = 22) of BC interspecies hybrids (Fig. 4). Therefore, rab3d serves as the only gene that can carry the R(Diff) function.

Fig. 4.

Gene expression profiling of the R(Diff) locus. Bar graph showing dnad5, tspan1, and rab3d expression in X. maculatus dorsal fin spots (Sd), benign hyperplasia (Bn), and malignant tumors (Mm) of BC hybrids. The y axis represents library size-normalized read counts. Error bars represent SD. Differential expression analyses between different genotypes were performed using DESeq2, and P values are corrected using FDR.

In addition, differential expression analysis of rab3d among the parental dorsal fin spots, interspecies hybrid melanocyte expansion, and in melanoma showed rab3d expressed at the same level between the melanocyte spots and hyperplastic pigmented cell lesions (adjusted P value = 0.81; Fig. 4), and expressed at a higher level in melanocyte hyperplasia (adjusted P value = 5.4 × 10−5) than in melanoma of BC fish.

The two parental rab3d alleles differ by an Asn residue in X. maculatus and Lys in X. hellerii at the C-terminus tail downstream of the P-loop domain that harbors the guanosine triphosphatase (GTPase) activity (i.e., Asn/Lys-204; Fig. 5). Although this Asn/Lys site locates in the hypervariable C terminus of rab3d, comparative genomics showed that different from X. maculatus the Lys is conserved in human and 60% of all fish species analyzed (i.e., the next dominating allele is Arg) that include another xmrk-null Xiphophorus species, Xiphophorus couchianus, and other Poeciliidae fish, suggesting the X. maculatus allele is a genetic outlier (Fig. 5 and Dataset S1).

Fig. 5.

Sequence alignment of rab3d genes. Protein sequence comparison between X. hellerii rab3d and X. maculatus, X. couchianus, and Poeciliid fish Poecilia mexicana. Different amino acids between X. hellerii and X. maculatus and between X. hellerii and P. mexicana are labeled in red and black, respectively. Functional domains are labeled with black lines underneath protein sequences, with thinner lines linked functional domain in tertiary structure. The only amino acid change between X. maculatus to other Xiphophorus is Lys-204 > Asn-204 mutation. The Lys is conserved, and only X. maculatus encodes an Asn.

Discussion

The finding that rab3d is the only expressed gene located in the tumorigenesis-determining locus on chromosome 5 in Xiphophorus genome, forwarded by both genetic mapping and transcriptomics, concludes a search for the hypothetical locus R(Diff) that has been ongoing for over five decades (22, 25, 28–31). As a viviparous fish where embryonic development occurs within the female, the technical hurdles involved in transgenesis have not yet been overcome for Xiphophorus. Therefore, genetic manipulation of rab3d cannot be performed in the Xiphophorus system until our current development of Xiphophorus transgenesis is proven successful and efficient. Although a readily available system for such a test is currently not available, the negative epistasis between xmrk and rab3d hallmarks cancer, in addition to hybrid lethality (9) and sterility (32), an innovative mechanism for decreasing hybrid fitness, and reinforcing speciation.

Epitasis underlying human disease can be elucidated by investigating mechanisms that “evolutionary mutant models” developed to cope with similar mutations as in human disease and/or produce adaptive phenotypes that are similar to human disease (19, 33). X. maculatus is one such species where oncogenicity of a mutant EGFR (i.e., xmrk) is compromised by a regulatory allele R(Diff). EGFR is one of the most prevalent oncogenes exhibiting mutation and/or dysregulation in many varied human cancers (34–53). Despite over 40 y of effort in attempting to inhibit EGFR by blocking its kinase activity, and development of four generations of small molecules, or monoclonal antibodies, success in disease control is very limited to three cancer subtypes, that is, nonsmall cell lung cancers with kinase-activating mutations in EGFR (54–56), ∼10% of metastatic colorectal cancers (57, 58), and a subcategory of advanced head and neck cancers (59, 60). Current methodologies to inhibit EGFR attempt to directly block the adenosine 5′-triphosphate (ATP) binding pocket or ligand binding pocket or by targeting acquired mutations that lead to acquired resistance have all turned out to be inefficient in providing promising therapeutic benefit. Therefore, it is of utmost importance to revolutionarily reattack the question of EGFR-associated cancer etiology and identify reliable next-generation treatment strategies that enable disease control with higher response rates and lower resistance. The xmrk, originated from a gene duplication event (61), encodes a mutant EGFR that autodimerizes and activates downstream proproliferative pathways in a ligand-independent manner (62). The xmrk gene is a bona-fide oncogene because its ectopic expression leads to transformation and tumorigenesis of melanocytes in murine cells and medaka fish (63, 64). However, X. maculatus does not exhibit tumorigenesis, suggesting X. maculatus harbors a mechanism [i.e., R(Diff)] in suppressing the driver oncogene. Therefore, the identification of rab3d as gene exhibiting R(Diff) regulatory function delineate mechanisms of how X. maculatus counteracts the deleterious EGFR mutant in its genome. The molecular mechanism underlying rab3d suppression of xmrk can lead to innovative strategies in developing next-generation EGFR inhibitors.

The rab3d genes showed both codon mutation between the two parental alleles and transcriptional differences between normal dorsal fin pigment cells, pigment cell hyperplasia, and melanoma tumor. These results suggest both structural and expression divergences between the two parental alleles of rab3d are essential to elucidate molecular interactions we can learn in order to advance our knowledge in control of EGFR and associated disease. There are two questions to answer regarding rab3d regulation of xmrk. First, what is the molecular nature of this xmrk–rab3d interaction? RAB3D is a ras-related small G protein with GTPase activity controlling exocytosis (65–67) and has been shown to regulate secretion of a broad range of molecules in different cell types (68–71). For example, RAB3D-dependent secretion of matrix metalloproteinase in macrophages is a prerequisite for macrophage recruitment to tumor cells (68), while secretion of the same molecule by tumor cells is a signal for tissue invasion and metastasis (72). We hypothesize that the X. maculatus allele of rab3d regulates cancer cell invasion by either mediating immune cell recruitment to tumor microenvironment or by hampering tumor cell secretion of molecules that facilitate metastasis. However, RAB3D function in cancer has not been clearly characterized despite efforts made investigating its function in cancer cell proliferation and metastasis in vitro (73–83). Currently there is no understanding of the cell population within the Xiphophorus tumor microenvironment (e.g., tumor cells, endothelial cells, cancer-associated fibroblasts, and immune cells) where expression of different parental rab3d alleles may change the fundamental outcome of xmrk expression (84, 85). RAB3D is also involved in cell membrane–associated protein dynamics. It has been shown that human RAB3D interacts with GOLM1 and selectively assist cytoplasmic EGFR recycling to cell membrane (83). Therefore, it is our second hypothesis that rab3d regulates xmrk function by directly controlling xmrk protein turnover and cellular localization where it displays its full activity. Second, what functionality change is associated with the amino acid difference between the two parental alleles? The change of Asn204 in X. maculatus to Lys204 in X. hellerii is located in the C terminus of RAB3D, where proper modifications (i.e., methylation and geranylgeranylation) are required for subcellular localization (86, 87). The Lys (X. hellerii) and Asn (X. maculatus) exhibit both physical differences (i.e., Lys is charged and Asn is uncharged) and varied posttranslational modifications. This amino acid change may alter protein hydrophobicity, affect RAB3D subcellular localization, and hinder efficient transportation of secretion granules and eventually affect on-site and dynamics of the above mechanisms. In summary, the conclusion of epistasis underlying the interspecific hybridization-induced tumorigenesis provides insights into a strategy in counteracting detrimental conditions.

Overall, vertebrate organisms that support the molecular mechanism proposed by the BDM model is only limited to mice (32) and now Xiphophorus fishes. The discovery of xmrk–rab3d genetic interaction underlying spontaneous tumorigenesis in interspecies hybrids poses an example showing that hybrid-induced disease can act as a mechanism that reduces hybrid fitness. Characterizing mechanism of RAB3D functional regulation of EGFR can lead to development of innovative EGFR regulation strategy.

Materials and Methods

Animal Model.

X. maculatus Jp163 A, X. hellerii (Rio Sarabia), and first-generation BC (BC1) animals used in this study were supplied by the Xiphophorus Genetic Stock Center (https://www.xiphophorus.txstate.edu/). X. maculatus Jp163A strain female fish were artificially inseminated with sperm from male X. hellerii (Rio Sarabia strain) to produce F1 interspecies hybrids. F1 hybrid males were then backcrossed to X. hellerii females to generate the BC1 animals. At dissection, all fish were anesthetized in an ice bath and upon loss of gill movement were killed by cranial resection. Organs were dissected into RNAlater (Ambion Inc.) and kept at −80 °C until use. All BC1 fish were kept and samples taken in accordance with protocols approved by Texas State University Institutional Animal Care and Use Committee (IACUC 2015107711).

The advanced BC interspecies hybrids (i.e., BC2 through BC8) were produced in an independent series of crosses in the Biocenter of the University of Wurzburg, Germany. F1 interspecies hybrids originated from the reciprocal cross: X. maculatus Jp163A males were mated to X. hellerii (Rio Lancetilla strain) females. The F1 hybrid females (cdkn2ab genotype) were then successively backcrossed to X. hellerii males to produce the advanced-generation of BC hybrids. For each generation of backcrossing, only the interspecies hybrid BC fish that exhibited benign pigment cell hyperplasia and had the cdkn2ab genotype were used to produce the next generation of BC progeny. Fin clips were collected from all advanced BC fish and stored in ethanol at 4 °C. All advanced BC fish were kept and samples taken in accordance with the applicable European Union and national German legislation governing animal experimentation. When needed, fish were killed by overanesthetization with MS222. These experiments were performed under authorization (568/300-1870/13) of the Veterinary Office of the District Government of Lower Franconia, Germany, in accordance with the German Animal Protection Law (TierSchG).

DNA and RNA Isolation.

Fin clip, or muscular tissue, was digested by Protease K at room temperature for 1 h. The lysate was then transferred to 2.0-mL collection tubes. DNA isolation was performed by a QIAcube HT (Qiagen) automated biosample isolation system, with reagent contained in the QIAamp 96 DNA QIAcube HT Kit. The isolation system is equipped with a robotic arm with eight pipettes. Each pipette is able to pick and eject pipette tips, self-clean, and transfer liquids between wells/columns, or between master reservoirs and wells/columns in standard 96-well plate formats. Each sample was independently maintained throughout the isolation process. Concentrations of DNA samples were measured using Qubit 2.0 fluorometer (Life Technologies) and adjusted for sequencing library preparation.

Dorsal fin spots, dorsal fin exhibiting benign hyperplasia, and melanoma tumors were excised from X. maculatus Jp163A, and BC interspecies hybrid. Tissue samples were homogenized in TRI-reagent (Sigma Inc.) followed by addition of 200 μL/mL chloroform, vigorously shaken, and subjected to centrifugation at 12,000 × g for 5 min at 4 °C. Total RNA was further purified using an RNeasy mini RNA isolation kit (Qiagen). Column DNase digestion at 25 °C for 15 min removed residual DNA. Total RNA concentration was determined using a Qubit 2.0 fluorometer (Life Technologies). RNA quality was verified on an Agilent 2100 Bioanalyzer (Agilent Technologies) to confirm that RNA integrity number scores were above 8.0 prior to subsequent gene expression profiling.

Genetic Variants Identification and Annotation.

To identify interspecies polymorphisms between the X. maculatus and X. hellerii, genomic DNAs of 4 X. maculatus and 4 X. hellerii were isolated. DNA samples were forwarded for genome shotgun sequencing library preparation using Illumina Nextera sequencing Library Prep Kit, followed by sequencing on HiSEq. 2000 (Illumina, Inc.) using 150-bp paired-end sequencing strategy. Raw sequencing reads were trimmed and filtered using a custom Perl script, and adapter sequences were removed from the sequencing reads. The reads were truncated based on similarity to library adaptor sequences using custom Perl scripts. Then, low-scoring sections of each read were removed, preserving the longest remaining sequencing read fragment (88). Filtered genome sequencing reads were mapped to the reference X. maculatus genome (GenBank assembly accession no. GCA_002775205.2) using Bowtie2 “head-to-head” mode to show mismatches between sequenced animals and reference genome (89). Alignment files were sorted using Samtools (90). Subsequently, pileup files were generated for each X. maculatus, and X. hellerii sample, and variant calling was processed by both BCFtools and VarScan for polymorphisms detection, with minimum variant locus coverage of 2 and a P value for variant detection of 0.05 for VarScan and variant genotyping call Phred score of 0 and alternative genotyping Phred score ≥ 20 for BCFtools (90–92). Only the variants that were identified by both pipelines were forwarded for further analyses.

To localize fixed variants between the X. maculatus and X. hellerii, homozygous loci of X. maculatus were compared to those of X. hellerii. Such loci were identified if all X. maculatus were homozygous for one allele and all X. hellerii homozygous for the alternative allele.

These fixed species-specific genetic variants were functionally annotated using snpEff (93). A genome database was created using the X. maculatus genome sequence and annotation files (GenBank assembly accession no. GCA_002775205.2). Each variant was queried to the genome database to determine if it was located in a genetic or intergenic region, and to determine what effect each variant may have on the peptide sequence structure.

Amplicon Sequencing, Data Filter, and Genotyping.

Variants between X. maculatus and X. hellerii were used as references to design specific capture probes for targeted genomic sequencing. Variants with very high sequencing depth were removed due to the possibility of locating them in repetitive sequences. Sequencing probes were designed to amplify regions surrounding genetic variants. To genetically map candidate R(Diff) loci in a region (chromosome 5: 10,000,000 to 16,000,000) identified in a previous study (22), 406 sets of probe were designed to reach a resolution of 14.8 kbp within the 6-Mbp locus; for sex-determining regions (chromosome 21: 23,750,000 to 26,250,000; and Dataset S2), 101 sets or capture probes were designed for 24.8-kbp definition genetic mapping; for the rest of the genome, 1,510 sets of probes were designed for genotyping and establish individual BC progeny haploid map at definition of 459 kbp ( and Dataset S2). Therefore, a total of 2,017 probe sets were produced for amplicon sequencing. Amplicons were custom-made using Illumina Genotype Ne library preparation kit, with i7 and i5 indices incorporated into adaptor sequences added to each end of PCR products amplified by capture probes. Pooled sequencing libraries were sequenced on Illumina MiSeq platform employing a 75-bp paired-end sequencing strategy (Illumina).

Sequencing adaptor contamination was first removed from raw sequencing reads using fastx_toolkit, followed by trimming of low-quality sections of each sequencing read. Low-quality sequencing reads were further removed from sequencing result (http://hannonlab.cshl.edu/fastx_toolkit/index.html). Processed sequencing reads were mapped to X. maculatus genome v5.0 (GenBank assembly accession no. GCA_002775205.2) using Bowtie2 (89). Mpileup files were made using samtools and genotyping was processed using both Bcftools and VarScan (90–92). Genetic variant call and genotype were required to be supported by both pipelines for further analyses (i.e., Bcftools: MAPQ ≥30, Phred score of genotype call = 0, with alternative genotype call Phred score ≥20; VarScan: MAPQ ≥ 30, P value <0.05, depth ≥20). Herein “genotype” refers to inheritance of ancestral alleles, with “heterozygous” meaning that a locus exhibited genetic material from both ancestors (i.e., X. maculatus and X. hellerii), and “homozygous” means that a locus exhibited genetic material from only the recurrent ancestor (i.e., X. hellerii). A haploid map was produced for each individual of BC progeny. To control the amplicon-sequencing-based genotyping result target specificity, only genotyping calls that locate less than 75 bp from the designed polymorphic sites and also supported by at least another variant genotyping call within a 75-bp range were kept. Qualified genotyping calls were subsequently ordered by chromosome name and their chromosomal location to produce the haploid map. For each BC hybrid individual, percentage of heterozygous loci is calculated as Heterozygous %= (number of heterozygous loci)/(total number of genotyped loci). Because BCn were selected for a marker on chromosome 5 (i.e., cdkn2ab) and chromosome 21 (i.e., Sd), percentage of heterozygous was calculated using all genotyped loci, or loci that are outside of chromosomes 5 and 21 ().

Calculation of Allele Frequency.

Genotyped variants loci of all BC hybrids were combined together first to yield a data table with the rows as chromosomal coordinates and columns as hybrid individuals. Heterozygous, designated as “0/1,” refers to inheritance of both X. maculatus and X. hellerii parental alleles, while homozygous, designated as “1/1,” refers to inheritance of only X. hellerii alleles. The X. maculatus allele frequency for each locus is calculated as

Linkage Analyses.

For advanced BC samples exhibiting benign melanocyte hyperplasia their genetic backgrounds are predominantly represented by the recurrent parental genome (X. hellerii). An average X. hellerii genome component per BC generation follows a rule determined by (1 − 0.5n+1), where n equals the BC generation. Therefore, every locus is expected to exhibit dominance of X. hellerii allele and therefore disequilibrium (). Since our previous studies had determined that the R(Diff) locus is heterozygous within the BC hybrids exhibiting melanocyte benign hyperplasia, we only plotted the −log10 P values of loci where X. maculatus allele frequency is higher and assigned −log10 P values of X. hellerii dominated loci arbitrarily to 0, in order to visualize dominantly heterozygous loci within the BC hybrids (). Because there is only one pigment cell phenotype of the BCn hybrids, numbers of heterozygous and homozygous individuals per variant site were used to form a one-dimensional contingency table and tested using a goodness-of-fit χ2 test. χ2 test P values were adjusted using Bonferroni correction.

For BC1 progeny, numbers of heterozygous and homozygous individuals were counted for each pigmentation phenotype (i.e., melanoma and pigment cell hyperplasia), and numbers of each genotype per phenotype group (i.e., tumor or pigment cell hyperplasia) were used to form a contingency table and subsequently tested using a χ2 contingency table test, with the null hypothesis that both genotypes’ distribution follows random assortment. χ2 test P values were adjusted using Bonferroni correction.

Gene Expression Profiling.

RNA sequencing was performed upon sequencing libraries construction using the Illumina TruSeq messenger RNA (mRNA) library preparation kit (Illumina, Inc.). RNA libraries were sequenced as 125-bp paired-end fragments using an Illumina Hi-SEq. 2000 system (Illumina, Inc.). Short sequencing reads were filtered using an in-house data processing pipeline (88). RNA-sequencing reads produced from three sets of pooled dorsal fin spots collected from X. maculatus, dorsal fin exhibiting melanocyte hyperplasia of 22 BC hybrid, and melanoma of 22 BC hybrid were produced. Sequencing reads were mapped to X. maculatus reference genome (GenBank assembly accession no. GCA_002775205.2) using STAR (94). Gene expression was subsequently profiled by counting number of sequencing reads that mapped to gene models annotated by Ensembl using RSEM (95). For data visualization, gene expression read counts were normalized to library size and were plotted as bar graph using R (v3.5.1). Differentially expressed genes were identified using R package DESeq2, with P value adjusted using the false discovery rate (FDR) method (96). FDR <0.05 was used to determine differential expression.

Annotation of X. Hellerii Allele of rab3d.

The genome sequences of X. maculatus chromosome 5 (GenBank assembly accession no. GCA_002775205.2) and X. hellerii chromosome 5 (GenBank assembly accession no. GCA_003331165.2) were aligned to each other using lastz (http://www.bx.psu.edu/miller_lab/dist/). The R(Diff) region alignment was extracted using genomic coordinates of R(Diff) for data visualization. To localize the X. hellerii rab3d gene the coordinates of rab3d exons and coding sequences for X. maculatus were transferred to the X. hellerii allele using the lastz alignment. Sequence comparisons between mRNA and genomic DNA, protein sequence comparisons between X. hellerii to other species were processed using BLAST (https://blast.ncbi.nlm.nih.gov/Blast.cgi).

Statistical Information.

Goodness-of-fit and contingency table χ2 tests were used to identify linkage disequilibrium for advanced BC (n = 90) and BC1 hybrids (n = 66), respectively. χ2 test P values were adjusted using Bonferroni correction.

Two-tailed t tests were used to test if sizes of BC7 (n = 23) and BC8 (n = 3) hybrid heterozygous loci are different from expected 5.5 and 2.7 Mbp, respectively.

Differential gene expression among parental dorsal fin melanocyte spots (n = 3 pooled samples), hybrid melanocyte hyperplasia samples (n = 22), and hybrid melanoma tumors (n = 22) were performed using the R package DESeq2 that implemented a modified Fisher’s exact test. Multiple test P values were corrected using the FDR method (96).