Animal use was approved by the Texas State University Institutional Animal Care and Use Review Board (IACUC protocol #2015107711). Fish utilized in this study were maintained in accordance with the applicable OLAW guidelines governing animal experimentation in the United States, and international legislation regulations governing animal experimentation.

Xiphophorus utilized in this study were bred and maintained in the Xiphophorus Genetic Stock Center¹. Xiphophorus maculatus (X. maculatus) JpA and JpB strains used for genome re-sequencing were at their 116th and 109th generation of inbreeding, respectively.

For DNA isolation, 4 male JpA and 4 male JpB fish were sacrificed by over-anesthetization with MS222 (0.06%). Fish tissues were digested with Proteinase K at 37°C for 1 hr. The lysate was then used for DNA isolation and purification using DNeasy Blood & Tissue Kit (Invitrogen).

For RNA isolation, fish were anesthetized in an ice bath and upon loss of gill movement were sacrificed by cranial resection. Skin and liver tissues were dissected directly into TRI reagent (Sigma Inc., St Louis, MO, United States) and flash frozen in an ethanol dry ice bath. For JpA, skin samples (n = 25), and liver samples (n = 6) were collected; For JpB, skin (n = 30) samples, and liver samples (n = 4) were collected for RNA isolation and sequencing library preparation. RNA isolation was performed following the Qiagen RNeasy RNA isolation protocol (Qiagen, Valencia, CA, United States). Tissue samples harvested from fish were first homogenized using a hand-held homogenizer in a 1.5 mL microcentrifuges tube while the sample remained frozen in TRI Reagent (Sigma Inc., St Louis, MO, United States). After homogenization, 300 μL of fresh 4°C TRI Reagent was added to the samples followed by room temperature incubation for 5 min. Chloroform extraction was performed by adding 120 μL chloroform and shaken for 15 s. Samples were centrifuged (16,100 rcf for 5 min at 4°C) for phase partitioning. The aqueous layer was transferred to a new 1.5 mL microcentrifuges tube and a second chloroform extraction performed (300 μL TRI Reagent, 60 μL chloroform). After extraction, nucleic acids in the aqueous phase were precipitated with 500 μL 70% EtOH in diethylpyrocarbonate (DEPC) treated water. The sample was then transferred to a Qiagen RNeasy mini spin column and on-column DNase treatment was performed for 15 min at 25°C. RNA samples were then washed and eluted in 100 μL RNase free water. RNA concentration was measured with a Qubit 2.0 fluorometer (Life Technologies, Grand Island, NY, United States). To further assess the RNA quality, an RNA integrity (RIN) score was determined using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, United States). All samples processed for RNA sequencing had a RIN score above 8.

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., San Diego, CA, United States) using 150 bp paired-end (PE) sequencing strategy with an average of 33 X genome coverage. Isolated RNA samples were forwarded for Illumina High-throughput Sequencing using the Illumina TruSeq mRNA Library Prep Kit on the HiSeq 2000 platform (Illumina, Inc., San Diego, CA, United States). Each RNA sample was used to construct a single library for sequencing (PE sequencing; See Supplementary Table 7 for details). For both genomic and transcriptomic sequencing, raw reads were trimmed and filtered using a custom Perl script and adapter sequences were removed from the sequencing reads (Garcia et al., 2012). The reads were truncated based on similarity to library adaptor sequences using custom Perl scripts (Garcia et al., 2012). Then, low-scoring sections of each read were removed, preserving the longest remaining sequencing read fragment. Sequencing statistics are included in Supplementary Table 7.

Filtered genome sequencing reads were mapped to the reference genome (X.mac V5.0; NCBI accession number: GCA_002775205.2) using Bowtie2 “head-to-head” mode (Langmead and Salzberg, 2012). Alignment files were further filtered by keeping reads alignments with a MAPQ score ≥30, and sorted using samtools (Li et al., 2009; Li, 2011). Subsequently, pileup files were generated for each X. maculatus sample and variant calling was processed by BCFtools and VarScan for SNP and Insertion/Deletion 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 (Koboldt et al., 2009). Only variants that were identified by both pipelines were forwarded for further analyses.

