Genomic deciphering of sex determination and unique immune system of a potential model species rare minnow (Gobiocypris rarus).

Gobiocypris rarus is sensitive to environmental pollution, especially to heavy metal and grass carp reovirus (GCRV). Hence, it has potential utility as a biological monitor. Genetic deciphering of its unique immune system will advance our understanding of its unique adaptive strategies, which provide cues for its better application. A de novo genome of rare minnow was obtained, and its sex determination mechanism is ZZ/ZW. We identified several specific mutation genes and specific lost genes of rare minnow, and these might be related to the sensitivity of rare minnow to environmental stimuli. We also analyzed the gene expression level of different organs/tissues and found that several IFIT genes may play key roles in GCRV resistance. In addition, knockout of the gene PCDH10L indicates that PCDH10L affects Pb2+-induced mortality in rare minnow. Rare minnow is ready for genetic manipulation and shows potential as an emerging experimental model.

INTRODUCTION

Rare minnow (Gobiocypris rarus) is a freshwater cyprinid fish that is small in size (adults are 38 to 85 mm), is easy to maintain, spawns frequently (every 4 days), has high fecundity, has short generation times (3 to 4 months), tolerates a range of temperatures (0° to 35°C), and is sensitive to environmental pollutants (). All these characteristics make it a suitable experimental model. It belongs to Cypriniformes, and the genetic background is close with that of Danio rerio; the closest order to Cypriniformes is Siluriformes (Ictalurus punctatus). Rare minnow is currently the only native Chinese fish species recommended by the national standard of China (GB/T29763-2013) in chemical testing methods and biological tests of wastewater detection (, ), and many studies have successfully used this species to evaluate the toxicological effects of both inorganic compounds () and organic pollutants (). Unlike the commonly used model fish such as D. rerio, the rare minnow is sensitive to the grass carp reovirus (GCRV), a double-stranded RNA virus () that causes hemorrhagic disease in grass carp (Ctenopharyngodon idellus) and results in major economic losses, such as 18.36% of the harvest of all freshwater Chinese fisheries in 2019 (). It is also an ideal model for investigating the molecular mechanisms of anti-GCRV infection immunity.

Sex determination mechanisms in fishes is a complex issue and has been triggering a lot of interest. The determination types include male-heterogametic gonochorism (XY), female-heterogametic gonochorism (ZW), hermaphroditism, and environmental sex determination. However, the sex determination mechanism in Gobiocypris remains unknown. Given the potential importance of rare minnow as a biological monitor, the sex determination mechanism is the important genetic foundation of rare minnow.

In this study, we performed the first genome sequencing, assembly, and annotation for rare minnow, in addition to a series of molecular evolutionary genomic analyses at the chromosome-level genome, to reveal the mechanism underlying the high sensitivity of these fish to environmental pollution, chemicals, and GCRV and identify genes related to environmental sensitivity through genome-wide analysis.

RESULTS AND DISCUSSION

Chromosome-level genome assembly and evaluation

To characterize the basic features of rare minnow genome, we generated 150.72 giga–base pairs (Gbp) of Illumina short–insert size reads with ~132.90× depth (table S1). The 17-mer analysis indicated that the genome size of rare minnow is ~1.04 Gb (fig. S1). Single-molecule sequencing (Oxford Nanopore Technologies), Illumina paired-end sequencing, and Hi-C sequencing were used (Illumina; tables S1 to S3) to further assemble the genome. Last, we obtained a 1.13-Gb genome with N50 value (the shortest sequence length at 50% of the genome size) of 42.38 Mb; a total of 26 chromosomes were anchored, and the mounting rate was 93.11% (figs. S2 and S3 and tables S4 and S5). The length of the chromosomes ranged from 31.70 to 60.35 Mb, comparable to that of the model organism (D. rerio; 37.50 to 78.09 Mb). The completeness, contiguity, connectivity, and accuracy of the reference genome were evaluated using different methods, such as BUSCO (97.7% in eukaryote and 97.3% in metazoa; tables S6 and S7), Illumina reads mapping ratio (90.05%; fig. S4 and table S8), assembled transcripts mapping ratio (96.12 to 98.93%; tables S9 to S11), and genome synteny between rare minnow and the two closest species (D. rerio and I. punctatus; Fig. 1A and fig. S5). We were able to find that all the 25 chromosomes in D. rerio have homologous chromosomes in rare minnow with strong genomic collinearity of the whole chromosomes, which means that the relationship and the genetic background between them are very similar. However, the Z chromosome (chromosome 26) in rare minnow has genomic collinearity with several chromosomes of D. rerio, including chromosomes 4, 3, 21, 8, 2, 15, 16, 12, and 9. All of these analyses indicated that the rare minnow genome had a high level of accuracy, continuity, and connectivity (fig. S6).

Fig. 1.

Genome assembly and comparative analysis.

(A) Genomic synteny between G. rarus, D. rerio, and I. punctatus. (B) Sequence collinearity between chr26 and scaf28. The blue module on the sequence represents the location and length of each syntenic block between the two sequences. The line shown between the two sequences indicates their syntenic relationships. The red module on the sequence represents Sox9 gene on Z chromosome. (C) Transposable element (TE) content of each chromosomes and each chromosome TE ratio in the genome. (D) Gene Ka/Ks value of each chromosome and gene length ratio of each chromosome in the genome. (E) Comparison of genome size and its composition among species. The x axis indicates the length of the nucleotide sequence of each genomic element, including the coding region, DNA element, long interspersed nuclear element (LINE), long terminal repeat (LTR), short interspersed nuclear element (SINE), and unknown TEs, among these species. (F) Comparison of the insertion history of TEs among these species. The x axis indicates the inferred insertion time (unit: million years ago) of TEs in the genome. The y axis represents the total/each length of the TE in each species.

Basic features of rare minnow genome

