Skip to main content
Download PDF

Reference genome provide insights into sex determination of silver aworana (Osteoglossum bicirrhosum)

BMC Biology volume 23, Article number: 29 (2025) Cite this article

Abstract

Background

Silver arowana (Osteoglossum bicirrhosum) is a basal fish species with sexual monomorphism, while its sex determination mechanism has been poorly understood, posing a significant challenge to its captive breeding efforts.

Results

We constructed two high-quality chromosome-level genome assemblies for both female and male silver arowana, with scaffold N50 values over 10 Mb. Combining re-sequencing data of 109 individuals, we identified a female-specific region, which was localized in a non-coding region, i.e., around 26-kb upstream of foxl2 gene (encoding forkhead box L2). Its strong interaction with the neighboring foxl2 on the same chromosome suggests foxl2 as a candidate sex-related gene in silver arowana. We subsequently propose a complex gene network in the sex determination process of silver arowana, with foxl2 acting as the central contributor. Transcriptome sequencing of gonads support our hypothesis that the regulation of foxl2 can be influenced by the spatial proximity of the female-specific fragment, thereby promoting ovarian function or inhibiting testicular function to stimulate gonadal differentiation. Furthermore, we found the sex chromosomes to be homomorphic with a potentially recent origin, as a linkage disequilibrium analysis proved minor recombination suppression.

Conclusions

These results taken together serve as a crucial foundation for conducting extensive investigations on the evolution and differentiation of sex-determining mechanisms, as well as the emergence and development of sex chromosomes in various fishes.

Background

Myriad of sex determination mechanisms evolved independently across the tree of life. In particular, the processes are exceptionally diverse among various teleost fishes, including hermaphroditism, environmental sex determination, and genetic sex determination [1, 2]. The genetic sex determination mechanisms of fishes are rather distinct from those of other vertebrates. While heteromorphic sex chromosomes (due to recombination suppression and genetic degeneration in sex-related regions) and common master sex-determining genes (SDGs) are typically identified in birds and mammals [3, 4], sex chromosomes of various fish species are generally different, and master SDGs vary even among closely related species [5, 6]. This indicates that the evolution of new sex-determining regions through “turnover events” might be frequent in fishes [7]; in fact, the sex chromosomes of many teleost species have recently evolved, and some are more resistant to recombination suppression than mammals and birds, perhaps by rarely undergoing chromosome rearrangements or due to slow accumulation of repetitive sequences [8]. Therefore, the practical identification of sex chromosomes or sex-determination regions in fish genomes has been a great challenge.

Some fish species, without heteromorphic sex chromosomes, may have sex-specific regions or express genes specific to determine their genetic sex [9]. Some species have proved that sex determination genes were located in sex-specific region, such as in Atlantic herring. Although diverse sex-determination processes have been reported in certain fishes, most identified SDGs thus far fall into only two groups, owing to the relatively conserved genetic network regulating gonadal development in various vertebrates. Fish SDGs are either involved in the transforming of growth factor-β signaling pathway (e.g., anti-Müllerian hormone and gonadal soma derived factor), or related to some transcription factors (such as dmrt1 and sox3) [10,11,12,13,14,15]. Typically, mutations in the coding regions of SDGs result in the promotion or suppression of gonadal development of specific sex. A previous study reported that although sex chromosome karyotypes are often largely similar among different fish species, it is difficult to identify genes related to sex determination and these genes were only identified in a few model species [16]. Nonetheless, recent advances in whole-genome sequencing have sped up identification of sex-linked or sex-determining regions in some non-model fish species, shedding light on the structural properties of sex chromosomes and their evolution [17,18,19,20].

Osteoglossum bicirrhosum, commonly known as silver arowana, is originally endemic to the Amazonian floodplains. It belongs to the family Osteoglossidae under the order Osteoglossiformes, which is one of the basal lineages of ray-finned fish. The silver arowana is highly exploited as a food source, and as a highly esteemed ornamental species due to its large size, graceful movements, and prehistoric appearance [21]. Like other members of the Osteoglossidae, silver arowana exhibits long lifespan, late sexual maturity, low ovulation rate, mouth-brooding, and a single functional gonad on the left side [22]. In addition, silver arowana has been a good model for studying the detailed mechanisms of growth, sex differentiation, and sex determination within the Osteoglossids due to the relatively mature propagation techniques developed for it. However, it is still difficult to determine the accurate sex of any silver arowana individual because of its sexual monomorphism (Fig. 1A). The difficulty of morphological sex identification for the pairing of reproducing individuals has been a major burden for effective captive breeding of this species. Yet, studies on the sex determination mechanisms of Osteoglossid have been very limited. Among Osteoglossid, the sex chromosome system has only been known for Asian arowana Scleropages formosus (ZW/ZZ) and arapaima Arapaima gigas (XX/XY) [23, 24]. However, sex-linked regions have not been discovered in these species. Such knowledge gaps not only hamper the practical aquaculture of these commercially important species, but also restrict our understanding of the detailed evolution of sex chromosomes and sex determination mechanisms in these fish species.

Fig. 1
figure 1

Structure and evolution of the silver arowana genome. A Photo of a silver arowana used in this study. B A Circos plot of the genome-wide chromosome alignments. The tracks from outside to inside represent the position of chromosomes, GC content, and gene density, respectively. The center of these circles is the chromosome karyotyping of silver arowana. C,D Hi-C heatmap views of female and male individuals, respectively. E Syntenic relationships between silver arowana (upper) and two other arowana species, Asian arowana (S. formosus; center), and African arowana (H. niloticus; lower) with coding genes

Full size image

To improve our understanding of the genome and sex-determination mechanisms of silver arowana, we constructed sex-differential chromosome-level genome assemblies by integration of Illumina short-read, PacBio long-read, and Hi-C (High-throughput Chromosome Conformation Capture) sequencing techniques, and then conducted a genome-wide association study (GWAS; using sex as the target phenotype) on genome re-sequencing data. We identified a female-specific region within a non-coding region and discovered its neighboring coding gene foxl2, which exhibits a high frequency of long-distance interactions with this female-specific region. Therefore, the regulation of foxl2 was predicted to be influenced by the spatial proximity of the female-specific region on the same chromosome, thereby facilitating ovarian development in females and offering insights into a candidate genetic sex-determination mechanism of silver arowana. These genomic resources of this basal fish will become fundamental for in-depth research on the evolution and diversification of sex-determining systems and the origin and evolution of sex chromosomes in various fish species.

Results

Summary of the chromosome-level genome assemblies

PacBio sequencing generated 75.6 Gb and 91.7 Gb of clean reads for the female and male silver arowana genomes with about 97-fold and 118-fold sequencing depth, respectively (Table 1). The assemblies of female and male silver arowana contained a total length of 839.9 Mb and 792.5 Mb (Fig. 1B), and included 2212 and 1620 scaffolds with an N50 of 10.6 Mb and 10.4 Mb respectively (Table 2). The Hi-C sequencing generated 173.4 Gb and 62.0 Gb of clean reads for anchoring both female and male assemblies. The final haplotypic assemblies of 28 chromosomes have a total length of 739.8 and 700.4 Mb, respectively (Fig. 1C, D). These chromosomes have a range of chromosome lengths from 13.6 to 40.4 Mb (Additional file 1: Table S1). The completeness of female and male assemblies based on the BUSCO of Actinopterygii taxonomy were predicted to be 96.70 and 96.21% (Additional file 1: Table S2).

Table 1 Statistics of PacBio sequencing for silver arowana genomes
Full size table
Table 2 Statistics of both genome assemblies for silver arowana
Full size table

A total of 30,221 and 29,199 protein-coding genes were predicted in the female and male silver arowana genomes, with an average of 9.1 and 9.4 exons per gene. The average length of gene, coding sequence (CDS) and introns were 13,824.8 bp and 14,818.0 bp, 1784.3 bp and 1880.0 bp, and 1382.5 bp and 1442.4 bp, respectively. Finally, 27,242 (female) and 25,652 (male) genes were predicted to be functional, accounting for 90.1 and 87.9% of the total genes.

Karyotyping and genomic identification of centromeric satellites

