Oreochromis niloticus

Using zebrafish as the reference genome, whole-genome alignments of six teleost fishes were generated. The soft-masked genome sequence for zebrafish (Zv9, April 2010) was downloaded from the Ensembl release-75 FTP site. The following soft-masked genome sequences were downloaded from the UCSC Genome Browser: stickleback (gasAcu1, February 2006), fugu (fr3, October 2011), medaka (oryLat2, October 2005), Nile tilapia (oreNil2, February 2012). The H. comes genome sequence (hipCom0) was repeat-masked using WindowMasker (from NCBI BLAST+ package v.2.2.28) with additional parameter “-dust true”. About 32% (158.1/501.6 Mb) of the H. comes genome was masked using this method.
Only chromosome sequences of zebrafish were aligned while unplaced scaffolds were excluded. The reference (zebrafish) genome was split into 21 Mb sequences with 10-kb overlap, while the percomorph fish genomes (H. comes, stickleback, fugu, medaka and Nile tilapia) were split into 10 Mb sequences with no overlap. Pairwise alignments were carried out using Lastz v.1.03.54 (ref. 55 ) with the following parameters: –strand = both–seed = 12of19–notransition–chain–gapped–gap = 400,30–hspthresh = 3000–gappedthresh = 3000–inner = 2000–masking = 50–ydrop = 9400–scores = HoxD55.q–format = axt. Coordinates of split sequences were restored to genome coordinates using an in-house Perl script. The alignments were reduced to single coverage with respect to the reference genome using UCSC Genome Browser tools ‘axtChain’ and ‘chainNet’. Multiple alignments were generated using Multiz.v11.2/roast.v3 (ref. 56 (link)) with the tree topology “(Zv9 (hipCom0 ((fr3 gasAcu1) (oryLat2 oreNil2))))”.
Fourfold degenerate (4D) sites of zebrafish genes (Ensembl release-75) were extracted from the multiple alignments. These 4D sites were used to build a neutral model using PhyloFit in the rphast v.1.5 package^{57 (link)} (general reversible “REV” substitution model). PhastCons was then run in rho-estimation mode on each of the zebrafish chromosomal alignments to obtain a conserved model for each chromosome. These conserved models were averaged into one model using PhyloBoot. Subsequently, conserved elements were predicted in the multiple alignments using PhastCons with the following inputs and parameters: the neutral and conserved models, target coverage of input alignments = 0.3 and average length of conserved sequence = 45 bp. To assess the sensitivity of this approach in identifying functional elements, the PhastCons elements were compared against zebrafish protein-coding genes. Eighty per cent of protein-coding exons (197,508/245,556 exons) were overlapped by a conserved element (minimum coverage 10%), indicating that the identification method was fairly sensitive.
A CNE was considered present in a percomorph genome if it showed coverage of at least 30% with a zebrafish CNE in Multiz alignment. To identify CNEs that could have been missed in the Multiz alignments due to rearrangements in the genomes, or due to partitioning of the CNEs among teleost fish duplicate genes, we searched the zebrafish CNEs against the genome of the percomorph using BLASTN (E < 1 × 10⁻¹⁰; ≥80% identity; ≥30% coverage). Those CNEs that had no significant match in a percomorph genome were considered as missing in that genome. To account for CNEs that might have been missed due to sequencing gaps, we identified gap-free syntenic intervals in zebrafish and the percomorph genomes, and generated a set of CNEs that were missing from these intervals. These CNEs represent a high-confidence set of CNEs missing in the percomorph fishes and thus were used for further analysis. Functional enrichment of genes associated with CNEs was carried out using the GREAT software^{58 (link)} with each CNE assigned to the genes with the nearest transcription start site and within 1 Mb in the zebrafish genome, and significantly enriched functional categories identified based on a hypergeometric test of genomic regions (false discovery rate (FDR) q value < 0.05). We identified the statistically significant gene ontology biological process terms, molecular function terms and zebrafish phenotype descriptions of the genes that are associated with CNEs.
We also predicted CNEs in the Hox clusters of H. comes and other representative teleost fishes using the global alignment program MLAGAN. Orthologous Hox clusters were aligned using MLAGAN with zebrafish as the reference sequence and CNEs were predicted using VISTA.