The average genome GC (guanine-cytosine) content of rare minnow (39.20%) is similar to that of D. rerio (36.60%) and slightly lower than that of other related species (Cynoglossus semilaevis: 40.81%, Hippocampus comes: 43.30%, I. punctatus: 39.70%, Labrus bergylta: 40.89%, Lepisosteus oculatus: 39.59%, Mola mola: 41.13%, Oreochromis niloticus: 40.73%, Oryzias latipes: 40.84%, Paramormyrops kingsleyae: 43.92%, Syngnathus scovelli: 42.95%, and Takifugu rubripes: 45.46%; fig. S7A). The average CDS GC content of rare minnow (50.67%) is also similar to that of D. rerio (49.46%) and slightly lower than that of the following species: C. semilaevis: 52.35%, H. comes: 55.21%, I. punctatus: 51.27%, L. bergylta: 53.10%, L. oculatus: 53.64%, M. mola: 53.34%, O. niloticus: 50.86%, O. latipes: 53.08%, P. kingsleyae: 55.51%, S. scovelli: 55.13%, and T. rubripes: 54.16% (fig. S7B). The codon usage in cyprinid fishes (rare minnow and D. rerio) differed from that of other species based on the codon adaption index (CAI), frequency of optimal codon (FOP), and frequency of guanine + cytosine at the synonymous third position of the codons (GC3s). The GC3 content of cyprinid fishes was slightly lower compared with that of other species, suggesting that the synonymous third position of the codons was always A or T. The distribution of the CAI (fig. S8A) and the FOP (fig. S8B) support this hypothesis. Cypriniformes (rare minnow and D. rerio) exhibits specific codon usage, which is more similar to Chloroforms (I. punctatus) and Coliforms (O. niloticus) (table S12).

Sex determination in rare minnow

Compared with most mammals and birds, sex determination mechanism in fish shows a high degree of plasticity and complexity, which may be related to genetic or environmental factors or both (). We attempted to determine the sex determination mechanism of rare minnow using both the assembled genome and sequenced reads. We found that the sequencing depth of chromosome 26 and scaffold 28 is nearly half of that in other chromosomes/scaffolds, suggesting that chromosome 26 and scaffold 28 might be sex chromosomes (fig. S4). Because the sequenced individual of rare minnow was female, the sex chromosome type should be ZZ/ZW, not XX/XY (fig. S9). We further assessed the candidate sex-determining gene in the genome and found Sox9 gene in chromosome 26 (chromosome Z; Fig. 1B). Sox9 gene is an important transcription factor expressed in embryonic and postnatal Sertoli cells of the testis, and its function is essential for the maintenance of testicular function (, ). Results of the quantitative real-time polymerase chain reaction (qRT-PCR) indicated that Sox9 gene is widely expressed in many organs/tissues (brain, liver, spleen, intestine, kidney, muscle, testis, and ovary), with specific high expression in testis and low expression in ovary (fig. S10). These results are consistent with the previous report (). Current theory of the evolution of sex chromosomes indicates that sex chromosomes are highly repetitive (, ). In particular, the W chromosome is considered to be degraded, containing much less genes, which causes it to be shorter than the Z chromosome, and more repetitive sequences. It is more difficult to completely assemble than the W chromosome. The transposable element (TE) content in Z chromosome (71.80%) is much higher than autosomes (50.02 to 58.47%), which is consistent with the current theory that sex chromosomes are highly repetitive (, ). The TE content in partial W chromosome (56.61%) is similar to the autosomes (50.02 to 58.47%); this may be caused by the incomplete assembly of W chromosome (Fig. 1C). Among all the types of TEs, the long terminal repeat (LTR) ratio (18.10%) and long interspersed nuclear element (LINE) ratio (10.69%) in Z chromosome are much higher than the LTR ratio (5.80 to 9.95%) and LINE ratio (3.61 to 6.58%) in autosomes (table S13); the higher TE content in Z chromosome is mostly caused by LTR and LINE types. The TE insertion time of all chromosomes showed that the TE insertion time of Z chromosome was much concentrated and happened much more recently, which suggested that the Z chromosome may have experienced rapid evolution in the ~recent 10 million years (Ma) (fig. S11). Furthermore, the gene density (coding sequence length/total chromosome coding sequences length) of each chromosome showed that the Z chromosome (1.77%) has a much lower value than other autosomes (2.86 to 5.28%; table S14), and the gene Ka/Ks value of Z chromosome (average Ka/Ks: 0.269688124) is significantly higher than those in other autosomes (average Ka/Ks: 0.153745172 to 0.205880801) by Wilcox test (P value, 2.2 × 10−16; Fig. 1D and table S15). All the chromosomes in rare minnow also show strong genomic synteny with those of D. rerio, and we can easily find chromosome homology between rare minnow and D. rerio, except for the Z chromosome (Fig. 1A). The Z chromosome in rare minnow has genomic collinearity with several chromosomes of D. rerio, including chromosomes 4, 3, 21, 8, 2, 15, 16, 12, and 9. All the evidence indicated that the sex chromosomes in rare minnow are specific for Gobioninae; this Gobioninae-specific sex chromosome should have emerged after the divergence between rare minnow and D. rerio and may have experienced rapid evolution.

Chromosome-level genome annotation

The genome size of cyprinid fishes (rare minnow and D. rerio) ranged from 1.13 to 1.68 Gb, much larger than that of the most closely related species I. punctatus (783.27 Mb). The DNA and LTR clearly had a burst in cyprinid fishes (Fig. 1E). The TE insertion times of D. rerio and rare minnow both occurred after their divergence time, which means that TE insertion events occurred independently in D. rerio and rare minnow (Fig. 1F and table S16). In the assembled genome, ~641.66 Mb of repeat sequences were identified, which accounts for 56.58% of the assembled rare minnow genome (table S17). DNA transposons (27.18%) and LTRs (7.86%) were the most enriched repeat elements, and short interspersed nuclear elements were the least prevalent (0.47%; table S18). The rare minnow genome was predicted to contain 25,146 protein-encoding genes, and the quality was comparable to published high-quality gene sets (fig. S12 and table S19). The gene function annotations revealed that most of the genes (25,337; ~95.09%) have homologs in common public databases, which indicates that the annotated genes were reliable (table S20). Although the genome size of D. rerio (1.68 Gb) is much bigger than rare minnow (1.13 Gb), the gene number (D. rerio: 26,533 and rare minnow: 25,146) and CDS length (D. rerio: 4.57 Mb and rare minnow: 4.15 Mb) between them are quite close, so the gene arrangement in rare minnow genome may be much compact than in D. rerio.

