The whole process of the experiment using the fish was carried out in accordance with the guidelines of the Institutional Animal Care and Use Committee (IACUC) of CellQua (CQ-IACUC-2023-1-2) and Gangneung-Wonju National University (GWNU-2019-26). All methodologies of the experiment were verified by CellQua and GWNU ethical committee. Also, the authors of the present study have acquired Animal Welfare & Ethics Course certification, research ethics and compliance training program.

Silvertip tetra (H. nana) used in the experiment was in the fully matured adult stage (F1), and their average weight and length were 1.0 ± 0.038 g, 3.54 ± 0.075 cm in female (n=5), and 0.9 ± 0.019 g, 3.58 ± 0.037 cm in male (n=5), respectively ( Figure 1 ). Fish (F0) were purchased at a local aquarium (Busan, Korea), and the F1 generation was produced by natural breeding F0 and used for experiments. Genomic DNA was extracted from the fin tissue of the silvertip tetra for species identification, and the Cytochrome c oxidase subunit I (COI) gene was analyzed ( Supplementary Table S1 ). Embryo development was observed and recorded each stage by microscope (Motic, Xiamen, China). Fish and fertilized eggs were kept in a small-scale recirculation aquaculture system with at 24 ± 1.0°C and supplied different extruded pellets or artemia according to the growth stage. For the measure of body-size and -weight, and tissue collection, fish were anesthetized with 200 ppm of clove oil (Thermo Fisher Scientific, MA, USA). Tissue samples (brain-pituitary, caudal fin, and gonad); from females and males were kept in an ultra-freezer at -80°C for total RNA extraction. Total RNA was isolated from each tissue using Ribo EX Reagent (GeneAll, Seoul, South Korea) and purified with Hybrid-RTM Kit (GeneAll) according to the instructions of the manufacturer. The quantity and quality of the total RNA were measured using the ND-1000 nanometer (Thermo Fisher Scientific) and Bioanalyzer (Agilent, Santa Clara, CA).

Primary cells of silvertip tetra were isolated from the caudal fin tissue of an adult male. Following the manufacturer’s instructions for the explant culture-based primary cell isolation kit, pCELL-Ex kit (CellQua, Seongnam, South Korea), the cells were isolated. The isolated cells were suspended in Leibovitz’s L-15 medium (Sigma-Aldrich, St. Louis, MO, USA), supplemented with 10% FBS (Thermo Fisher Scientific) and 1% Antibiotic-Antimycotic (GenDEPOT, Katy, TX, USA), and transferred to a T-25 cell culture flask (SPL, Pocheon, South Korea), then incubated at 26°C. The medium was replaced every three days, and the cell morphology was monitored using a fluorescence inverted microscope (Motic). When the primary cell confluency reached 70-80%, cells were subcultured using 0.25% trypsin-EDTA (Thermo Fisher Scientific). Continuous subculturing was performed for 3 to 5 days to develop a stable cell line. A total of 123 passages were performed, and the established cell line was designated as H. nana caudal fin, HNA FIN.

To determine the optimal culture temperature for the cell line, passage 115 cells were used to evaluate cell growth under different temperature conditions. The cells were dissociated using 0.25% Trypsin-EDTA and seeded at a concentration of 1×10⁴ cells/mL in a 24-well plate (Corning, Steuben, NY, USA). The cells were incubated for 5 days in L-15 medium supplemented with 10% FBS and 1% Antibiotic-Antimycotic at various temperatures: 20°C, 26°C, 32°C, and 38°C, and their viability was assessed. Cell viability was analyzed on days 3, 5, and 8 using the Quanti-Max WST-8 cell viability assay kit (Biomax, Guri, South Korea). Additionally, after culturing the cells at the optimal temperature of 32°C for 5 days, cell number was measured using the CellDrop-FL (DeNovix, Wilmington, DE, USA) to determine the doubling time.

The cryopreservation and post-thaw viability assay of the cell line were conducted using passage 111 HNA FIN cells at a concentration of 3×10⁶ cells/mL. CELLBANKER 2 (ZENOGEN PHARMA, Sakamachi, Japan) was used as the cryopreservation medium. The frozen cells were stored at -80°C for varying durations: passage 52 cells were stored for 474 days, passage 115 cells for 110 days, and passage 111 cells for 10 days. After storage, the cells were quickly thawed in a 37°C water bath, and the cryopreservation medium was removed by centrifugation. The cell viability after thawing was measured using the Max-View™ Live/Dead Staining Kit (Biomax), with the number of live cells determined by fluorescence using the CellDrop-FL. To confirm the attachment and proliferation of the thawed cells, they were suspended in L-15 medium supplemented with 10% FBS and 1% Antibiotic-Antimycotic and transferred to a 90-mm dish (SPL), followed by incubation at 32°C.