Consistent with the chromosome-level genome assembly (Fig. 1B) and a previous study [25], our cytogenetic analysis confirmed the total number of diploid chromosomes to be 2n = 56 (center of the Fig. 1B). We observed one pair of submetacentric (sm) chromosomes, 10 pairs of subtelocentric (st) chromosomes, and 17 pairs of telocentric (t) chromosomes. This observation is slightly different from some previous reports. For example, Cioffi et al. (2019) found all chromosomes to be acrocentric (a) chromosomes [25], while 4 st + 52 a chromosomes were determined in Suzuki’s study [26]. However, all cytogenetic analyses indicated a lack of metacentric (m) chromosomes in silver arowana, which is in contrast with other Osteoglossiformes species. For instance, black arowana (Osteoglossum ferreirai, in the same genus) has 2 m + 54 st-a, Asian arowana (Scleropages formosus, in the same family) has 8 m-sm + 42 st-a chromosomes, and African Arowana (Heterotis niloticus, in the same order) has 40 m-sm chromosomes [27].

Genome evolution and phylogenetic analysis of silver arowana

To explore structural evolution in chromosomes (Chr), we conducted a collinearity analysis among the genomes of silver arowana, Asian arowanan, and African arowana. Arapaima was excluded from this analysis since its reported genome has not been assembled to a chromosome level. The silver arowana genome shared more collinear fragments with the Asian arowana genome (2n = 50; Fig. 1E), consistent with their closer evolutionary relationship. Chr1, 4, 8, 9, 15, 20, and 21 of silver arowana corresponded to Chr6, 11, 8, 12, 21, 18, and 24 of Asian arowana, respectively. Meanwhile, we observed at least 13 chromosomal fusion and fission events between both arowana species, which is in accordance with a previous report of multiple chromosomal fusions during the evolution of these arowana species [25].

We then constructed a phylogenomic tree of 41 Actinopterigian species using 252 single-copy orthologs (Additional file 2: Table S3). The topology (Fig. 2) was generally consistent with previous reports [28]. These results revealed that Osteoglossomorpha and other species diverged at about 521 million years ago (Mya). Subsequently, Osteoglossoidei and Notopteroidei diverged about 387 Mya. Particularly, within the Osteoglossoidei, the divergence time of silver arowana and Asian arowana occurred about 117 Mya, and both diverged from arapaima (Arapima gigas) and African arowana (Heterotis niloticus) around 205 Mya.

Fig. 2
figure 2

Phylogenetic analysis and divergence time estimation of silver arowana with other representative vertebrates. Numbers on the branches represent the estimated divergence times in millions of years ago (Mya), and the red circles mark the bootstrap support. More details about the genome resources of other species are provided in Additional file 8

Full size image

Screening of structural variations within the sex-specific region between male and female individuals

The female genome was aligned onto the male counterpart to identify sex differences in chromosome structures. Aligned regions covered over 90% of the total chromosome length of both individuals. Seven potential chromosomal inversions were detected after all-against-all alignments (Additional file 3: Fig. S1). These inversions occurred on the terminal regions of Chr1, 4, 7, 9, 12, and 22 of the female chromosomes, which corresponds to Chr1M, 4M, 7M, 9M, 12M, and 22M of the male individual, respectively.

The sex-specific region of silver arowana was determined by conducting GWAS upon analyzing re-sequencing data from 51 females and 58 males. In fact, we successfully identified 998,933 SNP loci, which were further utilized in subsequent investigations to identify markers specific for each sex. The Manhattan plots (Fig. 3A) revealed a distinct segment spanning 4 kb at the position 28.0 ~ 28.4 Mb on Chr 28, which exhibited complete differentiation between both sexes with a p-value of 10−4 across the entire genomes. Notably, all sex-related SNPs surpassing the threshold were found to be concentrated within this specific genomic region (Fig. 3B). Nine sex-specific SNPs with a p value < 10−15 and a female-specific region were identified (Fig. 3C).

Fig. 3
figure 3

A genome-wide association study identifies SNPs associated with phenotypic sex in silver arowana. A A Manhattan plot of genome-wide sex associations in 109 individuals (including 51 females and 58 males). Two molecular markers were identified at the position 28.20 Mb on Chr28 (ChrZ/W) with 100% compatibility with sex and a p-value of 1e−4. B An enlarged plot of the Chr28 (ChrZ/W) from 28.0 to 28.4 Mb. The x-axis represents the SNP position on Chr28. The y-axis shows the significance of association (− log10 (p) value) in the examined population. Nine genes are localized within the window, with the direction of transcription (upstream/downstream) being annotated with arrows. C Nine SNPs associated with sex were identified, and a 511-bp female-specific region was discovered on the ChrW, exhibiting a p-value lower than 1e − 15. The homologous segments of the female-specific region on the ChrZ are only 35 bp in length. The nearest coding gene, foxl2, is located downstream of the female-specific region at a genomic distance of 26 kb

Full size image

We observed that females exhibit heterozygosity at these 9 loci while males display homozygosity, suggesting that silver arowana possesses female heterogamety with ZW sex determination. Consequently, the Chr28 of silver arowana can be split into ChrZ and ChrW. Our subsequent analysis of individually sequenced silver arowana samples revealed that male individuals exhibited no coverage in the regions specific to females. The female-specific region of ChrW spans approximately 511 bp; however, the homologous segments on the ChrZ are only 35 bp in length. These segments lack any coding gene, with the nearest neighboring gene being foxl2 located at 26-kb downstream.

Chromatin accessibility and 3D genomic interactions

We conducted a comparative study of the 3D genomes using silver arowana ovary and testis as experimental materials, and compared the Hi-C data obtained from each sample. We constructed the Hi-C complex network for Chr28 (Fig. 4A). The distribution of interactions between the ovaries and testes reveals distinct variations in interaction densities across different regions, thereby delineating topologically associating domains (TADs) that demarcate discrete regions within the chromosome. We enlarged the chromatin conformation map of the TAD containing the female-specific region and its neighboring regions, revealing a significantly higher interaction strength within this TAD for the ovary compared to the other part of the genome (Fig. 4B).

Fig. 4
figure 4

High-resolution chromatin interaction map of the ZW chromosome region (27.24–28.75 Mb) in ovarian and testicular tissues. A The KR-normalized contact matrix at 10-kb resolution. B The normalized interaction frequencies for individual Hi-C libraries from female and male samples. C Graphical representation of chromatin in ovary (left) and testis (right)

Full size image

To investigate the potential interactions between the female-specific region and its upstream and downstream coding genes (Fig. 3B), as well as to compare these interactions in ovaries versus testes, we utilized the ICE normalization method [29] to perform a statistical analysis of interaction frequencies within 10-kb genomic windows. There was a significant interaction between the female-specific region sequence and foxl2, as well as with pik3cb and upstream gene copb2. Our results demonstrated that the ICE-normalized Hi-C data between the female-specific region and foxl2 were 38.8 in the ovary, compared to 12.5 in the testis. Based on the complete sex chromosomes, a T-test was also conducted using the ICE normalization values with difference multiples between ovaries and testes. The test focused on the foxl2 gene region and the female-specific region. Parameters for the test included two-sided tests with a confidence level set at 95%. The calculated p-value was 4.21e − 04, indicating a statistically significant difference in interaction intensity between ovaries and testes within this genomic region. The histogram in Fig. 4B showed the proportion of the foxl2 reads that interact with the female-specific regions relative to the total reads within the interaction regions determined by Hi-C sequencing. Specifically, in the ovaries containing the female-specific regions and the foxl2 gene, this proportion was 1.62%; while in the testes, this proportion was only 0.21% (p-value was 1e − 05). The ICE-normalized Hi-C data between the female-specific region and pik3cb were 4.43 in the ovary and 3.30 in the testis (p = 0.42), and between the female-specific region and copb2 were 17.50 in the ovary and 10.03 in the testis (p = 0.12). Thus, no significant differences were observed in the interactions of the female-specific region with pik3cb or copb2 between the ovary and testis.

Figure 5 illustrates the proportional distribution of reads that interact with foxl2 within a 100-kb region upstream and downstream of the foxl2 gene, using an observation window of 1 kb. The y-axis labeled “Ratio” represents the proportion of reads interacting with foxl2 in this region relative to the total number of reads. Interestingly, our findings reveal that in females, foxl2 exhibits a significant degree of self-interaction and demonstrates robust interactions with the female-specific fragments along with their adjacent sequences. Conversely, in males, foxl2 primarily engages in self-interaction while displaying minimal or negligible interactions with other fragments within the 100-kb proximity.

Fig. 5
figure 5

Interaction between foxl2 and the neighboring regions. The interaction strength refers to the ratio of sequences that produce mutual interactions