To localize fixed variants between the two inbred X. maculatus lines, homozygous loci of JpA were compared to those of JpB. Such loci were identified if ≥75% of JpA is homozygous for one allele, and ≥75% of the JpB is homozygous for another allele. These fixed genetic variants were functionally annotated using snpEff. A genome database was created using the X. maculatus genome sequence and annotation files. Each variant was queried to the genome database to determine if it was located in a genic or intergenic region, and to determine what effect each variant may have on the peptide sequence structure. We focused on genetic variants that led to conservative in-frame insertion/deletion, disruptive in-frame insertion/deletion, frame-shift, missense, start codon change, stop codon loss/gain, splicing pattern/exon usage alteration for further biological function and pathway analyses.

For higher stringency in detecting genetic variants, RNA-Seq reads from skin and liver were also mapped to the X. maculatus genome using Tophat2 (Kim et al., 2013), followed by sequence variant detection using the same method described above. Only variants that were supported by both DNA and RNA sequencing were kept as inter-strain variants for further analyses. Therefore, all genetic variants were supported by both genome re-sequencing and transcriptome sequencing data.

Animal records for scheme A and B interspecies backcross hybrids were collected from the Xiphophorus Genetic Stock Center. For both crossing schemes, backcross interspecies hybrids exhibiting xmrk linked macromelanophore pigmentation patterns were recorded. From November 2008 to December 2019, there were a total of 290 scheme A hybrids at the XGSC that inherited xmrk-Sd (123 exhibited spontaneous tumorigenesis). From January 2010 to June 2019, there were a total of 203 scheme B hybrids that inherited xmrk-Sp (11 of them developed spontaneous tumors). A contingency table was tested using animal numbers for each cross scheme that developed tumors, and only exhibited benign pigment cell hyperplasia. A Chi-square test was used to test if the tumor incidences of both cross schemes are independent.

Processed transcriptomic sequencing reads were mapped to the X. maculatus genome version 5.0 using Tophat2 (Kim et al., 2013), and gene expression of gene models annotated by NCBI was quantified using FeatureCount (Liao et al., 2014). Differentially expressed genes (DEG) between the 25 JpA and 30 JpB were identified using R/Bioconductor package edgeR (Robinson et al., 2010). The DEG test was performed between X. maculatus JpA and JpB basal level gene expression. Genes with Log₂Fold Change (Log₂FC) ≥1, or ≤−1, False Discovery Rate (FDR) <0.05, with receiver operating characteristic (ROC) curve area under curve (auc) ≥0.8 were determined to be DEGs in skin between the two fish lines. The test was performed using JpA as the control to calculate relative gene expression between the two strains; therefore, if Log₂FC ≥1, JpB showed higher-expression of a particular gene, or, if Log₂FC ≤−1, JpA showed higher-expression of a gene.

Sequence homologies between Xiphophorus and human were identified using blastn (Shen et al., 2013). A Reciprocal Best Hit (RBH) method was used to identify human orthologs of Xiphophorus genes. DEGs between X. maculatus JpA and JpB skin, and strain-specific alleles were converted to human homologs and were further submitted to Ingenuity Pathway Analysis (IPA, Qiagen, Redwood City, CA) for functional analyses. Bioinformatics analysis was performed using IPA for clustering and assessing the biological function of DEGs. Herein, the term “pathways” is short for canonical pathways as assigned by IPA based on input genes. Pathway analysis was performed by testing the over-representation of genes belonging to a certain pathway in the input gene list using Fisher exact test. Pathways with an enrichment -log₁₀(p-value) score >3 (p-value < 0.001) were kept for further analysis.

Principle Component Analyses (PCA) was performed using R function prcomp. Library size normalized gene expression values were further scaled and subsequently used for PCA. Plot of the first two dimensions were made using custom R scripts.

Heatmaps, dot plots and chromosomal plots were made using custom R scripts. All scripts are available upon request.

Spontaneous tumor incidences of two interspecies crosses established between X. hellerii and X. maculatus JpA or JpB [i.e., X. hellerii x (JpA × X. hellerii), or scheme A; X. hellerii× (JpB × X. hellerii), or scheme B] strains were first compared. From November 2008 to December 2019, there were a total of 290 scheme A hybrids at the XGSC that inherited xmrk-Sd, and 123 of them exhibited spontaneous tumorigenesis (tumor incidence = 42.4%). From January 2010 to June 2019, there were a total of 203 scheme B hybrids that inherited xmrk-Sp. 11 of them developed spontaneous tumors (tumor incidence = 5.4%; Table 1). Scheme A tumor incidence follows Mendelian segregation of an unlinked xmrk oncogene modifier, R(Diff) (X² = 3.1, df = 1, p-value = 0.08). In contrast, scheme B tumor incidence does not follow Mendelian distribution (X² = 98.4, df = 1, p-value < 2.2 × 10^–16). In addition, the genetic mechanism underlying both cross schemes is statistically independent (X² = 98.5, df = 1, p-value ≤ 2.2 × 10^–16; Table 1).