To assess the transfection efficiency of foreign genes in the cell line, the pDsRed2-1 plasmid vector expressing red fluorescent protein (Takara Bio, Kusatsu, Japan) was used. An expression vector with the beta-actin 1 promoter from the Fathead Minnow (Pimephales promelas) species was inserted into passage 60 cells. The experiment was performed following the manufacturer’s instructions for the FUGENE 4K Transfection Reagent (Promega, Madison, WI, USA). The optimal transfection reagent-to-vector ratio was determined, and the ideal ratio was found to be 5:1. Fluorescence expression in the transfected cells was observed using a fluorescence inverted microscope.

The whole-transcriptome next-generation sequencing (NGS) libraries were manufactured using Truseq stranded mRNA prep kit (Illumina, San Diego, CA, USA). NGS libraries were sequenced using NovaSeq 6000 (Illumina) on 101 x 2 stranded (paired-end) read module, and raw sequencing data were generated using base-calling software, bcl2fastq2 v2.2 and fastq quality score, ASCII Q-score (offset 33). NGS library adaptor sequences, low-quality sequences (limit=0.05), and ambiguous nucleotides (maximal two nucleotides) were removed using OmicsBox v3.3.0 (BioBam, Valencia, Spain). De novo assembly was constructed using trimmed pair-end reads since silvertip tetra genome assembly was not available for whole-transcriptome mapping reference. OmicsBox was used for generation of de novo assembly under following default parameters: word size: 25, bubble size 50, mismatch cost: 2, insertion cost: 3, deletion cost: 3, length fraction: 0.5, similarity fraction: 0.8, and minimum contig length 200 bp. Coding regions of the de novo assembled reads were predicted using TransDecoder v5.5.0 with universal genetic code and default parameters (Haas et al., 2013). Transcript redundancy was reduced and clustered biological sequences of assembled sequences using CD-HIT v4.6.4 with default parameters (sequence identity type: global, sequence identity threshold: 0.95) (Fu et al., 2012). Completeness assessment was performed using the Benchmarking Universal Single-Copy Orthologs (BUSCO) v4.0.5 with default parameters in Vertebrata lineage dataset (number of species: 3354, number of BUSCOs: 67) (Seppey et al., 2019).