Full size image

The interaction between the female-specific region and foxl2 suggests that despite a linear separation of 26 kb, the physical distance between these two locations on the chromosome is closer. The absence of a significant interaction between these two positions on the ChrZ implies that the chromosomal structure exhibits a more relaxed configuration at this location (Fig. 4C).

Regulatory role of foxl2 in gender differentiation of silver arowana

We collected nucleotide sequences of 74 sex-related genes from 41 ray-finned fish species [30] and identified 36 of them in the silver arowana genomes, of which 30 were single copy (Additional file 4: Table S4). The overall copy number, gene expansion, and loss pattern showed considerable diversity among the examined fishes (Fig. 6). Gene copy number contributes to genetic diversity both between and within species [31, 32]. This variation can create extreme variance in fitness and promote adaptive evolution [33]. However, no marked expansion was observed in these sex-related genes within the silver arowana genome, and we did not detect any sex-specific SNPs within these genes.

Fig. 6
figure 6

A heatmap of sex-related genes in various vertebrates. Color boxes represent gene number in each species

Full size image

To further examine the regulatory mechanism of sex determination in silver arowana, we analyzed gene expression in testes and ovaries at three developmental stages (Fig. 7A, Additional file 5: Table S5). Validation by qPCR, the expression level of foxl2 in the ovary was significantly higher than that in the testis. Moreover, during ovary development, there was a significant increase in foxl2 expression; however, no difference was observed during testis development (Fig. 7B).

Fig. 7
figure 7

The critical sex-related candidate gene foxl2 in silver arowana. A Relative expression of sex-related genes at three developmental stages. See more details about classification of each stage in Table S9 of Additional file 9. B qPCR validation of target genes including foxl2. Asterisk (*) represents a significant difference (p < 0.05). C Proposal of the possible signaling pathways for foxl2 in silver arowana gonads. D A phylogenetic tree of foxl2, indicating its relative conservation among various vertebrates. E Alignments of foxl2 and flanking sequence among representative species within the order Osteoglossiformes

Full size image

We found that gene expression patterns at different stages of gonadal development showed certain correlation between each other. We therefore clustered these genes based on their expression patterns. Several highly correlative gene expression patterns were identified, and the correlation were determined by Pearson correlation coefficients (see more details in Additional file 6: Table S6), including the co-expression of some known genes involved in gonad development in various vertebrates, such as emx2, map3k4, and sox3. Notably, we observed that foxl2 transcribed negatively with sox8a, gadd45g, cyp26b1, and dhh. These are key genes for male sex determination in various vertebrates [34, 35]. By contrast, the expression of foxl2 did not strongly correlate or even negatively correlated with several known genes for female sex determination, such as rspo1, esr1, fsta, igf1r, and gfra1. These genes are usually at the upstream of foxl2 in the putative regulatory cascade [36], and hence might not be regulated by foxl2 or the associated transcription factories. According to report, both pdgfa and trpv4 genes demonstrate a regulatory relationship with the sox9 gene, indicating the potential for indirect regulation between foxl2 and these genes [37,38,39].

Based on the correlations of gene expression (Additional file 6: Table S6), and the regulatory relationships between genes as predicted by previous studies [39,40,41,42], we propose a potential gene regulatory network associated with foxl2 and other sex determination-related genes in silver arowana (see more details in Fig. 7C). This network includes 14 key genes, some of which have been reported to be directly or indirectly involved in the sex determination in other vertebrates, such as foxl2, sox3, sox9, dmrt1, esr1, and sox8a. In the dsRNA injection group, the expression level of foxl2 in ovary was significantly lower than that of the control group 24 h post-injection (p < 0.01) (Additional file 3: Fig. S2). Meanwhile, there was no significant change in sox9 expression compared to the control group (p > 0.05), whereas dmrt1, Sox8a, and esr1 exhibited a significant increase in their expression levels (p < 0.05), which is consistent with the observed correlation among these genes (Additional file 3: Fig. S3).

With the downregulation of foxl2 expression, there was a concomitant upregulation in the expression levels of dmrt1, sox8a, and esr1, while no significant changes were observed in sox9 mRNA expression. These findings were consistent with the results obtained from transcriptome sequencing and analyses of gene expression correlation. The predictions of gene regulatory pathways (Fig. 7C) were partially substantiated by this evidence.

Discussion

We performed whole-genome sequencing and resequencing of sexually monomorphic silver arowana, and obtained male and female chromosome-level genome assemblies. A female-specific region in non-coding region with strong interaction with foxl2 was identified. The concentrated distribution of sex-specific SNPs suggests the potential occurrence of recombination suppression in the vicinity of sex-specific markers on the sex chromosomes. In fact, the suppression of recombination is a crucial step in the process of sex chromosome evolution [43]. Recombination suppression in the heterogametic sex chromosome usually leads to an elevation of evolutionary rate and an accumulation of deleterious mutations as the chromosome “mature” [44]. However, silver arowana possesses homomorphic sex chromosomes. In many genetically sex-determined organisms without heterotypic sex chromosomes, the region of recombination suppression on the sex chromosome does not extend significantly beyond the sex-determining locus, possibly due to these young and undifferentiated sex chromosomes [45].

The three-dimensional architecture of chromosomes plays a pivotal role in the regulation of gene expression. At the DNA level, distal regulatory elements can modulate the expression of target genes through spatial proximity. TADs are considered as fundamental units governing chromosome structure. Interactions between regulatory elements and gene promoters predominantly occur within TADs [46]. In Drosophila, fixed elements spanning extensive distances within the same TAD regulate neighboring regulatory elements and target genes, thereby facilitating rapid activation of developmental genes and co-regulation of distant homologous genes [47]. The enhancer Enh13, located in the non-coding region, has been identified as a critical regulator of testis development in mice by controlling the expression of the Sox9 gene over long distances [48]. In humans, three cooperative regulatory elements have been identified that synergistically drive robust expression of the Sox9 gene across long genomic distances, thereby ensuring proper testis development and male sexual differentiation [49]. The utilization of fish non-coding DNA sequences for spatial regulation in sexual determination remains unexplored in most fish species. In our current study, however, we identified a robust interaction between the female-specific region on the ChrW and foxl2 through Hi-C sequencing of silver arowana. However, weak corresponding interaction was observed in male fish due to the absence of a female-specific region on the ChrZ. Meanwhile, we observed that the expression level of foxl2 in the ovary significantly higher than that in the testis (Fig. 7B). Therefore, it is postulated that the presence of the female-specific region induces chromatin folding, thereby facilitating spatial proximity and regulation of foxl2 expression in females. Conversely, in males lacking the female-specific region, the linear arrangement of the chromosome results in lower level of foxl2 expression.

It is well known that the foxl2, encoding a forkhead transcription factor, is a key gene for ovarian development with a highly conserved role in regulation of gonadal differentiation in various vertebrates [50,51,52]. Zebrafish foxl2 synergistically regulates ovary development, and deletion of both genes causes complete female-to-male sex reversal [53]. The knockdown experiment of foxl2 in Nile tilapia (Oreochromis niloticus) also generated a female-to-male reversal [54]. In our present study, we performed both qRT-PCR and transcriptome sequencing on foxl2. Both results indicated that the expression of foxl2 in the female gonads was significantly higher than that in the male silver arowana (p < 0.05) (Fig. 7A, B). Combining with its interaction relationship with the female-specific region, we propose that foxl2 may be one of the sex-determining genes in silver arowana. Additionally, we also observed that foxl2 transcribed negatively with several male sex determination genes, such as sox8a, gadd45g, cyp26b1, and dhh, but it has no strong correlation with the female sex determination genes. A previous study has indicated that the foxl2 promotes ovarian development partly by repressing expression of male-related genes [55]. Therefore, it is reasonable to hypothesize that the female-specific region may promote the foxl2 gene to be highly expressed, which subsequently represses male-related genes and finally leads to a female phenotype at an early stage.

Comparing the genomic sequences around foxl2 and the female-specific region in silver arowana with other species from the order Osteoglossiformes, we revealed that the foxl2 gene was highly conserved in most of the representative species (Fig. 7D), but the upstream region of this gene was more variable in other species (Fig. 7E). It is important to note that such homologous sequence around the sex-specific SNPs was absent in other Osteoglossiformes fishes, even in the close relative Asian arowana (S. formosus; see Fig. 7E). This suggests that the acquisition of these candidate sex-determining regions in silver arowana may have occurred after its divergence from other arowana species, such as Asian arowana and African arowana. In other words, silver arowana owns a potentially young sex chromosome.