Comparative genomic analysis

Thirteen representative species (C. semilaevis, D. rerio, G. rarus, H. comes, I. punctatus, L. bergylta, L. oculatus, M. mola, O. niloticus, O. latipes, P. kingsleyae, S. scovelli, and T. rubripes) from Lepidosteiformes (the basal lobe-finned fishes) to other common ray-finned fishes group were selected for phylogenetic analysis. The phylogenetic trees showed that rare minnow was most closely related to D. rerio, followed by I. punctatus. Divergence time analysis indicated that G. rarus and D. rerio share a common ancestor that lived ~58.0 Ma ago, and Cypriniformes and Siluriformes share a common ancestor that lived ~162.5 Ma ago (Fig. 2A), which was verified by TimeTree website fossil records (table S21). As a result, 3043 single-copy genes were identified, and a total of 20,614 gene families were constructed (Fig. 2B). Rare minnow and D. rerio shared 13,161 common gene families, rare minnow and I. punctatus shared 12,329 common gene families (Fig. 2C), and D. rerio and I. punctatus were most closely related to rare minnow. Although Siluriformes (I. punctatus) and Cypriniformes (rare minnow and D. rerio) are different orders, they still share lots of common genes and have a similar genetic background. Siluriformes is the closest order to Cypriniformes.

Fig. 2.

Comparative genomics of rare minnow and other closely related species.

(A) Phylogenetic relationships among Actinopterygii. The red dot at the node represents a fossil record that was used for the calibration of the divergence time. The blue number in each node represents its divergence time for species. The red and green numbers in each node/species represent the expanded/contracted gene families, respectively. The coordinate axis below the phylogenetic tree shows the corresponding scale of divergence time. The color of the Latin name for each species corresponds to the taxon of species on the right. (B) Statistics of all orthologous/paralogous gene numbers in these species. (C) A Venn diagram displaying the overlap in orthologous genes in rare minnow and three other fishes. (D) Relative evolutionary rate of these fishes. The analysis was performed by the single-copy genes among these species with rare minnow as the reference species and L. oculatus as the outgroup species. The y axis indicates the relative evolutionary rate of species. The black dots indicate the specific magnitude of the relative evolutionary rate for each species. (E) Annotation and comparison of the Hox clusters among species. The box with a solid line indicates that its full length was directly annotated, and the box with the dotted line indicates that only its homeobox domain was present at this location. TSGD, teleost-specific whole-genome duplication.

There are 246 expanded gene families and a total of 3,977 genes identified in rare minnow, and they were enriched in biological processes involved in regulation of gene translation and protein metabolic process, autoimmune thyroid disease, and several signaling pathways, specifically Ca2+ signaling pathway (tables S22 and S23). The Ca2+ signaling pathway consists of a series of molecular biological events that link external stimuli to appropriate responses within cells. Ca2+ is a very important second messenger in cells. Many hormones and neurotransmitters achieve downstream signal transduction through the activation of Ca2+ signaling by intracellular Ca2+ elevation, further causing gene expression, substance secretion, motility, muscle contraction, and the release of neurotransmitters. This pathway has been proven to involve a heavy metal effect in vivo and in vitro (), and the expanded genes in Ca2+ signaling pathway may partially explain the sensitivity of heavy metal pollution of rare minnow. We checked and found that several genes in this pathway have copy number variation in rare minnow compared with D. rerio. Atp2a1, atp2a2b, atp2a3, cacna1a, cacna1c, cacna1e, cacna1g, cacna1h, cacna1i, cacna1s, egfra, erbb3, erbb4, prkca, mylka, nos2a, par1, plcb1, plcb4, ryr1a, ryr2a, and ryr3 genes have remarkably expanded in rare minnow, and these genes are all involved in the Ca2+ signaling pathway (fig. S13). These results of expanded genes in rare minnow strongly suggested us that the expanded Ca2+ signaling pathway genes cause the heavy metal sensitivity difference between rare minnow and D. rerio.

Evolutionary status of rare minnow

Species in the Cypriniformes (rare minnow and D. rerio) have a lower genome-wide evolutionary rate than Siluriformes (I. punctatus; Fig. 2D and tables S24 and S25), and the evolutionary rates of rare minnow and D. rerio is quite similar. Although the Hox genes in teleosts have been extensively studied (, ), as a potential model species, all background information of rare minnow should be integrated and presented. To further check the genetic diversity between rare minnow and other teleosts, especially for the genetic diversity between rare minnow and D. rerio, we identified the Hox genes in the 13 previously described species. Obviously, after the teleost-specific whole-genome duplication event that happened ~227.8 Ma ago, the number of Hox gene cluster is doubled. The Osteoglossiformes (P. kingsleyae) concurrently lost the HoxDb cluster and HoxCb cluster, but most of the other teleosts only lost the HoxCb cluster. The Cypriniformes (D. rerio and G. rarus) and Siluriformes (I. punctatus) lost the HoxDb cluster and retained the HoxCb cluster (, ). Compared with D. rerio, rare minnow lost HoxBb4, HoxCa11, and HoxCa12 and gained HoxAa2, HoxAb8, and HoxBb11. The Hox gene clusters in Cypriniformes are quite stable, and no obvious differences between Hox gene cluster in rare minnow and D. rerio were observed (Fig. 2E).

Specific gene loss in rare minnow