Annotation of the de novo assembly contigs was carried out against the BLAST non-redundant protein sequences (NR) database in OmicsBox with the following conditions: E-value cut-off was 1.0E-3, the word size of BLAST parameters was 6, HSP length cut-off was 33, and HSP-Hit coverage was 0. In silico prediction was performed using InterProScan v5.45-80.0 (Finn et al., 2017), Gene Ontology (GO) (Götz et al., 2008), GO-Slim, and EggNOG-Mapper v1.0.3 with EggNOG v5.0.0 (Huerta-Cepas et al., 2019) for functional annotation. The functional annotation results were combined with those of the BLAST NR annotation to construct a whole-transcriptome reference of silvertip tetra, and the contigs from each sample were mapped to the whole-transcriptome reference assembly. Differentially expressed genes (DEGs) were investigated between females and males in each tissue (BRA, FIN, GON) with the FDR p-value below 0.05 using CLC Genomics Workbench v11.0. In addition, DEGs across various tissues were analyzed for relevant GO terms associated with pigmentation and sexual dimorphism using QuickGO (https://www.ebi.ac.uk/QuickGO/). The expression patterns of DEGs for functional analysis were classified based on fold change value based on up-regulation (fold change ≥ 2.0) and down-regulation (fold change ≤ -2.0). To analyze the expression of genes involved in pigmentation and identify their interactions, we employed STRING v11.0, utilizing the zebrafish (D. rerio) database with a minimum interaction score threshold set at > 0.5. This analysis allowed us to construct a comprehensive interaction network, incorporating multiple proteins. The resulting network was visualized using Cytoscape v3.8.1 software (Shannon et al., 2003).

qRT-PCR validation was performed on selected genes according to tissues and compared with the expression pattern of transcriptome data. Subsequent to normalize the concentration of total RNA in each tissue, cDNA was synthesized using HyperscriptTM RT Master Mix (Geneall) which included random hexamer and oligo (dT)₁₈. Additionally, to rectify the concentration of each sample in qRT-PCR, phosphatidylinositol-binding clathrin assembly protein-like isoform x3 (PICALM) were used as an internal control gene (Lee et al., 2024). The qRT-PCR primers were designed from whole-transcriptome reference and confirmed qRT-PCR amplification efficiency using serially diluted cDNA series acquired from a single sample ( Supplementary Table S1 ). qRT-PCR reactions were executed using QuantStudio™ 7 Flex Real-time PCR (Applied Biosystems™, Waltham, USA) with PowerUP™ SYBR™ Green Master Mix (Applied Biosystems) and gene-specific primers under the following cycling conditions: 2 steps of amplification (45 cycles), 95°C for 15 s and 60°C for 60 s. The Ct value of qRT-PCR was calculated on the delta-delta Ct method and expressed as mean ± SE and carried out by one-way ANOVA with a significance level of p < 0.05 using IBM SPSS 25.0 software (Livak and Schmittgen, 2001).

Information on embryonic development is essential in presenting the possibility as a new fish experiment model. Therefore, while securing the F1 generation of silvertip tetra, we tried to present information on the identified embryonic developmental stages and time. It took 22:10 hr from fertilization to hatching at 25°C ( Table 1 ). Fertilized eggs are small and are characterized by non-adhesive, and transparent chorion. 10 min after fertilization, the 1-cells were observed ( Figure 2A ). After 5 min of 2-cells stage, the second cleavage occurred and four blastomeres were shown ( Figure 2B ). Blastomeres occurred up to the morula stage at 1:15 hr ( Figure 2C ). The early blastula ( Figure 2D ) was observed at 1:40 hr after fertilization and at 2:40 hr the late blastulation process was shown ( Figure 2E ). The movement of yolk syncytial layer (epiboly) with the closure of the blastopore, 50% epiboly, 90% epiboly, and blastopore closure of the gastrula stage with 4:10 hr, 5:40 hr, 6:40 after fertilization were showed, respectively ( Figure 2F-H ). 10-myotomes were observed at 9:40 hr after fertilization ( Figure 2I ), the heartbeat has begun, tail became sharper and longer, with melanophore on the yolk at 12:40 hr ( Figure 2J ). After the tail reached to the head at 14:10 hr after fertilization ( Figure 2K ), embryos were fully formed at 15:40 hr ( Figure 2L ). Before hatching, there was an increase in embryo movement, and the embryo broke the chorion at 21:50 hr ( Figure 2M ) and hatched at 22:10 hr after fertilization ( Figure 2N ).

Primary cells cultured using the explant culture method were observed proliferating from the tissue by day 7 after seeding. In the early stages of culture, both fibroblast-like and epithelial-like cell morphologies were observed. By day 15 after seeding, the primary cells reached approximately 70-80% confluency on the surface of the T-25 flask, and subculturing was initiated. The cell line has been continuously subcultured for over 120 passages, and the cells have been morphologically confirmed as fibroblast-like cells ( Figure 3 ). After cryopreservation and thawing, the viability of HNA FIN cells was assessed by Live/Dead staining. The viability was 89.5 ± 2.7% for passage 52, 73.4 ± 2.5% for passage 111, and 78.8 ± 1.3% for passage 115. Additionally, after cell attachment, proliferation was confirmed, with cells reaching 70-80% confluency in the dish within 5 days of seeding. No morphological changes in the cells were observed after freezing and thawing.

To determine the optimal culture temperature for the cell line, cells were cultured at temperatures ranging from 20°C to 38°C, and the highest proliferation rate was observed at 32°C ( Figure 4A ). The growth rate of the cell line peaked on day 8 at 32°C. At 38°C, the cells showed a tendency to proliferate until day 3 after seeding, but the proliferation rate decreased from day 5 onward. To measure the doubling time of the cell line, the cells were cultured at the optimal temperature of 32°C for 5 days, and the cell count was used to calculate a doubling time of 28 hours ( Figure 4B ).

The observation of cells expressing fluorescent protein confirmed the successful transfection of the foreign expression vector. Fluorescent protein expression was detected 48 hours after transfection, and the number of expressing cells increased over time ( Figure 4C ).

A total of 246,799 contigs (160,853,660 bp) was generated from females and males of silvertip tetra with the Illumina Novaseq6000 platform ( Supplementary Table S2 ). The average length and N50 of assembled contigs was 652 bp and 1,012 bp. The result of the comparison with similar protein sequences using CD-HIT, a total of 245,900 clusters was identified at 95% identity. 17,243 contigs of the complete open reading frame (ORF), 7,509 contigs of the 5’-partial, 4,874 contigs of the 3’-partial, and 6,701 contigs of the internal ORF were predicted using TransDecoder 5.5.0. 2,779 contigs (81.63%) of the BUSCO groups have complete gene representation [single copy; 2,738 contigs ±(81.63%) and duplicated; 41 contigs (1.22%)], while 310 contigs (9.24%) were only partially recovered, and 265 contigs (7.9%) were missing in vertebrata_odb 10 lineage dataset. The high quality of the de novo assembly of silvertip tetra was verified based on the mapping rate of each tissue between 84.44% and 85.59% ( Supplementary Table S3 ).

A total of 48,876 contigs (19.8%) matched against the protein databases in BLASTx. Additional annotations were carried out against the EggNOG, GO, InterProScan, and KEGG database, resulting in a total of 37,154 contigs (15.1%), 20,465 contigs (8.3%), 54,881 contigs (22.2%), and 19,770 contigs (8.0%), respectively. 62% of the annotated genes matched the gene sequences from Astyanax mexicanus, which belongs to the same family (Characidae) as the silvertip tetra. Expression and clustering data were estimated using Principal Component Analysis (PCA) and Hierarchical clustering heat map for investigating the relationship between different gender and tissues ( Figure 5A ). As a result of PCA and heat map examination, each tissue-specific cluster was classified with different patterns according to gender within the tissue ( Figure 5B ). For gonad and caudal fin, the distribution varies along axis 1 (43.8% variability) and 2 (20.7% variability) between female and male, respectively. However, the brain-pituitary showed little difference in the distribution of axis 1 and axis 2 between female and male. In the heat map, the expression patterns of genes between female and male were similar broadly in the brain, whereas the expression patterns between genes were different as caudal fin and gonad were deviated clusters between female and male.

DEGs of each tissue were investigated for functional enrichment with GO and KEGG. The distribution of three GO domains (biological process; BP, molecular function; MF, and cellular component; CC) was characterized based on level 7 classification. For DEGs of each tissue, BP domains were represented by ‘nucleic acid-templated transcription’ and ‘regulation of nucleic acid-templated transcription’. MF domains mainly were comprised by ‘nucleoside-triphosphatase activity’ and ‘zinc ion binding’. In CC domains, ‘microtubule’ among GO terms were overwhelming occupied. Furthermore, several GO terms were identified by specific tissues or gender in all three GO domains ( Figure 6 ). Among the DEGs, several were identified as belonging to the Biological Process (BP) domain, specifically involved in pigmentation (GO:0043473), pigment metabolic process (GO:0042440), sex differentiation (GO:0007548), sex determination (GO:0007530), female sex differentiation (GO:0046660), and male sex differentiation (GO:0046661).

KEGG terms were dominated by ‘metabolic pathways’ in DEGs of BRA, FIN, and GON. Additionally, some KEGG terms have been annotated for specific tissue or gender, such as GO results. In particular, ‘Rap1 signaling pathway’, ‘Ras signaling pathway’, and ‘regulation of actin cytoskeleton’ mainly involved genes that showed higher expression in males than females in BRA, FIN, and GON ( Figure 7 ).

A Venn diagram was generated to investigate the types and distribution of DEGs in each tissue of females and males, with a maximum FDR p-value < 0.05 and a minimum absolute fold change > 1.5 ( Figure 8A ). A total of 49,605 DEGs were identified from the brain-pituitary, caudal fin, and gonad, comparing females and males. The distribution of DEGs across tissues revealed 3,715 DEGs in the brain, 28,053 in the caudal fin, and 39,931 in the gonad. Tissue-specific expression patterns were observed with 1,671 DEGs in the brain-pituitary, 7,600 in the caudal fin, and 19,159 in the gonad. Additionally, 919 DEGs were commonly expressed across the three tissues.

Among the DEGs showing common expression patterns across all three tissues, 13 DEGs in females and 19 in males exhibited a fold change ≥ 3 and an FDR p-value < 0.05, excluding undefined hypothetical proteins ( Figure 8B ).

In tissue-specific analyses, DEGs involved in reproduction, maturation, and growth were selected from the brain-pituitary, while DEGs involved in pigmentation and cation transport were identified in the caudal fin, and those related to gonadal development and sex differentiation were selected from the gonad.

In the brain-pituitary, the gene with higher expression in females was luteinizing hormone beta subunit (LHB), showing a fold change of 5.10. In contrast, the genes with higher expression in males were pituitary adenylate cyclase-activating polypeptide type I receptor-like isoform X4 (PACAP-R1 x4) with a fold change of -1.21, neuropeptide Y receptor type 1 (NPY1R) with -1.30, brain aromatase-like (BAA) with -1.33, insulin-like growth factor I (IGF-I) with -2.51, and prolactin-like (PRL) with -5.41 ( Figure 9A ).

In the caudal fin, the gene with higher expression in females was sodium/potassium/calcium exchanger 4 isoform X2 (NCKX4 x2), showing a fold change of 27.35, and beta-beta-carotene 15,15’-dioxygenase-like (BCO1) with a fold change of 7.73. In contrast, the genes with higher expression in males were melanocortin receptor (MC) with a fold change of -2.22 and melanocyte-stimulating hormone receptor (MSHR) with -2.92 ( Figure 9B ).

In the gonad, the gene with higher expression in females was growth/differentiation factor 9 (GDF9), showing a remarkable fold change of 482.68. Other genes with higher expression in females included doublesex- and mab-3-related transcription factor A2 (DMRT2) with a fold change of 23.41, follistatin isoform X1 (FST x1) with 9.58, and homeobox protein EMX2 (EMX2) with 9.50. In contrast, the gene with higher expression in males was aromatase (CYP19A1), showing a fold change of -123.47 ( Figure 9C ).

The interaction network of 16 proteins related to pigmentation was analyzed in the caudal fin. Among the 16 proteins, 3 were predominantly expressed in females, while the remaining 13 were predominantly expressed in males (Figure 10 ). The proteins with higher expression in females were identified as Endothelin-converting enzyme 2-like (ECE2), endothelin receptor type B-like (ednrb1a), and protein Mpv17 (mpv17), with fold change values of 10.63, 2.62, and 1.73, respectively. In contrast, the proteins with higher expression in males were melanocyte stimulating hormone receptor (mc1r), Tyrosinase-like (tyr), and microphtalmia-associated transcription factor-like isoform X1 (mitfa), with fold change values of -2.92, -3.44, and -6.36, respectively. Additionally, proteins such as melanocyte protein PMEL (pmela), membrane-associated transporter protein (slc45a2), transcription factor Sox-10 (sox10), and xanthine dehydrogenase/oxidase isoform x1 (xdh) exhibited higher expression in males. The protein-protein interaction (PPI) network analysis identified a total of 20 protein nodes and 60 edges. The average local clustering coefficient was 0.647, indicating a tendency for the proteins to form clusters within the network. Notably, the PPI enrichment p-value was less than 1.0e-16, strongly suggesting that these protein interactions were not random, but rather form a functionally significant protein network.

To validate the expression patterns of DEGs, qRT-PCR was performed on five genes selected from the BRA, FIN, and GON. In the brain-pituitary, the selected genes were LHB, Type II iodothyronine deiodinase (DIO2), IGF1, PRL, and Pro-neuropeptide Y (NPY). In the caudal fin, the selected genes were BCO1, PMEL, Dehydrogenase/reductase (SDR family) member (DHRS), Sodium/calcium/potassium exchanger 2 (NCKX2), and SRY (sex determining region Y)-box 7 (SOX7). In the gonads, the selected genes were Paired box 7 (PAX7), GDF9, (sex determining region Y)-box 9 (SOX9), Doublesex and mab-3 related transcription factor 1 (DMRT1), and Forkhead box L2 (FOXL2). qRT-PCR and RNA-seq results showed similar overall expression trends, despite variations in fold change among genes. Correlation analysis between qRT-PCR data and RNA-seq data revealed a correlation coefficient (R²) of 0.9097, indicating that 90.97% of the variation between the two datasets was consistent ( Figure 11 ).