An adult gold-morph Midas cichlid ( Figure 1A ), bred in the Pearl River Fisheries Research Institute in Guangzhou city of China, was collected. It was determined as a male since its testis was anatomized ( Figure S1 ). Its liver was sampled for Illumina next-generation sequencing (NGS), and its muscle was collected for Nanopore third-generation sequencing (TGS). Several tissues including muscle, heart, and brain were mixed for transcriptome sequencing, which enables a comprehensive coverage of mRNAs for a better gene annotation.

For the NGS, the quality of isolated genomic DNAs was verified by using an integrated strategy of the following methods: (1) DNA degradation and contamination were monitored on 1% agarose gels; (2) DNA concentration was measured by a Qubit DNA Assay Kit in Qubit 3.0 Flurometer (Invitrogen, Carlsbad, CA, USA). A total amount of 0.2 μg DNA per sample was used as the input material for library preparations. The sequencing library was generated using Next Ultra DNA Library Prep Kit (New England Biolabs, Ipswich, MA, USA) following the manufacturer’s recommendations, and index codes were added to each sample for paired-end sequencing in an Illumina HiSeq 2000 system (Illumina, San Diego, CA, USA). In brief, the genomic DNA sample was fragmented by sonication to a size of 350 bp. Sequence artifacts, including reads with adapter sequences, low-quality, and/or unrecognizable nucleotides, were removed. We used Fastp (version 0.19.7) (Chen et al., 2018) to perform basic statistics on the quality of these raw reads. The detailed steps of data processing were set as follows: (1) Discard any paired reads when either one contains adapter sequences; (2) Discard any paired reads when more than 10% of bases are uncertain in either one read; (3) Discard any paired reads if the proportion of low-quality (Phred quality < 5) bases is over 50% in either one read. For the TGS, a 20-kb library was constructed and sequenced on a Nanopore PromethION platform (Nanopore, Oxford, England) using the routine SMRT sequencing strategy.

For the transcriptome sequencing (RNA-seq), the integrity of total RNA was assessed using a Fragment Analyzer 5400 (Agilent Technologies, Santa Clara, CA, USA). Sequencing libraries were generated using NEBNext Ultra RNA Library Prep Kit (New England Biolabs) following the manufacturer’s instructions, and index codes were added to attribute sequences to each sample. Briefly, mRNAs were purified from the total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under an elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5x). First strand cDNAs were synthesized using random hexamer primers and M-MuLV Reverse Transcriptase (RNase H). Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase. After adenylation of 3’ ends of DNA fragments, NEBNext Adaptors with hairpin loop structures were ligated for hybridization. In order to select cDNA fragments of preferentially 250~300 bp in length, the library fragments were purified with an AMPure XP system (Beckman Coulter, Beverly, MA, USA). Subsequently, 3 μl of USER Enzyme (New England Biolabs) was used with size-selected, adaptor-ligated cDNAs at 37°C for 15 min followed by 5 min at 95°C before PCR. Then PCRs were performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers, and Index primer pairs. At last, PCR products were purified, and then their quality was checked on an Agilent Bioanalyzer 2100 (Palo Alto, CA, USA).

After obtaining the sequencing data, we firstly employed KAT (kmer analysis tools, version 2.4.2 (Mapleson et al., 2017)) to evaluate if experimental contamination was hidden in the sequencing library. Then, we used kmc version 3.0.0 (Kokot et al., 2017) to count 17-mer coverage depth and extracted heterozygous 17-mers to estimate the ploidy of the sequenced fish using Smudgeplot default version (Ranallo-Benavidez et al., 2020). Meanwhile, we applied GCE version 1.0.3 (Ranallo-Benavidez et al., 2020) to estimate the genome size and heterozygous rate based on the 17-mer coverage depth.

After ensuring the quality in sequencing and some basic genome features such as ploidy, genome size and heterozygous rate, we generated de novo assemblies using three genome-assemblers including NECAT version 0.0.1_update20200803 (Chen et al., 2021), wtdbg2 version 2.5 (Ruan and Li, 2020), and flye version 2.9-b1768 (Kolmogorov et al., 2019). We chose the best result with a high contiguity (represented by a high contig N50) and a close length to the genome-size estimate as an optimized outcome, and this assembly was used for subsequent analyses. The errors in the assembly generated from TGS were fixed by two rounds of self-correction by TGS reads, then further polished by NGS reads with three iterations using NextPolish version 1.4.0 (Hu et al., 2020).

The quality of our assembly was evaluated by several approaches. Firstly, we applied conserved single-copy genes of Actinopterygii to align our assembly, which is conducted in BUSCOs (version 3) (Seppey et al., 2019). Secondly, we used KAT to compare the k-mers that appeared in both our assembly and NGS; when a high rate of k-mers appeared in the assembly, it is proven to have a high completeness. Thirdly, we applied LAI (Ou et al., 2018) to evaluate the continuity of our assembly based on the ratio of whole LTR retrotransposons (LTR-RTs) using LTR_retriever version 2.9.0 (Ou and Jiang, 2018). Finally, we mapped previously published 37 transcriptomes from diverse tissues (see more details in Table S1 ) to our assembly for a quality evaluation. Meanwhile, we evaluated the two available versions of Midas cichlids genome in NCBI by comparing the indexes of contig N10, N50, N90, and BUSCOs under genome mode and LAI.