By genome synteny analysis, gene annotation, and sequenced read check, we found three genes—i.e., SNAP47, TIMP4, and TMEM244—that were widely present in the other 13 species (C. semilaevis, D. rerio, O. latipes, L. bergylta, O. niloticus, L. oculatus, P. kingsleyae, T. rubripes, M. mola, S. scovelli, H. comes, I. punctatus, and G. rarus) but lost in rare minnow. Synteny analysis of flanking regions containing the five upstream and downstream genes of these lost genes supports the loss of them (Fig. 3, A to C), and the five upstream and downstream genes of these lost genes were chosen for display. SNAP47 is involved in postsynaptic and presynaptic functions in the nervous system (). The heavy metals Pb, Cd, As, and Hg are widely distributed, and they all affect the synaptic function (–). Metal Pb disrupts the Ca2+ signaling in neuronal synapsis and causes the dysregulation of Ca2+-sensitive signaling pathways (). Metal Cd can produce memory loss and mental illness () that upset the balance between excitation and inhibition in synaptic neurotransmission (). Metal As exposure affects the synaptic activity of neurons localized in the hippocampus region of the brain in rat (). Metal Hg leads to the disruption of neurotransmitter and microtubules and the alteration of Ca2+ ion homeostasis (). All these common heavy metals can affect the synaptic function, and the loss of SNAP47 genes in rare minnow means that partial synaptic function may be lost in rare minnow and easily affected by the heavy metals. TIMP4 encodes a tissue inhibitor of metalloproteinase-4 and belongs to the TIMP gene family. It was reported to be involved in cancer () and plays important roles in the immune response of infection (, ). Matrix metalloproteinases facilitate leukocyte recruitment, cytokine and chemokine processing, and defensin activation (, ). Furthermore, there are series of three copies of TIMP4 gene in D. rerio, and all of these three copies were lost in rare minnow. The loss of all these copies of TIMP4 genes may degrade the immune function in rare minnow, which causes it to be susceptible and its immune system to be unique. TMEM244 is a transmembrane protein, and the TMEM gene family may be involved in physiological processes such as plasma membrane ion channels, activation of signal transduction pathways, autophagy (), and cancer (). Loss of TMEM244 gene in rare minnow may degrade the immune function in rare minnow and cause its immune system to be unique. Besides, all these three specific lost genes in rare minnow are related to Ca2+, directly or indirectly. These results indicated to us that these lost genes may have also affected the Ca2+ concentration and further affected the heavy metal metabolism in rare minnow, which is consistent with the Ca2+ signaling pathway genes expanded and further testifies that the Ca2+ signaling pathway may play important roles in the immune system function of rare minnow.

Fig. 3.

Lost genes and rare minnow–specific gene knockout.

(A) The absence of snap47 genes in rare minnow genome. Each exon is shown in this figure. (B) The absence of timp4 genes in rare minnow genome. (C) The absence of TMEM244 genes in rare minnow genome. (D) Schematic of the genomic structure of rare minnow PCDH10L. Sequence information of PCDH10L and PCDH10L. (E) Expression of PCDH10L mRNA in larvae from the WT (wild type), PCDH10L, and PCDH10L. **P < 0.01; ****P < 0.0001. (F) The dose-response curve of the death rate of PCDH10L and PCDH10L exposed to Pb2+ with concentrations ranging from 50 to 3200 μg/liter for 96 hours.

Positively selected genes and specific mutation genes in rare minnow

Positively selected genes of rare minnow using five cyprinid species (rare minnow, D. rerio, Cyprinus carpio, Carassius auratus, and Onychostoma macrolepis) were calculated, and the three genes PUSL1, RSRC1, and TMEM141 were identified. PUSL1 encodes an inner membrane–associated mitochondrial matrix protein required for efficient mitochondrial translation (). RSRC1 encodes a member of the serine- and arginine-rich–related protein family, associated with altered brain function in schizophrenia, and may also be associated with cancer cell proliferation and migration (, ). Genome-wide association studies have shown that single-nucleotide polymorphisms of RSRC1 are associated with schizophrenia (, ), and its mutation causes intellectual disability, aberrant behavior, hypotonia, and mild facial dysmorphism (). This gene mutation in rare minnow may also cause their behavior and brain function to be different from other cyprinid fishes. Unexpectedly, the TMEM141 gene and the previously described lost TMEM244 gene both belong to the TMEM gene family, which indicates to us that the TMEM gene family plays important roles in the immune function of rare minnow and causes its immune system to be unique. TMEM141 has six potential positively selected sites (fig. S14A). The three-dimensional (3D) structure of TMEM141 of the outgroup species and rare minnow differed at the ends of both sites, which demonstrates that the structures of TMEM141 in rare minnow and outgroup species are different (fig. S14B).

Furthermore, we identified rare minnow–specific mutation genes in these five cyprinid species and I. punctatus. There were 339 specific mutation genes in rare minnow; the enrichment analysis showed that these genes are related to the regulation of diverse biological processes, such as signaling regulation, regulation of signal transduction, and regulation of response to stimulus (tables S26 and S27). The signaling regulation item included several Ca2+ signal genes, which is associated with the expanded gene families in rare minnow. PCDH10L (Protocadherin 10) gene is the only typical calcium-dependent cell adhesion protein with rare minnow–specific mutations. It has been reported to be associated with autism spectrum disorder (), but its involvement in heavy metal toxicity has never been proven. Hence, we knocked out PCDH10L via CRISPR-Cas9 and generated two distinct mutant lines (PCDH10L: 5-bp deletion; PCDH10L: 17-bp deletion; Fig. 3D). qRT-PCR analyses revealed that PCDH10L mRNA was greatly reduced in these two mutants (Fig. 3E), indicating the nonsense-mediated decay of rare minnow PCDH10L mRNA. Knockout of PCDH10L did not cause obvious phenotypic abnormalities. Intriguingly, however, the mutants did show notably lower sensitivity to Pb2+. We used the semistatic bioassays to determine the median lethal concentration (LC50) of lead toxicant. We found that the LC50 values for exposure to Pb2+ of these two mutants were 507.98 μg/liter (PCDH10L) and 802.23 μg/liter (PCDH10L), respectively, much higher than that of the wild type (423.01 μg/liter; Fig. 3F). These results strongly suggest that PCDH10L plays an important role in the sensitivity of rare minnow to Pb2+. Further functional exploration of the key mutation in PCDH10L will help in understanding the high sensitivity of rare minnow to Pb2+.

GCRV resistance in D. rerio and rare minnow