We obtained 4,122 one-to-one orthologous genes from the gene family analysis (Supplementary Information, section 4.1). The protein sequences of one-to-one orthologous genes were aligned using MUSCLE^{48 (link)} with the default parameters. We then filtered the saturated sites and poorly aligned regions using trimAl (ref. 49 (link)) with the parameters “-gt 0.8 –st 0.001 –cons 60”. After trimming the saturated sites and poorly aligned regions in the concatenated alignment, 2,128,000 amino acids were used for the phylogenomic analysis. The trimmed protein alignments were used as a guide to align corresponding coding sequences (CDSs). The aligned protein and the fourfold degenerate sites in the CDSs were each concatenated into a super gene using an in-house Perl script.
The phylogenomic tree was reconstructed using RAxML version 8.1.19 (ref. 50 (link)) based on concatenated protein sequences. Specifically, we used the PROTGAMMAAUTO parameter to select the optimal amino acid substitution model, specified spotted gar as the outgroup, and evaluated the robustness of the result using 100 bootstraps. To compare the neutral mutation rate of different species, we also generated a phylogeny based on fourfold degenerate sites. The phylogenomic topology was used as input and the “-f e” option in RAxML was used to optimize the branch lengths of the input tree using the alignment of fourfold degenerate sites under the general time reversible (GTR) model as suggested by ModelGenerator version 0.85 (ref. 51 (link)). We calculated the pairwise distances to the outgroup (spotted gar) based on the optimized branch length of the neutral tree using the cophenetic.phylo module in the R-package APE^{52 (link)}. The Bayesian relaxed-molecular clock (BRMC) method, implemented in the MCMCTree program^{53 (link)}, was used to estimate the divergence time between different species. The concatenated CDS of one-to-one orthologous genes and the phylogenomics topology were used as inputs. Two calibration time points based on fossil records, O. latipes–T. nigroviridis (~96.9–150.9 million years ago (Mya)), and D. rerio–G. aculeatus (~149.85–165.2 Mya) (http://www.fossilrecord.net/dateaclade/index.html), were used as constraints in the MCMCTree estimation. Specifically, we used the correlated molecular clock and REV substitution model in our calculation. The MCMC process was run for 5,000,000 steps and sampled every 5,000 steps. MCMCTree suggested that H. comes diverged from the common ancestor of stickleback, Nile tilapia, platyfish, fugu, and medaka approximately 103.8 Mya, which corresponds to the Cretaceous period.

We applied a second approach to identify sex-specific sequences and their location in the genome as described in Akagi et al. (2014) (link) and Böhne et al. (2019) (link) for 10 of the 11 populations. We excluded the Ruzizi River population due to the low number of female samples (Fig. 1). Starting from trimmed sequencing reads (see above), we generated k-mer catalogs per population of all possible k-mers starting with “AG” and a length of 37 bp present in at least five specimens per population using a Python script provided in Akagi et al. (2014) (link). We divided k-mer catalogs into four categories: Y-k-mers = male-specific, Z-k-mers = male-biased, X-k-mers = female-biased, and W-k-mers = female-specific. To this end, we applied a linear regression to the k-mer counts of each population and retained outliers from the general distribution by calculating studentized residuals from a linear model (i.e. jack-knifed residuals). Outliers were defined as all k-mers with an absolute studentized residual value equal to or bigger than 3, as an observation with an absolute value of 3 is deemed to be an outlier (Belsley et al. 1980 ; Hettmansperger 1987 (link); Atkinson 1994 (link)). Subsequently, sex-specific k-mers (i.e. Y- or W-k-mers) were defined as k-mers having zero counts in one sex but not in the opposite sex. Sex-biased k-mers were obtained based on the ratio of counts between males and females, expecting larger counts for the homogametic sex (e.g. X-k-mers = female count/male count > 4, depending on the population analyzed). In summary, we retained outlier k-mers from the linear regression and from there we took (i) sex-specific k-mers (either Y- or W-k-mers) and (ii) sex-biased k-mers with ratios bigger than four for all populations but not in Kalambo River 1, (i.e. ratio threshold set to 12 for Z-k-mers due to the lower number of female samples for this population, see Fig. 1). Next, we tested for an increased amount of sex-specific k-mers per population with a Wilcoxon test, aiming to detect the heterogametic sex of the population. Additionally, we identified k-mers shared among populations in each category with UpSetR (Conway et al. 2017 (link)) in R.
For the Kalambo River (Ka2) and Chitili River (Ch1) populations, we extracted sequencing reads and their mates containing Y-k-mers of each population. Next, we assembled the extracted reads with MEGAHIT (Li et al. 2015 (link)) with –k-max 12. We also placed the resulting contigs onto the Nile tilapia reference genome with BWA and compared the contig data sets using blastX to the NR database in Blast2GO (Gotz et al. 2008 (link)) to retrieve functional annotations.