We firstly employed RepeatModeler version 2.0.1 (Flynn et al., 2020) calling RepeatScout version 1.0.6 (Price et al., 2005), TRF version 4.09 (Benson, 1999), RECON version 1.08 (Benson, 1999), LtrHarvest version 1.6.2 (Ellinghaus et al., 2008), and Ltr_retriever (Ou and Jiang, 2018) to build a de novo repeat library. Secondly, we searched the repeats belonging to lineages of Cichliformes including Midas cichlids from Repbase database (Bao et al., 2015) as an known library. The two libraries were combined and then searched in the assembled Midas cichlid genome with optimized parameters “-nolow -gff -poly -a -inv -e rmblast” to obtain a complete list of repeats in RepeatMasker version 4.1.1 (Chen, 2004).

We applied three distinct strategies to locate protein-coding regions in the assembled Midas cichlid genome, including de novo, homolog-based and transcripts-based predictions. Two programs, Augustus version 3.3.3 (Stanke et al., 2006) and Snap version 2006-07-28 (Korf, 2004), were used for de novo prediction. The proteins from seven representative species, including six Cichliformes and zebrafish ( Table S2 ), were used for homolog-based prediction that was performed by Tblastn version 2.11.0 (Gish and States, 1993) and Exonerate version 2.2.0 (Slater and Birney, 2005). For the transcripts-based prediction, we firstly generated a de novo assembly of transcriptome using Trinity version 2.11.0 (Haas et al., 2013), and then the produced raw transcripts were mapped onto the assembled genome by Blat version 36 (Kent, 2002) and Gmap version 2017-11-15 (Wu and Watanabe, 2005), followed by verification and integration to assemble spliced alignments (PASA, (Haas et al., 2008)) pipeline. Finally, we employed EvidenceModeler version 1.1.1 (Haas et al., 2008) to integrate these results generated from the above three strategies, and the false predictions with stop codon(s) inside coding regions were removed. To understand the biological functions of deduced proteins, we aligned their sequences to those that were reviewed in Uniprot database using Tblastn.

To explore gene transcriptional difference between the dark and gold phenotypes of Midas cichlid, we downloaded available scale transcriptomes from both groups to perform a comparative transcriptome analysis; the dark- and gold-morph adult individuals were chosen, and four repeated samples were checked for both color types (See Table S3 ). All the transcriptomes were mapped to our Midas cichlid genome assembly by Hisat2 version 2.2.1 (Kim et al., 2019) to obtain coverage depth of every gene. Then, reads mapping counts were calculated based on gene annotation by Stringtie version 2.1.4 (Pertea et al., 2015), and gene transcription levels were represented by TPM (transcripts per kilobase per million mapped reads). DEG analysis between dark- and gold-morphs was conducted using “edgeR” (version 3.15) in the R package (Robinson et al., 2010), with 0.05 as the FDR (false discovery rate) cut-off and a log2FC (fold change) cut-off of 1. Down- and upregulated DEGs were further summarized in a heatmap using TBtools version 1.098725 (Chen et al., 2020a) and enriched by GO (Gene Ontology) analysis using GOATOOLs version 1.2.3 (Klopfenstein et al., 2018).