Except for the heavy metals, rare minnow is also very sensitive to GCRV. Compared with the control group, significant hemorrhagic symptoms were observed in the GCRV-infected rare minnow, and there were no obvious hemorrhagic symptoms in both the GCRV-infected and control groups in D. rerio (fig. S15). The mortality curves of the rare minnow revealed that the first death was recorded as early as 5 days after infection in the GCRV-II–infected group, and 100% total mortality was reached at 6 days after infection; however, mortality remained at 0% in the control group. In the D. rerio group, although the first death was recorded at 5 days after infection in the GCRV-II–infected group, total mortality was 42.5% at 8 days after infection, much lower than that in the rare minnow (Fig. 4A). These results suggested that rare minnow is more sensitive to the GCRV than D. rerio.

Fig. 4.

GCRV infection in rare minnow and D. rerio.

(A) Cumulative mortality of rare minnow and D. rerio infected with GCRV. The number of dead fishes in the GCRV-infected group was counted daily. Phosphate-buffered saline was the negative control. (B) Heatmap of highly expressed genes after GCRV treatment in rare minnow. (C) Heatmap of highly expressed genes after GCRV treatment in D. rerio. (D) CIK cells were seeded in six-well plates overnight and transfected with GrIFIT8, 10, 14-Flag, or pCMV-Flag. After 24 hours posttransfection, the CIK cells were infected with GCRV (multiplicity of infection = 0.01) and incubated at 28°C. The total RNAs were extracted to examine the mRNA levels of VP2 and VP7 transcripts of GCRV by qRT-PCR. *P < 0.05; **P < 0.01; ***P < 0.001.

We chose five different immune-related tissues/organs (gill, head kidney, kidney, liver, and spleen) (, ) and compared the transcriptomes of these tissues/organs between D. rerio and rare minnow. The principal components analysis results showed that three replicates of each tissue in both species were clustered (fig. S16). We further detected 84 up-regulated genes of the GCRV treatment group of all tissues/organs in rare minnow (Fig. 4B and fig. S17A) and 107 up-regulated genes of the GCRV treatment group of all tissues/organs in D. rerio (Fig. 4C and fig. S17B). In addition, coexpressed modules were constructed using these highly expressed genes in D. rerio and rare minnow by WGCNA with a soft threshold of β = 14 and 8, respectively (fig. S18). In total, 11 modules and 13 modules in rare minnow and D. rerio were identified (figs. S19 and S20). To further understand the GCRV resistance between them, we performed correlation analyses between modules and tissues (fig. S21). The results indicated that different organs/tissues have different core modules, and different organs/tissues have different genes participating in GCRV resistance.

In these genes, the IFIT gene families show a different expression pattern between the control group and GCRV-infected group. IFIT8, IFIT10, and IFIT14 were all up-regulated in the GCRV-infected group of rare minnow, which indicates that the IFIT gene families may play important roles in GCRV resistance. The IFIT gene family is a group of antiviral RNA binding proteins that were the most highly expressed genes during antiviral immune responses; they also interact with other cellular proteins to regulate the host antiviral response by modulating innate immune signaling and apoptosis (). This gene family is conserved in all vertebrates and evolved together with the interferon system (). The expression of VP2 and VP7 RNA in CIK cells (grass carp kidney cells) transfected with GrIFIT8-Flag was significantly inhibited compared with that in CIK cells transfected with pCMV-Flag (Fig. 4D); similarly, the obviously repressed RNA level of VP2 and VP7 in CIK cells transfected with GrIFIT10-Flag or GrIFIT14-Flag was also observed, which suggests that GrIFIT8, GrIFIT10, and GrIFIT14 can inhibit GCRV replication. In addition, compared with the control group, the cytopathic effect of GrIFIT8, GrIFIT10, and GrIFIT14 was significantly reduced 24 hours after GCRV injection (fig. S22). In summary, our data indicate that IFIT8, IFIT10, and IFIT14 may play crucial roles in the host antiviral innate immune response.

Major histocompatibility complex gene mutation in rare minnow

The genomic organization of the major histocompatibility complex (MHC) plays various important roles in immune function; however, genome organization varies greatly among vertebrates, especially in fishes. The MHC cluster in D. rerio was present on three different chromosomes: 19, 16, and 15. The MHC cluster in rare minnow was also present on three chromosomes (10, 17, and 2), and the order of MHC genes is different from D. rerio. However, the MHC gene cluster in I. punctatus was present on six different chromosomes (14, 26, 12, 1, 23, and 17; Fig. 5). The human MHC is a cluster of more than 200 genes, broadly characterized as class I, II, and III genes, that is present on a single chromosome (chromosome 6) (, ). The ancestral organization of the ray-finned MHC has remained unresolved because of the low-quality assembly of the MHC cluster in spotted gar, but at least one pair of class I and class II genes is linked in spotted gar (). These findings indicated that the Cypriniformes (D. rerio and rare minnow) may retain the ancestral configuration. I. punctatus may differ from the ancestral state. In addition, lypc in rare minnow shows copy number variation. The lypc gene has only one copy in D. rerio, but in rare minnow, there are three copies, distributed in chromosomes 2 and 10. More gene copies mean that much more functions are enhanced, and this gene is related to circadian rhythm (), which indicates to us that the circadian rhythm of rare minnow may be different from that of D. rerio. Overall, rare minnow is an excellent emerging model organism for aquatic toxicity testing and chemical safety assessment and is particularly valuable for biological and evolutionary research and wastewater detection.

Fig. 5.

Distribution of the MHC genes in the three species.

Each black line represents one MHC gene. Orthologs are connected with a green/orange line.

METHODS

Materials

The adult female (3 months old) rare minnow used for library construction and sequencing was obtained from the Institute of Hydrobiology, Chinese Academy of Sciences (Wuhan, China). The adult (3 months old) Tübingen strain of D. rerio used for the experiment was obtained from the China Zebrafish Resource Center (Wuhan, China).

Library construction and sequencing