Our present work could be substantially enhanced with more investigations. Specifically, a functional experiment should be undertaken to validate the contribution of the female-specific region to the gene function of foxl2. Additionally, an interesting analysis focused on predicting transcription factor binding sites may be imperative for further elucidating the potential functions of this region. We aim to address these limitations in our future research.

Conclusions

In this study, we reported the first chromosome-level genome assemblies for both male and female silver arowanas. Their scaffold N50 values reached over 10 Mb, suggesting that they are of high quality. Our karyotype and genomic analyses indicated that the centromeres of almost all chromosomes are close to the terminal regions (i.e., mostly acrocentric, subtelecentric, or telecentric), a pattern quite different from other members in the same family of Osteoglossidae. Multiple chromosomal fusion and fission events may have occurred during the evolution of this species. Furthermore, a female-specific region was identified to be strictly associated with the phenotypic sex. The female-specific region was localized in the non-coding region at 26-kb upstream of foxl2, and a strong interaction between this female-specific region and the foxl2 region on the same chromosome suggests that foxl2 may be a candidate sex-related gene in silver arowana. Transcriptomic data revealed that the foxl2 may play a predominant role in sex determination and sex differentiation via a cross-regulation with sox8a, cyp11a1, and some other auxiliary genes. In summary, our chromosome-level genome assemblies, multiple tissue transcriptomes, and genome resequencing data provide a valuable genetic resource to improve our understanding of the evolution of sexual system in silver arowana.

Methods

Sample collection and chromosome karyotyping

All procedures for handling and treatment of experimental fishes during this study were approved by the Laboratory Animal Ethics Committee of the Pearl River Fisheries Research Institute (PRFRI), Chinese Academy of Fishery Sciences (CAFS), China. Experiments were carried out according to the Guide for the PRFRI in Guangzhou, China. The number of fish samples used in this study is summarized in Additional file 7: Table S7.

To assess the accurate chromosome number, we conducted a karyotype analysis [56] using kidney tissues of 20 randomly selected 1-year-old silver arowanan individuals. Ten micrograms per gram phytohaemagglutinin (PHA; Yeasen, Shanghai, China) was injected into the abdominal cavities 12 h prior to fish dissection. Fishes were treated with colchicine (4 μg/g; Apexbio, Houston, TX, USA) 60 min before euthanasia and chromosome preparation. The kidney tissues were dissected and grounded in 0.75% NaCl solution, followed by hypotonic treatment using 75.0 mmol/L KCl for 30 min, and fixation in 3:1 methanolacetic acid (three changes) for 30 min. Chromosome slides were stained with Giemsa stain solution (Merck KGaA, Darmstadt, Germany). For each specimen, images of 50 metaphases were collected for calculating chromosome numbers.

Illumina sequencing and genome size estimation