The JpA and JpB strains were in their 116th and 109th generation of inbreeding, respectively (Figure 1). Using both genome re-sequencing data on 4 JpA and 4 JpB and RNA-Seq data of 25 JpA and 30 JpB samples, we observed 3,292 homozygous polymorphic sites between the JpA and JpB genomes (Supplementary Table 1). The majority of these genetic variants are located within non-coding regions, or located in a coding region, but are not expected to lead to codon changes (i.e., synonymous variants; Supplementary Table 1). However, we identified 393 variants that are expected to produce conservative in-frame insertion/deletion (amino acid change predicated to have minimal effects on protein product), disruptive in-frame insertion/deletion (amino acid change that is predicated to have large effect on protein product), frame-shift, missense, start codon change, stop codon loss/gain, splicing pattern/exon usage alterations, and are therefore predicted to change encoded polypeptide sequences and potentially molecular functions. These 393 variants are located within 244 genes (Figure 2 and Supplementary Table 2). In addition, genetic variants are not randomly distributed throughout the genome, but instead show defined clustering into chromosome “blocks” (Figure 2).

The inter-strain genetic variants are hypothesized to be heterozygous among the ancestral population that gave rise to both strains. These variants, shaped by evolutionary mechanisms (e.g., natural selection, genetic drift and demographic processes), may lead to functional divergence of strain-specific alleles (Wray et al., 2003; Wittkopp, 2006; Wray, 2007; Emerson and Li, 2010). Therefore, we sought to assess transcriptional differences between these two Xiphophorus strains in order to infer functional divergence caused potentially by the genetic variants. When comparing transcriptional profiles between the two strains, 412 genes (237 highly expressed in JpA, 175 highly expressed in JpB) were identified to be differentially expressed in the skin (Log₂FC) ≥1, or ≤−1, FDR < 0.05, ROC curve AUC ≥0.8; Figure 3 and Supplementary Figures 1–3 and Supplementary Table 3).

Using a 3 Kbp length up- or down-stream of the DEGs to estimate cis-effects of polymorphisms of gene expression, we found 163 polymorphisms (5%) are adjacent to 32 DEGs (7.8%) in the skin (Figure 3 and Supplementary Table 4).

Next, we investigated functional differences between JpA and JpB strains using observed genomic and transcriptomic divergence between the two fish strains. Both genetic variants (i.e., genes exhibiting coding sequence differences between JpA and JpB), and DEGs between the two strains, were used together to test over-representation of specific signaling pathways (Supplementary Figure 4; pathway enrichment analyses p-value < 0.001). This strategy was employed to identify pathways that may include genes exhibiting similar expression levels, but differentiated function due to codon changes, and genes that showed diverged expression patterns. Signaling pathway enrichment analyses were performed by comparing each dataset to databases consisting of common signaling pathways and genes.

A total of 3 major pathways (GP6 signaling pathway, Synaptogenesis signaling pathway and Base Excision Repair pathway; Figure 4) were identified (for specific pathways, see Supplementary Table 5). In GP6 signaling pathways, JpA and JpB genetic differences are represented by laminin and multiple collagen genes. Products of these genes play fundamental roles in establishing the extracellular microenvironment, and differences in expression levels and sequences between JpA and JpB may suggest fundamental microenvironmental differences exist. Synaptogenesis signaling pathways are enriched by genes that are involved in Ca²⁺ outflux (cacna2d1), neurotransmitter exocytosis (syn2 and snap25), synapse organization (thbs1 and nign1), and RAC1 mediated effect (e.g., cytoskeletal organization, cell growth, glucose uptake) regulated by efna2, itsn2, kalrn, and farp1 (Ridley, 2006). Another over-represented functional category involves core enzymes of BER genes (i.e., pole, lig1, ogg1), suggesting JpA and JpB exhibited diverged DNA repair functionality or repair efficiency (Figure 4).