Whole genomic DNA of rare minnow was extracted from muscle using a DNeasy Blood & Tissue kit (Qiagen). Next, the quality of the extracted genomic DNA was examined using an Agilent 2100 Bioanalyzer (Agilent Technologies). RNA was extracted for all samples using TRIzol (Invitrogen) and subsequently purified using an RNeasy Mini kit (Qiagen). The libraries for long read sequencing were sequenced on the Nanopore platform. For short read sequencing, a paired-end library with an insert size of ~300 bp and a 100-bp paired-end read length was constructed and sequenced with the Illumina HiSeq 4000 platform with a 100-bp paired-end read length. For Hi-C sequencing, the libraries were sequenced on the Illumina NovaSeq 6000 platform with a 150-bp paired-end read length.

Quality control of raw reads

Two different strategies were used for short read and long read filtering. For Nanopore long read data, those with short length (read length of <1 kb) or low quality (read average quality of <7) were discarded by Perl scripts. For Illumina data, the raw sequencing reads were filtered with in-house Perl scripts according to the evaluation reports. Next, the adaptor sequences, any read with more than 10% unknown bases (represented by “N”), and reads with low-quality bases greater than 30% were discarded.

Estimation of genome size

The genome size of rare minnow was estimated using all of the filtered short reads with the K-mer–based method (). The K-mer number used was 17, and the genome size was estimated using the total number of 17-mers divided by the peak 17-mer frequency. Because 17-mers with low frequency mostly resulted from sequencing errors, we discarded all of the 17-mers with counts less than 2.

De novo genome assembly