Muscle tissues of one female and one male silver arowana were separately collected for genomic DNA (gDNA) extraction using a Marine Animal DNA Kit (CoWin Biosciences, Taizhou, Jiangsu, China). Gel electrophoresis (1.0% agarose) and NanoDrop 2000 spectrometer (Thermo Scientific, Waltham, MA, USA) were employed to evaluate the quality and quantity of these DNA samples. To perform paired-end sequencing on an Illumina NovaSeq 6000 (Illumina Inc., San Diego, CA, USA), high-quality libraries of gDNAs with the insert size of 350 bp were constructed following the manufacturer’s instructions. Trimmomatic v0.39 [57] was applied to filter adaptors and low-quality reads. The estimation of genome sizes for both male and female specimens was conducted through a kmer frequency analysis (k = 21) using the jellyfish software (v2.2.10; http://www.cbcb.umd.edu/software/jellyfish/index.shtml) as described before [58].

PacBio and Hi-C sequencing

A total of 50 μg of high-quality gDNAs from the muscle tissues of the female or male fish were used to generate a 40-kb SMRTbell library, respectively. Library preparation and sequencing were performed in accordance with the manufacturer’s protocols (Pacific Biosciences, Menlo Park, CA, USA).

To build Hi-C libraries, muscle samples of the female or male fish were fixed using formaldehyde. Cell lysis, chromatin digestion, proximity-ligation treatments, DNA recovery, and subsequent DNA manipulations were performed according to the methods described in a previous study by using MboI as the restriction enzyme [59]. Both Hi-C libraries were paired-end sequenced on an Illumina Novaseq 6000 in accordance with the manufacturer’s protocols.

Genome assembly and evaluation

The draft genomes were assembled using Canu v2.1.1 (https://github.com/marbl/canu) with PacBio sequencing reads of male and female, respectively [60]. Hi-C reads were collected to build Hi-C maps using Juicer v1.6 (https://github.com/aidenlab/juicer) [61]. The assembled contigs were then ordered and oriented using the 3D-DNA pipeline v201008 (https://github.com/aidenlab/3d-dna) [62]. Juicebox Assembly Tools v2.12.00 (https://github.com/aidenlab/Juicebox) was applied [61] to manually identify and remove assembly errors.

To generate Hi-C contact maps, we analyzed the Hi-C sequencing data using the Hic-Pro2.11.3 pipeline with default parameters and LIGATION_SITE = GATCGATC [63]. Valid ligation products were merged for interaction matrix construction. Chromosomal contact map was plotted in HiCPlotter v0.6.6 (https://github.com/akdemirlab/HiCPlotter) based on the interaction matrix.

Completeness of both genome assemblies was evaluated by using BUSCO (Benchmarking Universal Single-Copy Orthologs) v5.2.2 (https://busco.ezlab.org/) [64] with the actinopterygii_odb10 database (https://busco-data.ezlab.org/v5/data/lineages/) as the reference under the genome mode.

Genome annotation

To generate transcriptome reads, total RNAs were extracted by using the TRIzol kit (Invitrogen, Carlsbad, CA, USA) from muscle, gonad, liver, skin, and fins from the male or female fish, and were then pooled in equal amount. These pooled RNAs were reversely transcribed to cDNAs and were subsequently amplified with 14 cycles of PCR by SeqAmp DNA polymerase (TaKaRa, Beijing, China). Full-length cDNAs of different sizes were selected and corresponding cDNA libraries were constructed using BluePippin (Sage Science Inc., Beverly, MA, USA). The size-selected products were then used for a 12-cycle large-scale PCR to generate double-stranded cDNAs for SMRTbell library construction, according to the Iso-Seq protocol (Pacific Biosciences, SMRT Analysis 2.3.0; https://smrt-analysis.readthedocs.io/en/latest/SMRT-Analysis-Software-Installation-v2.3.0/).

Eight SMRT cells were sequenced on a PacBio Sequel platform. The Iso-Seq sequencing data were analyzed with IsoSeq v3.4.0 (https://github.com/PacificBiosciences/IsoSeq). Augustus v.3.3.1 and GlimmerHMM [65] were employed to perform de novo gene prediction. We generated a final set of gene annotation using the software MAKER2 [66]. Coding genes were aligned against the NCBI GenBank NR database using blastx with an e-value of 1e−5 (other parameters were set as the defaults).

RepeatModeler 2.0 [67] was applied to predict repeat families in the silver arowana genomes. Tandem repeats were searched using Tandem Repeat Finder v4.09.1 (https://github.com/Benson-Genomics-Lab/TRF) with the parameters “2 5 7 80 10 50 2000” (for matching weight, mismatching, indel penalty, match indel probability, minimum alignment score and maximum period size to report). Above results were filtered by pyTanFinder (https://github.com/Kirovez/pyTanFinder) to remove redundancy of the repeat families. These predicted repeat families and tandem repeats were combined with the existing fish repeat library from Dfam and RepBase (20,170,127) as the input library for repeat annotation and the program RepeatMasker v4.0.7 [68].

Genome evolution analysis

Genome sequences of silver arowana (generated in this study) and 40 representative fish species (downloaded from the NCBI database; see Additional file 8: Table S8) were retrieved. Single-copy orthologue genes were detected using OrthoFinder v2.5.4 (http://www.stevekellylab.com/software/orthofinder) with deduced protein sequences and default parameters. Finally, a data set of 252 single-copy orthologs was obtained. The nucleotide sequences of these orthologs were individually aligned using MAFFT v7.475 [69]. A phylogenetic tree was constructed using IQtree2 v2.0.5 [70] with edge-linked proportional partition model, 1000 ultrafast bootstrap replicates, and 1000 replicates for SH-like approximate likelihood ratio test.

Identification of sex-specific regions

Synteny identification between the chromosome-scale genome assemblies of male and female silver arowana was performed by Nucmer module of MUMmer v4.0beta1 [71]. Sequences with alignment identity > 0.9 and alignment length < 2 kb were retained. The chromosome synteny regions and structural variations were visualized using RectChr (https://github.com/BGI-shenzhen/RectChr). To confirm the identified structure variations, we employed Minimap2 [72] with default parameters to align PacBio HiFi reads to the female or male assembly correspondingly, and then applied IGV (Integrative Genomics Viewer) [73] to check alignments and to exhibit the detailed read coverage of these critical regions.

To identify sex-determining regions, we conducted whole-genome re-sequencing of 51 female and 58 male individuals. gDNAs were extracted from muscle tissues of these individuals (a total of 109) separately using the Marine Animal DNA Kit (CoWin Biosciences). Corresponding libraries with an insert size of 350 bp were sequenced on an Illumina Novaseq 6000 platform. Sequencing reads were aligned on the male or female chromosome-level genome assembly using bwa (v0.7.17, http://bio-bwa.sourceforge.net/) [74] with default parameters. SNPs were called using GATK (v3.8.0; https://software.broadinstitute.org/gatk/) [75], and then were filtered with the following parameters: variant quality vs. depth of non-ref samples less than 2.0 folds, Phred score Fisher’s test p-value for strand bias greater than 60.0, and mapping quality less than 40.0.

A genome-wide association study (GWAS) using a generalized linear model was conducted in TASSEL (v5.0, https://tassel.bitbucket.io/) [76] to identify genomic regions significantly correlated with sex (MAF < 0.05). Manhattan plots at both genomic and chromosomal levels were generated using a 5-kb window with a mean value of − log10(p). LocusZoom (http://csg.sph.umich.edu/locuszoom/) was employed to produce regional association plots. Sex-specific regions were identified based on comparisons of sex-specific SNPs and genomic alignment results.

Chromatin interaction analysis

Hi–C technology was utilized to perform high-resolution genome-wide chromatin interaction mapping between the gonads of both female and male silver arowana. We performed a HiC-Pro analysis (version 3.1.0, https://github.com/nservant/HiC-Pro) on the Hi-C sequencing data [63]. The Hi-C sequencing utilizes the MobI restriction enzyme with a GATC cut site. The alignment of Hi-C reads to the reference genome assembly was conducted using the bowtie2 method, and then based on read positions a Hi-C interaction matrix was generated at resolutions of 10 and 1 kb. Finally, an interaction heatmap depicting the Hi-C relationships was plotted at a resolution of 10 kb.

To further elucidate the interactions between target genes and their neighboring regions, we conducted an analysis of Hi-C sequencing reads within a 100-kb range upstream and downstream of the target gene. The strength of each interaction was quantified using a statistical window of 1000 bp, wherein we calculated the proportion of reads that interacted with the target gene within each window of 100-kb region. This proportion exhibited a positive correlation with the strength of interaction, which was subsequently visualized through a bar chart.

Gene expression analysis

From three male and three female silver arowana individuals, we collected gonads at three developmental stages (see Table S9 of Additional file 9 for stage classification) [77, 78]. Three biological replicates, each with three technical repeats, were generated. Total RNA was extracted with the TRIzol kit (Invitrogen), and then mRNAs were isolated and purified from each sample with an Oligotex mRNA Midi Kit (Qiagen GmbH, Hilden, Germany). RNA concentration and purity were measured by a Nanodrop spectrophotometer, and the RNA integrity was verified using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Transcriptome libraries were constructed using NEBNext Ultra RNA Library Prep Kit (NEB, Ipswich, MA, USA) according to the manufacturer’s protocol, and then were sequenced on an Illumina Novaseq 6000 platform with the 150-bp paired-end protocol.

Sequencing reads from ovary tissues were aligned on the female genome assembly and sequencing reads from testis tissues were aligned on the male genome assembly by using hisat2 v2.1.0 [79]. Gene expression levels in terms of gene count, transcript count, and FPKM (fragments per kilobase of exon model per million mapped fragments) values were calculated using StringTie v1.3.5 [80]. Differentially expressed genes (DEGs) were identified using EdgeR [81]. The Retrieval of Interacting Genes/Proteins database (STRING v11.0) [82] was searched to construct a protein–protein interaction (PPI) network.

Pearson’s correlation was applied to examine the relationship among gene expression patterns using the R function cor.test. We validated the expression pattern of the target gene foxl2 using quantitative RT-PCR (qRT-PCR). Total RNA was isolated from gonad tissues at the three developmental stages, from each of the three male and three female individuals, using Total RNA Kit I (Omega Bio-Tek, Norcross, GA, USA). A total of 1 µg of RNA from each sample was reversely transcribed by using PrimerScript 1st Strand cDNA Synthesis Kit (TaKaRa, Dalian, Liaoning, China). qRT-PCRs were performed on an ABI QuantStudio6 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA) using TB Green Premix Ex Taq (TaKaRa). Each 20-μL PCR reaction includes 10 μL of TB Green Premix Ex Taq, 1 μM of each primer, and 1 μL of cDNA. GAPDH was used as the internal control reference. Melting curve analyses were performed following the amplifications. Quantification of the abundance of selected mRNA transcripts was conducted based on PCR amplification efficiencies and crossing point (CT) differences [83].

RNA interference

Design-specific primers (Additional file 9: Table S10) for foxl2 dsRNA amplify the target fragment of 290 bp. Subsequently, ligate the amplified fragment into a T7 vector (Vazyme, China) and transfect this construct into competent cells (TOLOBIO, China). Screen PCR products 1 and 2 using primers Foxl2-ds5-F/Foxl2-T7ds5-R and Foxl2-ds5-R/Foxl2-T7ds5-F (Additional file 9: Table S10), respectively, to isolate PCR products containing the sense T7 adapter. Combine PCR products 1 and 2 in a ratio of 1:1, then utilize the T7 RNAi transcription Kit (Vazyme, China) to synthesize the final dsRNA.

The dsRNA was administered intravenously through the caudal vein into the silver arowana’s bodies. The injected dose was 6 μg dsRNA/g body weight (BW). A comparable volume of deionized distilled water (ddH2O) was administered to the control group. The various groups were housed in tanks, where their conditions were monitored, and ovaries from both experimental and control female fish were harvested 24 h post-injection. Changes in the expression levels of related genes dmrt1, sox8, sox9, and esr1 following interference with foxl2 expression were validated using specific primers (Additional file 9: Table S10).

Data availability

Sequence data that support the findings of this study were deposited at NCBI under the Bioproject PRJNA881487 [84] and PRJNA881492 [85].

Abbreviations

SDGs:

Sex-determining genes

Hi-C:

High-throughput chromosome conformation capture

GWAS:

Genome-wide association study

CDS:

Coding sequence

sm:

Submetacentric

st:

Subtelocentric

t:

Telocentric

a:

Acrocentric

m:

Metacentric

Chr:

Chromosomes

Mya:

Million years ago

TAD:

Topologically associating domain

PHA:

Phytohaemagglutinin

gDNA:

Genomic DNA

DEGs:

Differentially expressed genes

PPI:

Protein-protein interaction

References

  1. Godwin J. Neuroendocrinology of sexual plasticity in teleost fishes. Front Neuroendocrinol. 2010;31(2):203–16.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  2. Kobayashi Y, Nagahama Y, Nakamura M. Diversity and plasticity of sex determination and differentiation in fishes. Sex Dev. 2013;7(1–3):115–25.

    Article  CAS  PubMed  Google Scholar 

  3. Nam K, Ellegren H. The chicken (Gallus gallus) Z chromosome contains at least three nonlinear evolutionary strata. Genetics. 2008;180(2):1131–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. Veyrunes F, Waters PD, Miethke P, Rens W, Mcmillan D, Alsop AE, et al. Bird-like sex chromosomes of platypus imply recent origin of mammal sex chromosomes. Genome Res. 2008;18(6):965–73.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Roberts RB, Ser JR, Kocher TD. Sexual conflict resolved by invasion of a novel sex determiner in Lake Malawi cichlid fishes. Science. 2009;326(5955):998–1001.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Ross JA, Urton JR, Boland J, Shapiro MD, Peichel CL. Turnover of sex chromosomes in the stickleback fishes (gasterosteidae). Plos Genet. 2009;5(2):e1000391.

    Article  PubMed  PubMed Central  Google Scholar 

  7. Vicoso B. Molecular and evolutionary dynamics of animal sex-chromosome turnover. Nat Ecol Evol. 2019;3(12):1632–41.

    Article  PubMed  Google Scholar 

  8. Charlesworth D. When and how do sex-linked regions become sex chromosomes? Evolution. 2021;75(3):569–81.

    Article  PubMed  Google Scholar 

  9. Fu B, Zhou Y, Liu H, Yu X, Tong J. Updated genome assembly of bighead carp (Hypophthalmichthys nobilis) and its differences between male and female on genomic, transcriptomic, and methylation level. Front Genet. 2021;12:728177.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Kamiya T, Kai W, Tasumi S, Oka A, Matsunaga T, Mizuno N, et al. A trans-species missense SNP in Amhr2 is associated with sex determination in the tiger pufferfish, Takifugu rubripes (fugu). Plos Genet. 2012;8(7):e1002798.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Matsuda M, Nagahama Y, Shinomiya A, Sato T, Matsuda C, Kobayashi T, et al. DMY is a Y-specific DM-domain gene required for male development in the medaka fish. Nature. 2002;417(6888):559–63.

    Article  CAS  PubMed  Google Scholar 

  12. Myosho T, Otake H, Masuyama H, Matsuda M, Kuroki Y, Fujiyama A, et al. Tracing the emergence of a novel sex-determining gene in medaka, Oryzias luzonensis. Genetics. 2012;191(1):163–70.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Nanda I, Kondo M, Hornung U, Asakawa S, Winkler C, Shimizu A, et al. A duplicated copy of DMRT1 in the sex-determining region of the Y chromosome of the medaka, Oryzias latipes. Proc Natl Acad Sci U S A. 2002;99(18):11778–83.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Smith CA, Roeszler KN, Ohnesorg T, Cummins DM, Farlie PG, Doran TJ, et al. The avian Z-linked gene DMRT1 is required for male sex determination in the chicken. Nature. 2009;461(7261):267–71.

    Article  CAS  PubMed  Google Scholar 

  15. Takehana Y, Matsuda M, Myosho T, Suster ML, Kawakami K, Shin-I T, et al. Co-option of Sox3 as the male-determining factor on the Y chromosome in the fish Oryzias dancena. Nat Commun. 2014;5:4157.

    Article  CAS  PubMed  Google Scholar 

  16. Bai Z, Liu F, Li J, Yue GH. Identification of triploid individuals and clonal lines in Carassius auratus complex using microsatellites. Int J Biol Sci. 2011;7(3):279–85.

    Article  PubMed  PubMed Central  Google Scholar 

  17. Kijas J, Mcwilliam S, Naval SM, Kube P, King H, Evans B, et al. Evolution of sex determination loci in Atlantic salmon. Sci Rep. 2018;8(1):5664.

    Article  PubMed  PubMed Central  Google Scholar 

  18. Lin A, Xiao S, Xu S, Ye K, Lin X, Sun S, et al. Identification of a male-specific DNA marker in the large yellow croaker (Larimichthys crocea). Aquaculture. 2017;480:116–22.

    Article  CAS  Google Scholar 

  19. Ou M, Yang C, Luo Q, Huang R, Zhang AD, Liao LJ, et al. An NGS-based approach for the identification of sex-specific markers in snakehead (Channa argus). Oncotarget. 2017;8(58):98733–44.

    Article  PubMed  PubMed Central  Google Scholar 

  20. Zhang A, Huang R, Chen L, Xiong L, He L, Li Y, et al. Computational identification of Y-linked markers and genes in the grass carp genome by using a pool-and-sequence method. Sci Rep. 2017;7(1):8213.

    Article  PubMed  PubMed Central  Google Scholar 

  21. Marie-Annick M, Oliver TC. Potential threat of the international aquarium fish trade to silver arawana Osteoglossum bicirrhosum in the Peruvian Amazon. Oryx Int J Conserv. 2006;40(2):152–60.

    Google Scholar 

  22. Koenig LA, Gallant JR. Sperm competition, sexual selection and the diverse reproductive biology of Osteoglossiformes. J Fish Biol. 2021;99(3):740–54.

    Article  PubMed  Google Scholar 

  23. Bian C, Hu Y, Ravi V, Kuznetsova IS, Shen X, Mu X, et al. The Asian arowana (Scleropages formosus) genome provides new insights into the evolution of an early lineage of teleosts. Sci Rep. 2016;6: 24501.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Du K, Wuertz S, Adolfi M, Kneitz S, Stock M, Oliveira M, et al. The genome of the arapaima (Arapaima gigas) provides insights into gigantism, fast growth and chromosomal sex determination system. Sci Rep. 2019;9(1):5293.

    Article  PubMed  PubMed Central  Google Scholar 

  25. Cioffi MB, Rab P, Ezaz T, Bertollo L, Lavoue S, Oliveira EA, et al. Deciphering the Evolutionary History of Arowana Fishes (Teleostei, Osteoglossiformes, Osteoglossidae): Insight from Comparative Cytogenomics. Int J Mol Sci. 2019;20(17):4296.

    Article  PubMed  Google Scholar 

  26. Suzuki A, Taki Y, Urushido T. Karyotypes of Two Species of Arowana, Osteoglossum bicirrhosum and O.ferreirai. Japanese J Ichthyol. 1982;29(2):220–2.

    Google Scholar 

  27. de Oliveira EA, Bertollo L, Rab P, Ezaz T, Yano CF, Hatanaka T, et al. Cytogenetics, genomics and biodiversity of the South American and African Arapaimidae fish family (Teleostei, Osteoglossiformes). PLoS ONE. 2019;14(3):e0214225.

    Article  PubMed  PubMed Central  Google Scholar 

  28. Betancur-R R, Wiley EO, Arratia G, Acero A, Bailly N, Miya M, et al. Phylogenetic classification of bony fishes. Bmc Evol Biol. 2017;17(1):162.

    Article  PubMed  PubMed Central  Google Scholar 

  29. Imakaev M, Fudenberg G, Mccord RP, Naumova N, Goloborodko A, Lajoie BR, et al. Iterative correction of Hi-C data reveals hallmarks of chromosome organization. Nat Methods. 2012;9(10):999–1003.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Nagahama Y, Chakraborty T, Paul-Prasanth B, Ohta K, Nakamura M. Sex determination, gonadal sex differentiation, and plasticity in vertebrate species. Physiol Rev. 2021;101(3):1237–308.

    Article  CAS  PubMed  Google Scholar 

  31. Lemay DG, Lynn DJ, Martin WF, Neville MC, Casey TM, Rincon G, et al. The bovine lactation genome: insights into the evolution of mammalian milk. Genome Biol. 2009;10(4): R43.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Schrider DR, Hahn MW. Gene copy-number polymorphism in nature. Proc Biol Sci. 2010;277(1698):3213–21.

    CAS  PubMed  PubMed Central  Google Scholar 

  33. Brynildsrud O, Gulla S, Feil EJ, Norstebo SF, Rhodes LD. Identifying copy number variation of the dominant virulence factors msa and p22 within genomes of the fish pathogen Renibacterium salmoninarum. Microb Genom. 2016;2(4): e000055.

    PubMed  PubMed Central  Google Scholar 

  34. Freyaldenhoven BS, Freyaldenhoven MP, Iacovoni JS, Vogt PK. Avian winged helix proteins CWH-1, CWH-2 and CWH-3 repress transcription from Qin binding sites. Oncogene. 1997;15(4):483–8.

    Article  CAS  PubMed  Google Scholar 

  35. Sullivan SA, Akers L, Moody SA. foxD5a, a Xenopus winged helix gene, maintains an immature neural ectoderm via transcriptional repression that is dependent on the C-terminal domain. Dev Biol. 2001;232(2):439–57.

    Article  CAS  PubMed  Google Scholar 

  36. Herpin A, Schartl M. Plasticity of gene-regulatory networks controlling sex determination: of masters, slaves, usual suspects, newcomers, and usurpators. Embo Rep. 2015;16(10):1260–74.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Zhang L, Fang Y, Shi M, Ren K, Guan X, Younas W, et al. Gonadal expression profiles reveal the underlying mechanisms of temperature effects on sex determination in the large-scale loach (Paramisgurnus dabryanus). Anim Reprod Sci. 2024;272: 107661.

    Article  PubMed  Google Scholar 

  38. Wu K, Li Y, Pan P, Li Z, Yu Y, Huang J, et al. Gestational vinclozolin exposure suppresses fetal testis development in rats. Ecotoxicol Environ Saf. 2020;203: 111053.

    Article  CAS  PubMed  Google Scholar 

  39. Martinez-Juarez A, Lopez-Luna MA, Porras-Gomez TJ, Moreno-Mendoza N. Expression of the Sox9, Foxl2, Vasa, and TRPV4 genes in the ovaries and testes of the Morelet’s crocodile, Crocodylus moreletii. J Exp Zool B Mol Dev Evol. 2018;330(3):148–64.

    Article  CAS  PubMed  Google Scholar 

  40. Elzaiat M, Todeschini AL, Caburet S, Veitia RA. The genetic make-up of ovarian development and function: the focus on the transcription factor FOXL2. Clin Genet. 2017;91(2):173–82.

    Article  CAS  PubMed  Google Scholar 

  41. Eggers S, Ohnesorg T, Sinclair A. Genetic regulation of mammalian gonad development. Nat Rev Endocrinol. 2014;10(11):673–83.

    Article  CAS  PubMed  Google Scholar 

  42. Migale R, Neumann M, Mitter R, Rafiee MR, Wood S, Olsen J, et al. FOXL2 interaction with different binding partners regulates the dynamics of ovarian development. Sci Adv. 2024;10(12):eadl0788.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Natri HM, Shikano T, Merila J. Progressive recombination suppression and differentiation in recently evolved neo-sex chromosomes. Mol Biol Evol. 2013;30(5):1131–44.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Graves JA. Evolution of vertebrate sex chromosomes and dosage compensation. Nat Rev Genet. 2016;17(1):33–46.

    Article  CAS  PubMed  Google Scholar 

  45. Wright AE, Dean R, Zimmer F, Mank JE. How to make a sex chromosome. Nat Commun. 2016;7: 12087.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Dixon JR, Selvaraj S, Yue F, Kim A, Li Y, Shen Y, et al. Topological domains in mammalian genomes identified by analysis of chromatin interactions. Nature. 2012;485(7398):376–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Mohana G, Dorier J, Li X, Mouginot M, Smith RC, Malek H, et al. Chromosome-level organization of the regulatory genome in the Drosophila nervous system. Cell. 2023;186(18):3826-3844.e26.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Gonen N, Futtner CR, Wood S, Garcia-Moreno SA, Salamone IM, Samson SC, et al. Sex reversal following deletion of a single distal enhancer of Sox9. Science. 2018;360(6396):1469–73.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Croft B, Ohnesorg T, Hewitt J, Bowles J, Quinn A, Tan J, et al. Human sex reversal is caused by duplication or deletion of core enhancers upstream of SOX9. Nat Commun. 2018;9(1):5319.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Cocquet J, Pailhoux E, Jaubert F, Servel N, Xia X, Pannetier W, et al. Evolution and expression of FOXL2. J Med Genet. 2002;39(12):916–21.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Schmidt D, Ovitt CE, Anlag K, Fehsenfeld S, Gredsted L, Treier AC, et al. The murine winged-helix transcription factor Foxl2 is required for granulosa cell differentiation and ovary maintenance. Development. 2004;131(4):933–42.

    Article  CAS  PubMed  Google Scholar 

  52. Uhlenhaut NH, Jakob S, Anlag K, Eisenberger T, Sekido R, Kress J, et al. Somatic sex reprogramming of adult ovaries to testes by FOXL2 ablation. Cell. 2009;139(6):1130–42.

    Article  CAS  PubMed  Google Scholar 

  53. Yang YJ, Wang Y, Li Z, Zhou L, Gui JF. Sequential, Divergent, and Cooperative Requirements of Foxl2a and Foxl2b in Ovary Development and Maintenance of Zebrafish. Genetics. 2017;205(4):1551–72.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Li MH, Yang HH, Li MR, Sun YL, Jiang XL, Xie QP, et al. Antagonistic roles of Dmrt1 and Foxl2 in sex differentiation via estrogen production in tilapia as demonstrated by TALENs. Endocrinology. 2013;154(12):4814–25.

    Article  CAS  PubMed  Google Scholar 

  55. Zhang X, Li M, Ma H, Liu X, Shi H, Li M, et al. Mutation of foxl2 or cyp19a1a Results in Female to Male Sex Reversal in XX Nile Tilapia. Endocrinology. 2017;158(8):2634–47.

    CAS  PubMed  Google Scholar 

  56. Poletto AB, Ferreira IA, Cabral-De-Mello DC, Nakajima RT, Mazzuchelli J, Ribeiro HB, et al. Chromosome differentiation patterns during cichlid fish evolution. Bmc Genet. 2010;11:50.

    Article  PubMed  PubMed Central  Google Scholar 

  57. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Marcais G, Kingsford C. A fast, lock-free approach for efficient parallel counting of occurrences of k-mers. Bioinformatics. 2011;27(6):764–70.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Shao C, Li C, Wang N, Qin Y, Xu W, Liu Q, et al. Chromosome-level genome assembly of the spotted sea bass, Lateolabrax maculatus. Gigascience. 2018;7(11):1–7.

    Article  CAS  Google Scholar 

  60. Koren S, Walenz BP, Berlin K, Miller JR, Bergman NH, Phillippy AM. Canu: scalable and accurate long-read assembly via adaptive k-mer weighting and repeat separation. Genome Res. 2017;27(5):722–36.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Durand NC, Shamim MS, Machol I, Rao SS, Huntley MH, Lander ES, et al. Juicer Provides a One-Click System for Analyzing Loop-Resolution Hi-C Experiments. Cell Syst. 2016;3(1):95–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Dudchenko O, Batra SS, Omer AD, Nyquist SK, Hoeger M, Durand NC, et al. De novo assembly of the Aedes aegypti genome using Hi-C yields chromosome-length scaffolds. Science. 2017;356(6333):92–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Servant N, Varoquaux N, Lajoie BR, Viara E, Chen CJ, Vert JP, et al. HiC-Pro: an optimized and flexible pipeline for Hi-C data processing. Genome Biol. 2015;16:259.

    Article  PubMed  PubMed Central  Google Scholar 

  64. Manni M, Berkeley MR, Seppey M, Simao FA, Zdobnov EM. BUSCO Update: Novel and Streamlined Workflows along with Broader and Deeper Phylogenetic Coverage for Scoring of Eukaryotic, Prokaryotic, and Viral Genomes. Mol Biol Evol. 2021;38(10):4647–54.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Majoros WH, Pertea M, Salzberg SL. TigrScan and GlimmerHMM: two open source ab initio eukaryotic gene-finders. Bioinformatics. 2004;20(16):2878–9.

    Article  CAS  PubMed  Google Scholar 

  66. Holt C, Yandell M. MAKER2: an annotation pipeline and genome-database management tool for second-generation genome projects. BMC Bioinformatics. 2011;12: 491.

    Article  PubMed  PubMed Central  Google Scholar 

  67. Flynn JM, Hubley R, Goubert C, Rosen J, Clark AG, Feschotte C, et al. RepeatModeler2 for automated genomic discovery of transposable element families. Proc Natl Acad Sci U S A. 2020;117(17):9451–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. arailo-Graovac M, Chen N. Using RepeatMasker to identify repetitive elements in genomic sequences. Curr Protoc Bioinformatics. 2009;Chapter 4:Unit 4.10.

    Google Scholar 

  69. Katoh K, Standley DM. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol Biol Evol. 2013;30(4):772–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. Minh BQ, Schmidt HA, Chernomor O, Schrempf D, Woodhams MD, von Haeseler A, et al. IQ-TREE 2: New Models and Efficient Methods for Phylogenetic Inference in the Genomic Era. Mol Biol Evol. 2020;37(5):1530–4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  71. Marcais G, Delcher AL, Phillippy AM, Coston R, Salzberg SL, Zimin A. MUMmer4: A fast and versatile genome alignment system. Plos Comput Biol. 2018;14(1): e1005944.

    Article  PubMed  PubMed Central  Google Scholar 

  72. Li H. Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics. 2018;34(18):3094–100.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Thorvaldsdottir H, Robinson JT, Mesirov JP. Integrative Genomics Viewer (IGV): high-performance genomics data visualization and exploration. Brief Bioinform. 2013;14(2):178–92.

    Article  CAS  PubMed  Google Scholar 

  74. Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009;25(14):1754–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Mckenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, et al. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 2010;20(9):1297–303.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  76. Bradbury PJ, Zhang Z, Kroon DE, Casstevens TM, Ramdoss Y, Buckler ES. TASSEL: software for association mapping of complex traits in diverse samples. Bioinformatics. 2007;23(19):2633–5.

    Article  CAS  PubMed  Google Scholar 

  77. Brown-Peterson NJ, Wyanski DM, Saborido-Rey F, Macewicz BJ, Lowerre-Barbieri SK. A Standardized Terminology for Describing Reproductive Development in Fishes. Mar Coast Fish. 2011;3(1):52–70.

    Article  Google Scholar 

  78. Duponchelle F, Ruiz Arce A, Waty A, Garcia-Vasquez A, Renno JF, Chu-Koo F, et al. Variations in reproductive strategy of the silver Arowana, Osteoglossum bicirrhosum Cuvier, 1829 from four sub-basins of the Peruvian Amazon. J Appl Ichthyol. 2015;31:19–30.

    Article  Google Scholar 

  79. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  80. Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  81. Robinson MD, Mccarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.

    Article  CAS  PubMed  Google Scholar 

  82. Szklarczyk D, Franceschini A, Wyder S, Forslund K, Heller D, Huerta-Cepas J, et al. STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 2015;43(Database issue):447–52.

    Article  Google Scholar 

  83. Pfaffl MW. A new mathematical model for relative quantification in real-time RT-PCR. Nucleic Acids Res. 2001;29(9): e45.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  84. Osteoglossum bicirrhosum strain: prfri_Ob_female_2019 | breed:silver-arowana-female (silver arawana). Genebank. https://www.ncbi.nlm.nih.gov/bioproject/PRJNA881487 (2022).

  85. Osteoglossum bicirrhosum strain: prfri_Ob_male_2019 (silver arawana). Genebank. https://www.ncbi.nlm.nih.gov/bioproject/PRJNA881492 (2022).

Download references

Acknowledgements

The authors express their gratitude to Professor Jianren Luo, Professor Yinchang Hu, and Professor Xuejie Wang for their invaluable technical supports. The authors also thank all the students of the laboratory for their assistance.

Funding

This study was supported by Central Public-interest Scientific Institution Basal Research Fund, CAFS (nos. 2022SJ-XT2 and 2019ZD0503), Guangdong Provincial Special Fund for Modern Agriculture Industry Technology Innovation Team (no. 2022KJ150), China-ASEAN Maritime Cooperation Fund (no. CAMC-2018F), National Freshwater Genetic Resource Center (no. FGRC18537), and Guangdong Rural Revitalization Strategy Special Provincial Organization and Implementation Project Funds (2022-SBH-00–001).

Central Public-interest Scientific Institution Basal Research Fund,Chinese Academy of Fishery Sciences,2022SJ-XT2,2019ZD0503,Guangdong Provincial Special Fund for Modern Agriculture Industry Technology Innovation Team,2022KJ150,National Freshwater Genetic Resource Center,FGRC18537,Guangdong Rural Revitalization Strategy Special Provincial Organization and Implementation Project Funds,2022-SBH-00-001

Author information

Author notes

  1. Yi Liu, Chao Bian and Ka Yan Ma contributed equally to this work.

Authors and Affiliations

  1. Key Laboratory of Prevention and Control for Aquatic Invasive Alien Species, Ministry of Agriculture and Rural Affairs, Guangdong Modern Recreational Fisheries Engineering Technology Center, Pearl River Fisheries Research Institute, Chinese Academy of Fishery Sciences, Guangzhou, China

    Yi Liu, Yexin Yang, Yuanyuan Wang, Chao Liu, Guochang Ouyang, Meng Xu & Xidong Mu

  2. Laboratory of Aquatic Genomics, College of Life Sciences and Oceanography, Shenzhen University, Shenzhen, China

    Chao Bian & Qiong Shi

  3. State Key Laboratory of Biocontrol, School of Ecology, Sun Yat-Sen University, Guangzhou, China

    Ka Yan Ma

  4. Shanghai Ocean University, Shanghai, China

    Guochang Ouyang

  5. College of Fisheries, Tianjin Agricultural University, Tianjin, 300384, China

    Jinhui Sun

  6. Key Laboratory of Sustainable Development of Marine Fisheries, Ministry of Agriculture, Yellow Sea Fisheries Research Institute, Chinese Academy of Fishery Sciences, Laboratory for Marine Fisheries Science and Food Production Processes, Qingdao National Laboratory for Marine Science and Technology, Qingdao, China

    Changwei Shao

  7. Science Corporation of Gene (SCGene), Guangzhou, China

    Jiehu Chen

  8. Guangdong Provincial Key Laboratory of Aquatic Animal Immunology and Sustainable Aquaculture, Guangzhou, China

    Yexin Yang

Authors
  1. Yi Liu
    View author publications

    You can also search for this author inPubMed Google Scholar

  2. Chao Bian
    View author publications

    You can also search for this author inPubMed Google Scholar

  3. Ka Yan Ma
    View author publications

    You can also search for this author inPubMed Google Scholar

  4. Yexin Yang
    View author publications

    You can also search for this author inPubMed Google Scholar

  5. Yuanyuan Wang
    View author publications

    You can also search for this author inPubMed Google Scholar

  6. Chao Liu
    View author publications

    You can also search for this author inPubMed Google Scholar

  7. Guochang Ouyang
    View author publications

    You can also search for this author inPubMed Google Scholar

  8. Meng Xu
    View author publications

    You can also search for this author inPubMed Google Scholar

  9. Jinhui Sun
    View author publications

    You can also search for this author inPubMed Google Scholar

  10. Changwei Shao
    View author publications

    You can also search for this author inPubMed Google Scholar

  11. Jiehu Chen
    View author publications

    You can also search for this author inPubMed Google Scholar

  12. Qiong Shi
    View author publications

    You can also search for this author inPubMed Google Scholar

  13. Xidong Mu
    View author publications

    You can also search for this author inPubMed Google Scholar

Contributions

X.M., Y.L., C.B., and Q.S. designed research. Y.L., C.L., Y.W., X.W., and X.M. provided and collected samples. Y.L., Y.Y., YW and G.O. performed research. C.B, K.Y.M., J.S., C.S., M.X, R.L., and J.C. analyzed data. Y.L., K.Y.M. and M.X. wrote the paper with input from all authors. C.B., Q.S. and X.M. revised article. All authors contributed to discussions of the results and manuscript revision. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Qiong Shi or Xidong Mu.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

12915_2025_2139_MOESM1_ESM.docx

Additional file 1: Tables S1-S2. Table S1. Summary of each chromosome in both female and male genome assemblies of silver arowana. Table S2. Summary of the female and male genomes completeness assessment with BUSCO.

Additional file 2: Table S3. Genes for construction of the phylogenetic tree.

12915_2025_2139_MOESM3_ESM.docx

Additional file 3: Figures S1-S3. Fig.S1 Chromosomal alignments between male and female chromosomes. Fig. S2. The relative expression of foxl2 in the ovary of silver arowana after foxl2-dsRNA injection. Fig. S3. The relative expression levels of dmrt1, sox8a, sox9, and esr1 were evaluated 24 hours after the injection of foxl2-dsRNA.

Additional file 4: Table S4. Copy number of sex-related genes in various fishes.

Additional file 5: Table S5. Expression profiles of sex-related genes at three gonad developmental stages.

Additional file 6: Table S6. Correlations and R values for gene expression patterns.

Additional file 7: Table S7. Number of the samples used in this study.

Additional file 8: Table S8. Representative fish species used for genome evolution.

12915_2025_2139_MOESM9_ESM.docx

Additional file 9: Tables S9-S10. Table S9. Classification of the three gonadal developmental stages. Table S10. The specific primers employed in this study.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Liu, Y., Bian, C., Ma, K.Y. et al. Reference genome provide insights into sex determination of silver aworana (Osteoglossum bicirrhosum). BMC Biol 23, 29 (2025). https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12915-025-02139-5

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12915-025-02139-5

Keywords

BMC Biology

ISSN: 1741-7007

Contact us