To explore functional networks of these DEGs, we imported the proteins encoded by DEGs into STRING (with zebrafish as the representative species), a search tool for the retrieval of interacting genes (http://string-db.org) (Szklarczyk et al., 2010), to reveal their protein-protein interactions. The significantly enriched terms (p-value > 0.05) about KEGG pathways or GO biological processes were used to annotate the main functions of these DEGs. The predicted PPI networks were further imported into Cytoscape version 3.9.1 (Shannon et al., 2003) for visualization, and the inside key functional modules were determined by the molecular complex detection technology (MCODE; version 1.32 (Bader and Hogue, 2003)) plug-in of Cytoscape with the following parameters “K-score = 2, degree cutoff=2, max depth = 100, score cutoff = 0.2 “. From each of the key functional modules, we identified hub genes by using the cytoHubba version 0.1 (Chin et al., 2014) plug-in of Cytoscape, and seven common algorithms (MCC, MNC, Degree, Closeness, Radiality, Stress, and EPC) were applied for practical prediction of hub genes.

After quality control, a total of 71.9 Gb of Illumina NGS reads and 77 Gb of Nanopore TGS reads for Midas cichlid were obtained. By breaking NGS data into k-mers (length = 17, 17-mers), a matrix containing distinct 17-mer counts at varying frequency and GC content was produced. We made a density plot that highlights genuine frequency between 40-85 folds with a GC content spreading from 3 to 17 ( Figure 1B ). No unexpected contents beyond the main density blocks were observed, indicating no contamination in the sequencing library of Midas cichlid.

We estimated the sample of Midas cichlid ( Figure 1A ) as a diploid by extracting and analyzing heterozygous 17-mer pairs ( Figure 1C ). Subsequently, we applied heterozygous mode to estimate heterozygous rate that is about 0.23% ( Figure 1D ). By using the routine k-mer frequency distribution analysis (Li et al., 2010a), we estimate the genome size of Midas cichlid to be about 957 Mb ( Figure 1D ).

After comparing genome size and contig N50 value of various genome assemblies generated from different assemblers (NECAT, wtdbg2, flye), we chose the assembly from NECAT as the best one. After polishing sequencing errors produced in TGS, the contig N50 was improved, reaching 10.6 Mb, which is much higher than any version of previously reported draft ( Table 1 ). Based on the final assembly size of 933.5 Mb, we calculated the NGS and TGS data used in this study as about 77- and 82.5-fold depth, respectively. The total sequencing coverage depth is about 159.5-fold.

The estimated assembly completeness from k-mer (k=31, 31-mer) spectrum using kmer analysis tools (KAT) is 99.87% ( Figure 2A ). A Benchmarking universal single-copy orthologs (BUSCOs) estimate indicated a 98.2% completeness, which is slightly higher than the two previous versions of GCA_000751415.1 and GCA_013435755.1 ( Figure 2B ). When we referred to LTR assembly index (the ratio of complete LTR-RTs, LAI), the LAI estimation was failed in the GCA_000751415.1 as complete LTR-RTs are not enough for the analysis. The LAI was evaluated as 6.4 in the GCA_013435755.1. However, it was estimated as 10.6 for our assembly (see Table 1 ), reaching a reference genome level according to the previously established LAI standards (Ou et al., 2018). In addition, the mapping ratio of 35 transcriptomes from different samples and replicates ranged from 94% to 98% (except for one at 87.45%; see Figure 2C ), suggesting a high completeness of our current assembly. Taken together, these results support achievement of a high-quality genome assembly of Midas cichlid in this study.

For genome elements, we totally annotated 41.03% repeats in the assembled genome of Midas cichlid. Among these repeats, DNA transposons were the most (18.2%) in abundance, Retroelements occupied the second place with 12.54%, and the rest are unknown (see more details about the landscape of repeats in Figure S1 ). For gene structures, we totally annotated 25,911 protein-coding genes after removal of possible false prediction such as premature stop codons within the preliminary list. The average length of these genes was 1,684.7 ( Table 1 ). We assessed the quality of the final gene set and corresponding deduced proteins using BUSCOs under protein and transcript modes, and their completeness was estimated to be 90.6% and 90%, respectively. These results suggest establishment of a high-quality list of gene annotation, although it is unable to perform a quality comparison between our assembly and the previous versions (GCA_000751415.1 and GCA_013435755.1) due to non-disclosure of gene annotations in the latter datasets.

Based on our gene annotation and mapping comparison between dark- and gold-morph transcriptomes, we obtained the transcription level of each gene that was represented by transcripts per kilobase million (TPM). A total of 97 upregulated and 180 downregulated DEGs were obtained ( Figures 3A, B ). Interestingly, these downregulated DEGs were enriched in many biological processes (such as melanin metabolic process and melanin biosynthetic process), cellular components (such as pigment granule membrane, melanosome membrane, keratin filament, and intermediate filament), and molecular functions (such as serine-type peptidase activity and serine-type endopeptidase activity; Figure 3C ). However, no significant GO terms were enriched for the upregulated DEGs.

After construction of protein-protein interaction (PPI) networks using those proteins derived from DEGs ( Figures 4A, C ), we identified one ( Figure 4B ) and two ( Figures 4D, E ) key functional modules for up- and downregulated DEGs, respectively. The PPI network derived from upregulated DEGs was enriched into a GO term of “regulation of transcription from RNA polymerase II promoter” and three KEGG (kyoto encyclopedia of genes and genomes) pathways including mitogen-activated protein kinases (MAPK) signaling, Toll-like receptor signaling, and gonadotropin-releasing hormone pathways. The inside key functional module contains eight genes, including fosb, jund, cebpb, jun, junbb, fos, junba, and egr1; among them, junbb, fos, junba, and egr1 were further identified as hub genes ( Figure 4B ).

Similarly, the first key functional module derived from downregulated DEGs ( Figure 4D ) contains 14 genes, including tnni2b.2, tnnt2a, mylpfa, myl1, ckma, ckmb, tnnt3b, pgam2, aldoaa, mybphb, oca2, atp2a1l, actc1a, and myh11a; among them, mylpfa, myl1, ckma, ckmb, tnnt3b, pgam2, and mybphb are identified as hub genes. The second key functional module derived from downregulated DEGs ( Figure 4E ) contains eight genes, including tyrp1b, oca2, pmela, mybphb, tyr, slc24a5, pnp4a, and gpnmb; among them, tyrp1b, oca2, pmela, tyr, and slc24a5 are identified as hub genes.