De novo genome assembly using the filtered long reads was carried out using NextDenovo (v2.0-beta.1; https://github.com/Nextomics/NextDenovo) with default parameters as follows: read_cutoff = 1k; seed_cutoff = 16,363; blocksize = 4g; pa_correction = 15; seed_cutfiles = 30; sort_options = -m 30g -t 3 -k 50; minimap2_options_raw = -x ava-ont -t 3; and correction_options = -p 8. To further improve the assembled draft genome quality and reduce the ratio of single-base errors, BWA software (v0.7.17, BWA-MEM module) was used to map short reads to the genome with default parameters; Pilon software (v1.22) was then used for error correction. After error correction, all Hi-C reads were used for chromosome construction. The Hi-C reads were aligned to these contigs, and anchoring was performed with 3D DNA assembly software (v180419). The chromosomes were generated by connecting two adjacent sequences containing 500-bp N bases. The assembly integrity of the final assembled genome was determined using the final genome size divided by the 17-mer–estimated genome size.

Genome quality evaluation

To evaluate the complement and accuracy of the assembled genome, several strategies from the core gene complement, reads mapping ratio, and expressed transcript mapping ratio were used. First, the eukaryota core gene and metazoa core gene complement were both identified using BUSCO (v2.0). Second, all filtered short Illumina reads were mapped to the assembled genome using BWA software (v0.7.12) with default parameters to determine genome integrity and accuracy. The depths of positions were calculated using the “depth” command in Samtools software (v1.3.1). The rare minnow genome was bound into nonoverlapping windows (window size of 10 kb), and any window with N content of >0.1 was filtered. GC content and the average sequencing depth in each window were then calculated and plotted. Third, all of the RNA sequencing (RNA-seq) reads of five different tissues/organs (gill, head kidney, kidney, liver, and spleen) with three repeats for each tissue/organ were assembled using Bridger (r2014-12-01), and the assembled transcripts were mapped to the genome using blastn (v2.3.0) with the e value set as 10−5. Blast results with identity of more than 85% were selected, and the mapping ratio of all transcripts was calculated and plotted using R.

Tandem repeats and TE annotation

The rare minnow repeat sequences were predicted using RepeatMasker (v.4.0.7) and RepeatProteinMask (v1.36) with the Repbase TE library, and a de novo repeat library was constructed using RepeatModeler (v.1.0.10) with default parameters. After the TE sequences were obtained, the kimura value was extracted from the RepeatMasker results, and the evolutionary rate was calculated using the r8s model. The insertion time of TE sequences was calculated using the divergence rate divided by the neutral mutation rate (, ). As a complement, the tandem repeat sequences were predicted using Tandem Repeat Finder (v4.09) with default parameters. We also applied the same protocol to annotate the genome sequences of D. rerio and I. punctatus for further comparison.

Protein-coding gene annotation

Protein-coding genes were predicted by combining ab initio prediction, homology-based prediction, and evidence from mapped RNA-seq data, separately. First, for de novo prediction, the assembled transcripts described before, along with the rare minnow genome, were used for Augustus training to obtain an accurate Augustus annotation species model; Augustus software (v2.5.5) was then used for gene annotation with default parameters. Second, for homology-based annotation, proteins from D. rerio (GCF_000002035.6) (), Branchiostoma floridae (GCF_000003815.1) (), Callorhinchus milii (GCF_000165045.1) (), T. rubripes (GCF_000180615.1) (), Latimeria chalumnae (GCF_000225785.1) (), Larimichthys crocea (GCF_000972845.1) (), Rhincodon typus (GCF_001642345.1) (), I. punctatus (GCF_001660625.1) (), Paralichthys olivaceus (GCF_001970005.1) (), and Homo sapiens (GCF_000001405.38) () were downloaded from the National Center for Biotechnology Information (NCBI), and the longest transcript of each gene was selected. All of the protein sequences for the remaining genes were aligned using tblasn with an e value of 1 × 10−5. The blast result formats were changed and prepared for gene structure prediction using GeneWise. Third, the assembled transcripts were mapped to the genome using BLAT (v34, with identity of >90% and coverage of >90%), and the gene structure was inferred using PASA (v2.1.0). At last, these three strategies were combined using EvidenceModeler (v1.1.1).

Gene functional annotation

Several different public protein databases were used for functional annotation, including InterProScan, Kyoto Encyclopedia of Genes and Genomes, Swiss-Prot, TrEMBL, Cog, and NR. InterProScan (v4.8) was used to screen proteins against five databases (Pfam, release 27.0; prints, release 42.0; prosite, release 20.97; ProDom, 2006.1; and smart, release 6.2). The other databases were all used for annotation by BLAST software (v2.3.0) with the e value set as 10−5.

Genome synteny

Using the rare minnow genome as the reference genome, the genomes of D. rerio and I. punctatus were aligned to the rare minnow genome using LAST (“lastal” command in LAST with parameters set as: -P 5 -m100 -E 0.05; v802), and the alignment files in maf format were sorted (“maf-sort” command in LAST). The one-to-one aligned sequences were plotted using Circos (v0.69-6).

Phylogenetic inference

To clarify the phylogenetic relationships of rare minnow, 13 species (C. semilaevis, D. rerio, O. latipes, L. bergylta, O. niloticus, L. oculatus, P. kingsleyae, T. rubripes, M. mola, S. scovelli, H. comes, I. punctatus, and G. rarus) were used for reciprocal BLAST (v. 2.3.0) analysis, and OrthoMCL software (v2.0.9) was used to determine the homology relationships among these species. Self-to-self alignments were performed using BLASTP with an e value of 1 × 10−5, and low-quality hits (identity of <30% and coverage of <30%) were removed. Orthologous groups were constructed by OrthoMCL using the default settings based on the previous BLASTP results. After obtaining these orthologous genes, all of these genes were connected into one super-gene in each species. Next, MUSCLE (v3.8.31) was used to align these sequences, and the corresponding CDS alignments were back-translated from the corresponding protein alignments using PAL2NAL. For phylogenetic reconstruction, both protein and CDS alignments of each single-copy gene were input into RAxML (v8.2.9) to infer a maximum likelihood tree with 100 bootstrap replicates; L. oculatus was used as the outgroup species.

Molecular clock analysis

To estimate the divergence time among species, all of the 4d sites were extracted from the super-sequence by Perl scripts. The divergence times were inferred using a relaxed molecular clock with autocorrelated rates, which was implemented in MCMCTree within the PAML package (v4.9) with samples drawn every 2000 steps and a total of 100,000 samples. Fossil record data were derived from TIMEtree database (www.timetree.org).

Relative evolutionary rate

The evolutionary rates of these 13 species were calculated using the CDS sequences of all the single-copy genes. All genes were first connected into one super-sequence and aligned using MUSCLE (v3.8.31). We used L. oculatus as the outgroup species; LINTRE software and MEGA (Tajima’s relative rate test) software were used for this analysis. In Tajima’s test, the higher number of lineage-specific substitutions corresponds to a faster evolutionary rate. In the LINTRE method, the evolutionary rate of each species was assessed using Z statistics and the tpcv module with –distance=7.

Expansion and contraction of the gene family

The gene birth and death rates (lambda) and gene families showing significant dynamics were determined by CAFE (v3.1) with the following parameters: -p 0.01, -r 1000, -s. Three results—the phylogenetic relationships, divergence time, and fossil records—were used as inputs in this analysis.

Lost genes in rare minnow

All of the genes in these 13 species were annotated using the Swiss-Prot database. Genes that were present in all 12 species and not present in rare minnow were used for further analysis. The protein sequences of these lost candidate genes were blasted to the rare minnow genome using tblastn with the e value set as 1 × 10−3. Only genes with an identity less than 50% and alignment length less than 50 bp were considered lost genes. Then, the Illumina reads, Hi-C reads, and RNA-seq reads of all the tissues were used to identify lost genes using the same parameters and procedures. Last, the lost genes were checked using the genome synteny.

Positive selection analysis

Four cyprinid species (C. carpio, C. auratus, O. macrolepis, and D. rerio) and rare minnow were used to assess the positively selected genes. To estimate the lineage-specific evolutionary rate of each branch, the free model (“model = 1, NSsites = 0”) in the codeml program was used in Paml (v4.9). After a general evolutionary pattern of selective pressure along the lineages was obtained, a branch-site model was used. We used the likelihood ratio test and chi-square test to calculate the P value, and positively selected genes with a P value less than 0.05 were conserved.

Identification of sex chromosomes

After the chromosome-level genome was assembled using Hi-C data, the depths of the final chromosomes and scaffolds were assessed using Illumina reads, and the average depth of each chromosome/scaffold was calculated. Because sex chromosomes may have almost half the depth of other chromosomes, the Z chromosome and W chromosome should show significant genomic synteny. In addition, unlike the autosomes, the X chromosome typically features male-beneficial genes (“demasculinization”), and the Y chromosome features male-related genes (“masculinization”) (). We thus also characterized the core sex-determining genes on these two sex chromosomes. TE content, gene density, and Ka/Ks of all genes in all chromosomes were all checked as well.

Gene expression

Differentially expressed genes among species in the GCRV treatment were identified. The gene expression value [fragments per kilobase of exon per million mapped fragments (FPKM)] was calculated using the classic Tophat-Cufflinks pipeline. At first, the RNA-seq reads were mapped to the genome using Tophat; the transcript assembly was conducted using Cufflinks. All GTF files were then merged with the reference GTF by Cuffmerge, and the expression levels were calculated by Cuffdiff.

Specific mutation genes of rare minnow in cyprinid

The longest transcript of the genes of six species (rare minnow, D. rerio, C. carpio, C. auratus, O. macrolepis, and I. punctatus) was selected and annotated using the Swiss-Prot database. The gene sequences were then aligned using MUSCLE (v3.8.31), and rare minnow mutation sites were detected using Perl scripts. Gaps in the alignment files of rare minnow more than 30% were removed to ensure accuracy. Only rare minnow genes with specific mutation sites different from other species were used for further analysis.

Experimental sample maintenance, Pb2+ treatment solution, and LC50 assessment

For the lead toxicity analysis, the rare minnows were maintained in flow-through aquariums at 26° ± 1°C with a photoperiod of 12/12 hours (light/dark) and were fed to satiety twice daily with Chironomus larvae. The Pb2+ treatment solution and LC50 assessment were performed as reported previously (). Briefly, main stock toxicant of Pb2+ solution (i.e., 1000 parts per million) was prepared by dissolving 1.599 g of lead (II) nitrate (Merck, Germany) in 1000 ml of ultrapure deionized water from Millipore instruments (Hi-tech Instruments Co., Shanghai, China) and mixed thoroughly for 1000 mg/liter. The concentrations of Pb2+ solution ranged from 100 to 1000 μg/liter for 96 hours. Then, triplicates of 5 liters of experimental liquid volume were obtained with seven groups of nominal concentration gradients using geometric series as interval. There were seven fish per glass tank for the experiment, and a blank control was also set up. During acute lead toxicity tests, physicochemical parameters (including pH, water temperature, specific conductivity, dissolved oxygen, and oxygen saturation) were analyzed every 24 hours. Their mean length and weight (±SD) were 3.68 ± 0.52 cm and 0.67 ± 0.22 g, respectively. The death rate of the experimental fish, which were reared for 14 days and fasted for 24 hours before metal exposure, was less than 2%.

For the GCRV infection analysis, adult rare minnows weighing from 0.69 to 1.35 g (average length of 4.50 cm) and D. rerio weighing from 0.45 to 0.80 g (average length of 3.99 cm) were used. For the qRT-PCR experiments of Sox9, adult male rare minnows (n = 3) weighing from 0.56 to 0.80 g (average length of 3.57 cm) and female rare minnows (n = 3) weighing from 0.91 to 1.03 g (average length of 4.00 cm) were used. All of the experimental procedures were approved by the Animal Care and Use Committee of the Institute of Hydrobiology, Chinese Academy of Sciences (approval protocol no. Y21201-1-101), and all experiments were conducted as per the guidelines of the committee. The Pb2+ treatment and LC50 assessment were performed as described previously ().

CRISPR-Cas9 target prediction for guide RNA selection

The disruptions of PCDH10L were generated by CRISPR-Cas9 technology as described in a previous study (). The genome datasets of the rare minnow are provided as resources. For each selected gene, the conserved regions in the coding sequence were identified using multiple vertebrate orthologs via http://genome.ucsc.edu/. The guide RNAs (gRNAs) of each gene were designed using the tools provided in the Massachusetts Institute of Technology online CRISPR software (http://crispr.mit.edu). Because PCDH10L has two exons and one intron, gRNAs were designed to target a region in the first exon immediately downstream from the translation start sites in this gene. Two distinct mutant lines (PCDH10L: 5-bp deletion; PCDH10L: 17-bp deletion) were identified using the primers listed in table S28.

Viral challenge experiment and sample collection

Adult rare minnows and D. rerio (both 3 months old) were randomly divided into two groups (approximately 60 per group). Each fish in group I was infected with 50 μl of GCRV-HZ08 (2.97 × 103 RNA copies/μl) by intraperitoneal injection, and fish in group II were injected with 50 μl of phosphate-buffered saline as a control. At 120 hours after injection, samples of gill, spleen, liver, kidney, and head kidney tissues were harvested from both groups (n = 3 biological replicates). The remaining fish were monitored carefully, the cumulative mortality of the fishes infected with GCRV-II (GCRV-HZ08) was determined, and the number of dead fish in the GCRV-infected group and control group was counted daily. The total mortality was calculated when no mortality was recorded for seven consecutive days, and then the experiment was terminated.

Plasmid construction, transient transfection, and viral infection

The open reading frames of IFIT8, IFIT10, and IFIT14 were subcloned into the pCMV-Flag vector (Clontech). All recombinant plasmids were verified by DNA sequencing. Next, CIK cells were seeded in 24-well plates overnight and transfected with 2 μg of GrIFIT8-Flag or the pCMV-Flag vector using GeneTwin Transfection Reagent (Biomed) as per the manufacturer’s protocol. At 24 hours after transfection, the CIK cells were infected with GCRV at a multiplicity of infection = 0.01 and incubated at 28°C. The plaque analysis and the detection of virus titers were performed as described previously (). The primers used in this section are listed in table S28.

qRT-PCR analysis

The genome of GCRV consists of 11 segments (termed s1 to s11), and these 11 genome segments encode seven structural proteins (VP1 to VP7) and five nonstructural proteins (). We designed primers to amplify VP2 and VP7 using qRT-PCR, which was performed with FastSYBR Green PCR Master Mix (Bio-Rad) on the Applied Biosystems StepOne Real-Time PCR System. In addition, the total RNA of brain, liver, spleen, intestine, kidney, muscle, testis, and ovary was extracted from three randomly selected rare minnow specimens, and qRT-PCR was performed to investigate mRNA expression patterns of the Sox9 gene in vivo. All the primers used in this section are listed in table S28. PCR conditions were as follows: 95°C for 5 min, followed by 45 cycles of 95°C for 20 s, 60°C for 20 s, and 72°C for 20 s. There were three replicates for each sample, and the β-actin gene was used as an internal control to normalize the gene expression. The length of PCR products was approximately 100 to 200 bp. Relative expression levels were calculated using the 2−ΔΔCt method ().

Identification of MHC gene clusters

The MHC gene positions and sequences of D. rerio were downloaded from NCBI (GCF_000002035.6) and used as reference sequences; all of these sequences were used for gene annotation in rare minnow and I. punctatus. Blast and GeneWise were used for gene annotation, and annotated genes with more than 30% identity and 80 alignment ratios with MHC genes of D. rerio were used. The MHC figure of D. rerio was drawn using the file downloaded from NCBI; the MHC figure of rare minnow and I. punctatus was drawn using annotated genes.

Inference of demographic history

The standard pairwise sequentially Markovian coalescent (PSMC) pipeline () was used to infer the demographic history of rare minnow. BWA (v0.7.12) was first used to align the Illumina short reads to the genome, and the mpileup module in bcftools (version using htslib 1.9) was used to generate the bcf file for each sequence in the genome with the key parameters “-d 150 -q 20 -Q 20.” After merging all of the results, the fq2psmcfa (psmc-master) subprogram was used with the parameters “-q 20.” On the basis of these results, a final PSMC analysis with 100 bootstrap replicates was performed with the parameters “-N25 -r5 -p 4+25*2+4+6,” where -N represents the maximum number of iterations, and -p defines the atomic time intervals. The results were delivered to the plot script using values g = 0.5, which means the number of years per generation. We calculated u, which defines the absolute mutation rate per nucleotide, using r8s (v1.71) with default parameters.