Identification of Sexually Dimorphic Genes in Pectoral Fin as Molecular Markers for Assessing the Sex of Japanese Silver Eels (Anguilla japonica)

© 2023 Academia Sinica, Taiwan Open Access Identification of Sexually Dimorphic Genes in Pectoral Fin as Molecular Markers for Assessing the Sex of Japanese Silver Eels (Anguilla japonica) Hsiang-Yi Hsu1, Chia-Hsien Chuang2, I-Hsuan Lu2, Chung-Yen Lin2, and Yu-San Han1,* 1Institute of Fisheries Science, College of Life Science, National Taiwan University, 1, Sec. 4, Roosevelt Rd., Taipei 10617, Taiwan. *Correspondence: E-mail: [email protected] (Han). E-mail: [email protected] (Hsu) 2Institute of Information Science, Academia Sinica, Taipei 115, Taiwan. E-mail: [email protected] (Chuang); [email protected] (Lu); [email protected] (Lin) Received 2 June 2022 / Accepted 25 October 2022 / Published 10 March 2023 Communicated by John Wang The Japanese eel (Anguilla japonica) is an important species in East Asian aquaculture. However, the production of seedlings for this purpose still depends on natural resources, as the commercial production of glass eels is not yet possible. Confusion about the sex of silver eels is one of the factors affecting the success rate of artificial maturation. This study sought to devise a harmless method to precisely assess the sex of silver eels. Partial pectoral fins were collected from females and males and the total RNA was extracted for transcriptomic analysis to identify sexually dimorphic genes as molecular markers for sex typing. An online database was constructed to integrate the annotations of transcripts and perform comparative transcriptome analysis. This analysis identified a total of 29 candidate sexually dimorphic genes. Ten were selected for a real-time quantitative polymerase chain reaction (RT-qPCR) to validate the transcriptomic data and evaluate their feasibility as markers. The transcriptomic analysis and RT-qPCR data implicated three potential markers (LOC111853410, kera, and dcn) in sex typing. The expression of LOC111853410 was higher in females than in males. In contrast, the expression of kera and dcn was higher in males than in females. The ΔCT values of three markers were analyzed to determine their inferred thresholds, which can be used to determine the sex of Japanese eels. The results suggested that if a silver eel had a pectoral fin with the pectoral fin having the ΔCT of LOC111853410 < 11.3, the ΔCT of kera > 11.4, or the ΔCT of dcn > 6.5 can be assessed it could be assessed as female. Males could be assessed by the ΔCT of LOC111853410 > 11.3, the ΔCT of kera < 11.4, or the ΔCT of dcn < 6.5 in their pectoral fins. The molecular functions of these markers and the biological significance of their differential expression require further exploration. Key words: Comparative transcriptome analysis, Japanese eel, Pectoral fin, Sex typing, Sexually dimorphic genes. Citation: Hsu HY, Chuang CH, Lu IH, Lin CY, Han YS. 2023. Identification of sexually dimorphic genes in pectoral fin as molecular markers for assessing the sex of Japanese silver eels (Anguilla japonica). Zool Stud 62:02. doi:10.6620/ZS.2023.62-02. BACKGROUND The Japanese eel, Anguilla japonica, is a catadromous fish that spends most of its growth phase in freshwater habitats but spawns in the open ocean. Its spawning site is located in the west of the Mariana Ridge, as confirmed by the eggs and leptocephali collected in the water column (Tsukamoto 1992 2006). In recent years, overfishing, habitat destruction, and changes in the ocean environment have led to the Zoological Studies 62: 2 (2023) doi:10.6620/ZS.2023.62-02 1 © 2023 Academia Sinica, Taiwan drastic depletion of natural resources (Chen et al. 2014). Consequently, the Japanese eel has been listed on the International Union for Conservation of Nature Red List of Threatened Species as an endangered species since 2014 (Pike et al. 2020). Conserving the natural resources of the Japanese eel and ensuring a stable supply of seedlings for eel culture requires the development of techniques for artificial propagation. The artificial propagation of Japanese eels has been studied for a long time (Okamura et al. 2014; Tanaka 2015). However, techniques for the mass production of glass eels have not yet been established. The reasons include unclear mechanisms of sex determination and differentiation, unrecognizable sexes for silver eels before artificial maturation, imperfect techniques of artificial maturation and an unstable supply of parent eels, unstable quality of fertilized eggs, unintelligible mechanisms of nutrient utilization and metabolism of eel larvae, and inappropriate feed and cultivation methods for eel larvae (Hsu et al. 2015; Okamura et al. 2014 2018 2020; Tanaka 2015; Tesch 2003; Yamada et al. 2019). An increasing number of studies are addressing these difficulties, with the goal of establishing techniques for the mass production of glass eels. Karyological analyses of anguillid eels have suggested that they possess heteromorphic sex chromosomes (Park and Kang 1979; Park and Grimm 1981). However, the presence of heteromorphic sex chromosomes is not a determinant of sex determination and differentiation in eels (Wiberg 1983). Sex determination and differentiation of eels begins at the juvenile (yellow eel) stage (Colombo and Grandi 1996). Environmental factors, such as water temperature, population density, and body size of juvenile eels, have been suggested to override the genetic effects on sex determination and differentiation (Colombo et al. 1984; Colombo and Grandi 1996; Davey and Jellyman 2005; Oliveira and McCleave 2000; Tesch 2003). Higher water temperature and population density, and larger body size of juvenile eels may increase the proportion of male eels. Sex differentiation in eels has been associated with the expression of several sexually dimorphic genes/hormones in the brain and gonads (Geffroy et al. 2016; Inaba et al. 2021; Jeng et al. 2018). The expression of amh, dmrt1, gsdf, sox9a, and vasa significantly increased during testicular differentiation of eels. The expression of cyp19a1, figla, foxl2a, foxn5, sox3, zar1, and zp3 was high during the ovarian differentiation of eels. Treatment with estrogens had a feminization effect on fish, including eels (Chiba et al. 1993; Colombo and Grandi 1995; Guiguen et al. 2010). Observations have indicated a higher proportion of females in wild eel populations compared to that in farmed eels due to the low population density. However, the detailed mechanisms controlling the expressions of sexually dimorphic genes remain unclear (Chu et al. 2006; Han and Tzeng 2006). Eel biologists have speculated that appearance characteristics, such as asymptotic length, eye size, and shape of the pectoral fin, vary between sexes (Tesch 2003). However, even after long-term observation, we found it difficult to precisely identify the sex of Japanese silver eels that completed sex differentiation based on their appearance characteristics, which may be affected by genetics, age, and environment. In a previous study on the reproduction of European eels, du Colombier et al. (2015) introduced ultrasonography as a noninvasive tool to assess the sex and maturity of European silver eels. The accuracy of this technique was insufficient. Additionally, a recent report mentioned that performing anesthesia surgery to sample a piece of gonadal tissue for histological analysis can determine the sex of Japanese silver eels, but it involves risks of causing injury, inflammation, and even death to the eels (Huang et al. 2020). Therefore, the development of a new technique for assessing the sex of silver eels with no harm and high accuracy at the same time is essential and would also be helpful for the artificial propagation of eels. Because of the discrepancy in the procedures of artificial maturation for females and males, such as the required time and hormone combination, knowing the sex of silver eels beforehand would improve the efficiency, success rate, and costs of artificial maturation (Okamura et al. 2014; Tanaka 2015). Subsequently, controlling the sex ratio of each batch of parent eels would assist researchers in arranging the start point of artificial maturation for each eel and the mating time among parent eels, which would increase the success rate of fertilization and the number of fertilized eggs to obtain a large number of leptocephali for cultivation (Okamura et al. 2020). This study aimed to devise a method for precisely and harmlessly assessing the sex of silver eels before artificial maturation. In addition to the brain and gonads, sexually dimorphic genes may also exist in other body tissues of the eels. Thus, a piece of the pectoral fins of two female and two male silver eels was collected for RNA-sequencing (RNA-seq). A comparative transcriptome analysis identified candidate sexually dimorphic genes, which may serve as markers for sex typing in Japanese silver eels. RT-qPCR experiments validated the transcriptomic data and assessed the feasibility of these genes as molecular markers. If there are available molecular markers for sex typing, it will improve the artificial maturation of Japanese eels. Moreover, introducing a similar research strategy to discover the molecular markers for assessing the page 2 of 16Zoological Studies 62: 2 (2023) © 2023 Academia Sinica, Taiwan gonadal maturity of Japanese eels during artificial maturation can also be possible. MATERIALS AND METHODS Sample collection Ten Japanese silver eels comprising five females (body weight 675–900 g) and five males (body weight 460–575 g) were purchased from a local eel farmer in Pingtung in southern Taiwan. On the basis of our previous study (Han et al. 2003), silver eels (adult stage) could be differentiated from yellow ones (juvenile stage) by their morphometric characteristics, such as skin coloration, total length, ocular index (OI), and fin-index (FI). The eels were acclimated to seawater in a two-ton tank at the laboratory for approximately two weeks for subsequent artificial maturation. Samples of pectoral fins were collected from all ten anesthetized silver eels before the artificial maturation. The biopsy samples were immediately stored in a -80°C freezer. Biopsy samples of two females and two males were randomly selected for RNA-seq, and the remaining biopsies were reserved for subsequent RT-qPCR analysis. Additional pectoral fin tissue biopsied from other silver eels was sampled for RT-qPCR analysis. Two females (body weights of 658 g and 705 g) and two males (body weights of 489 g and 618 g) were purchased from a commercial eel dealer in Tainan, southern Taiwan. Two female eels (body weights of 420 g and 560 g) were reared from elvers and fed an E2-containing diet (Chiba et al. 1993; Colombo and Grandi 1995; Guiguen et al. 2010) at the culture station of National Taiwan University for 24 months. Pre-sampling handling of silver eels and post-sampling storage of biopsies was performed as described above. Detailed information on all the silver eels is shown in table S1. In addition, ten cultured eels at the juvenile stage were bought from an eel farmer in Pingtung to sample the pectoral fins for RT-qPCR analysis. All of them were male, and no female juvenile eel was obtained in this study. These juvenile eels were reared in a twoton tank with freshwater for three days after arriving at the laboratory to ease their stress. Subsequently, they were sacrificed for the sampling and measurement of morphometric characteristics. Detailed information on these juvenile eels is displayed in table S2. All procedures and investigations were approved by the Institutional Animal Care and Use Committee of the National Taiwan University (NTU107-EL-00059) and executed in accordance with standard guidelines. RNA extraction, library construction, and sequencing Total RNA from the entire mass of the biopsied tissue was extracted using TRIzol® Reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. Purified RNA was quantified using an ND1000 spectrophotometer (Nanodrop, Wilmington, DE, USA) and characterized using a Bioanalyzer 2100 with an RNA 6000 Labchip kit (Agilent Technologies, Santa Clara, CA, USA). After the Bioanalyzer 2100 analysis, the RNA integrity numbers of all RNA samples prepared in this study exceeded 7.0. Sequencing libraries were constructed using the TruSeq RNA Sample Prep Kits v2 (Illumina, San Diego, CA, USA) according to the manufacturer’s instructions. They were subsequently sequenced using the Illumina HiSeq-2500 platform. Processing of sequence data, gene prediction, and gene annotation Four datasets included two female pectoral fins (FPF1: ERS13509836 and FPF2: ERS13509837) and two male pectoral fins (MPF1: ERS13509838 and MPF2: ERS13509839), which have been deposited in European Nucleotide Archive (ENA) of European Bioinformatics Institute (EBI). Raw RNA-seq data were filtered using the TrimGalore program (Babraham Bioinformatics, Cambridge, UK) to discard adaptors and low-quality reads (Q < 13). Low-complexity reads (repeat sequences) were removed using the prinSeq program. Finally, the general read properties were generated using the FastQC program (Babraham Bioinformatics). The MAKER2 pipeline (https://www.yandelllab.org/software/maker.html) was used to predict genes and generate the gene transfer format (GTF). The input datasets included the four RNA-seq datasets obtained in this study, previously published RNA-seq data (Hsu et al. 2015), PacBio Iso-Seq data, newly established Japanese eel genome (unpublished data), available protein sequences of A. japonica, all teleosts from the UniProt database, and the pre-existing zebrafish gene model. The new version of the Japanese eel genome established by our team had 1.06 Gb with 4,861 scaffolds and N50 over 4.6 Mb. By taking 448 scaffolds with lengths larger than 100 Kb, these scaffolds can cover 98% of the Japanese eel genome. The predicted open reading frames were annotated with homologs in the updated NR database (April 2019) using GhostX (Suzuki et al. 2014). Protein sequence features, including signal peptides, transmembrane domains, and domains described in the Pfam database, were detected using SignalP (Petersen et al. 2011), TMHMM page 3 of 16Zoological Studies 62: 2 (2023) © 2023 Academia Sinica, Taiwan (Krogh et al. 2001), and HMMER (Mistry et al. 2013), respectively. Gene ontology annotation was performed for genes with detectable Pfam domains according to Pfam2GO. Protein-coding genes were mapped to the canonical Kyoto Encyclopedia of Genes and Genomes (KEGG) database pathway using the KEGG Automatic Annotation Server (Moriya et al. 2007). The parameters were set to map to only related eukaryotic databases in the single-directional best-hit mode. Construction of database framework and estimation of gene expression An online database (MOLAS) for the Japanese eel was constructed using the LAMP system architecture (Ubuntu 14.04, Apache 2.04, PHP 5.1 and MySQL 8.0) with the Bootstrap 3 CSS framework (http:// getbootstrap.com/), jQuery1.11.1, and jQuery Validation v1.17 to provide an intuitive user experience. The entire system runs on a virtual machine on the cloud infrastructure of the Institute of Information Science, Academia Sinica, Taiwan. The analytical process was implemented using scripts written in R (3.4.2). According to the assembly with GTF, the gene expression profiles can be estimated using raw reads generated from RNA-seq via an intuitive graphical interface in Docexpress (https://hub.docker.com/r/lsbnb/ docexpress_fastqc) with a built-in process of Hisat2 -> StringTie -> Ballgown. Finally, the gene expression profiles in fragments per kilobase of exon per million fragments mapped (FPKM) values were submitted to MOLAS implemented by our team. Comparative analysis between female and male pectoral fin libraries for gene expression profiles In this study, the “Pairwise Comparison” function on the online database was first used to compare the datasets of the female and male pectoral fins to find the differentially expressed genes (highly expressed in females but minimally in males, and vice versa). Two female pectoral fin datasets (FPF1 and FPF2) and two male pectoral fin datasets (MPF1 and MPF2) were individually pooled for this analysis. According to previous studies on sexually dimorphic genes in Japanese eels (Geffroy et al. 2016; Inaba et al. 2021; Jeng et al. 2018), the expression levels of most sexually dimorphic genes were significantly different between female and male eels, which was over five fold-change. Therefore, screening parameters were set such that the expression levels (FPKM value) of transcripts were ≥ 0, and the differences in expression levels of transcripts were at least six-fold between the two pools. Then, the “Clustering” function was used to categorize the expression patterns of the transcripts identified from the above step and remove the transcripts that were highly expressed in only one dataset of female or male pool and that were minimally expressed in four datasets of female and male pool. Thereafter, the transcripts with higher expression in females and males were separately filtered out to identify candidate genes for assessing the sex of Japanese eel. RT-qPCR The specific primer pairs of candidate genes for RT-qPCR were designed using the Primer 3 software and listed in table 1. The specificity of each primer pair was verified using PCR and electrophoresis, and amplicons were validated using DNA sequencing. In addition, the qPCR efficiency of each primer pair was tested by serially diluting cDNA samples from eel pectoral fin, as listed in table 1. The cDNA of each specimen was reverse transcribed from 1 μg of total RNA with four repeats using iScript cDNA Synthesis kits (Bio-Rad, Hercules, CA, USA), according to the manufacturer’s protocol. RT-qPCR was performed using a MyiQ realtime PCR system (Bio-Rad) with a 25 μL mixture containing 12.5 μL of 2× SYBR Green Supermix (Bionova, Fremont, CA, USA), 1 μL of each primer, 1 μL of cDNA template (10× dilution), and 9.5 μL of distilled deionized water. The running conditions followed were pre-incubation at 95°C for 10 min and 40 amplification cycles of 95°C for 30 s, 58°C for 45 s, and 72°C for 45 s. Acidic ribosomal phosphoprotein P0 (arp), which is consistently expressed at the same level in eels, was used as a housekeeping gene (Sudo et al. 2013). RT-qPCR data were analyzed in two steps. First, arp was used for normalization by subtracting the cycle threshold (CT) value from the CT value of the target gene in every sample to get the ΔCT. Second, in each gene category, two female samples (F1 and F2) were used as the control for normalization by subtracting the average ΔCT value of them from the ΔCT value of each sample to get the ΔΔCT. The relative mRNA expression levels were calculated as 2-(ΔΔCT) and were presented as the Vertical Box Plot with two groups (S-Female and S-Male) and the Vertical Bar Chart with all the samples in each gene category. Besides, in another analysis for defining the criteria for sex typing of silver eels, the ΔCT (normalization with arp) values of candidate sexually dimorphic genes were compared between the S-Female and S-Male groups. Then, the threshold of ΔCT of each candidate gene was inferred by summing the lowest ΔCT of the group having higher ΔCT values with the highest ΔCT of the other group and taking the average. The ΔCT values of every candidate gene in the S-Female and page 4 of 16Zoological Studies 62: 2 (2023) © 2023 Academia Sinica, Taiwan S-Male groups were presented as the Vertical Point Plot. Statistical analyses For the statistical test on the differentially expressed transcripts that have six fold-change between the female pool (FPF) and male pool (MPF), we applied the R package, limma, to perform linear modeling on the expression levels (FPKM) and used the empirical Bayes method to moderate gene-wise variability, generating P values (Ritchie et al. 2015). The P values were adjusted by the Benjamini-Hochberg (BH) method (Benjamini and Hochberg 1995). The expression values of ten genes in all samples from RT-qPCR and the ΔCT values of three candidate sexually dimorphic genes were statistically analyzed using the SAS version 9.1 statistical software (SAS Institute Inc., Cary, NC, USA). The expression values of ten genes in all samples from RT-qPCR were analyzed using one-way ANOVA followed by a Student-Newman-Keuls multiple comparison test, with P < 0.05 indicating a significant difference. Student’s t-test evaluated the differences in the expression values of ten genes between the S-Female and S-Male groups in the Vertical Box Plot and the differences in the ΔCT values of three candidate sexually dimorphic genes between the S-Female and S-Male groups in Vertical Point Plot, with P < 0.05 indicating a significant difference. RESULTS Next-generation sequencing statistics and establishment of online database An average of approximately 11.54 million raw reads, excluding those for ribosomal RNA and mitochondrial RNA, of nearly 150 base pairs (bp) in length were obtained for each library (see Table S3). The raw data of each library were cleaned up and mapped to the Japanese eel genome (see Table S3), followed by gene prediction, functional annotation, and estimation of expression levels of transcripts. A gene model containing 27,350 Japanese eel transcripts was constructed. The transcripts were blasted against the NR database, and a total of 27,336 best hits were identified. In addition, the annotated information and expression levels of all transcripts were obtained after a series of functional annotations and calculations and integrated into the MOLAS (https://molas.iis.sinica.edu.tw/2022_ jpeel/login.php) database established by our research team (Fig. 1). Table 1. Specific primer pairs used for real-time quantitative PCR analyses (F, forward strand; R, reverse strand) Transcript ID Hit_seq name Gene name on NCBI Primer sequences Amplicon size (bp) qPCR efficiency AJA4564.p1 septin-2-like isoform X2 LOC108249696 F: 5'- AACAGGAACAGGGAAGTCTG -3' R: 5'- TCGGAGCTGATCTTGTCTTC -3' 84 99.9% AJA12577.p1 heat shock protein 30-like LOC110515957 F: 5'- TGAGAAGACCCTGCAGTTCC -3' R: 5'- TGCTGGTGATTCTGTGCTTC -3' 69 99.8% AJA24273.p1 early growth response protein 3 egr3 F: 5'- TCCTGGGAAGCGGGAGAG -3' R: 5'- CCAGAGAAGTACCACGGTCG -3' 137 99.5% AJA16680.p1 CLOCK-interacting pacemaker cipc F: 5'- GAGTCAGAGAGGGACTCTGG -3' R: 5'- GTAGGACCCCGTCATCACC -3' 136 99.6% AJA21056.p1 mucin-19-like LOC111853410 F: 5'- GCAACTGTGAAGGCAACAAG -3' R: 5'- ACTGTTACAAGGTACCACGC -3' 128 99.6% AJA20026.p1 src substrate cortactin cttn F: 5'- ATGTGTACCAGTCGGAGCC -3' R: 5'- CACTGCCGAACCTCCTCTAC -3' 229 98.8% AJA14877.p1 keratocan isoform X1 kera F: 5'- TCCCAATGCTCGTACTTTCG -3' R: 5'- TGCTTCAGGTTCCTGTTGTC -3' 214 99.1% AJA20844.p1 decorin precursor-like dcn F: 5'- ACAACAAGCTCCAGACCATC -3' R: 5'- AGGTGCCTCAGTTTGGAGTA -3' 139 99.5% AJA7588.p1 cbp/p300-interacting transactivator 1 cited1 F: 5'- GCATCATCGATTCGGACCC -3' R: 5'- CTCAGCAATGGTCAAACGGA -3' 188 99.4% AJA1691.p1 amine sulfotransferase-like LOC104575574 F: 5'- GTTCTTGGAGAGAGGCATGG -3' R: 5'- TCGCTGATCTTCACAACGAC -3' 159 99.4% -arp F: 5'- GTGCCAGCTCAGAACACTG -3' R: 5'- ACATCGCTCAAGACTTCAATGG -3' 107 99.8% page 5 of 16Zoological Studies 62: 2 (2023) © 2023 Academia Sinica, Taiwan Comparative analysis between female and male pectoral fin libraries for gene expression profiles The “Pairwise Comparison” and “Clustering” functions on the database were used to compare the libraries of pectoral fins of female and male eels to discover differentially expressed transcripts. After parameter setting and operation of the platform, a total of 15 transcripts displayed high expression levels in the pectoral fins of female eels but low or almost no expression in the pectoral fins of male eels. The detailed Fig. 1. Scheme of the MOLAS database for the Japanese eel. Annotated information of the assembled transcripts, including DNA sequence, amino acids sequence, open reading frames, signal peptides, transmembrane domains, gene ontology, hit KEGG pathway, and FPKM values, can be searched by the names of genes or the transcript ID in the Full-text search. A sequence of DNA can also be used to find the transcript with high similarity through the Sequence Search/BLAST. Moreover, the Pairwise Comparison can compare the differences in expression levels of the transcripts between two different libraries or two library groups, and then the differentially expressed genes can be summarized to a gene list. Furthermore, the Import Genelist can analyze the protein function, gene ontology enrichment, heatmap of expression, and hit terms on a KEGG pathway for a gene list. Additionally, the Clustering can be used to categorize the expressional patterns of transcripts between two different libraries or two library groups. Finally, the KEGG GlobalView and Gene List Analysis can see the hit terms of transcripts on the map of KEGG pathway and compare the different gene lists by Venn diagrams, respectively. page 6 of 16 Zoological Studies 62: 2 (2023) © 2023 Academia Sinica, Taiwan annotations of these transcripts are listed in table 2. Fourteen other transcripts displayed higher expression levels in the pectoral fins of male eels but low or almost no expression levels in the pectoral fins of females. The detailed annotations of the 14 transcripts are presented in table 3. These twenty-nine differentially expressed transcripts all pass the statistical test (P < 0.05) which implies that they indeed have significant differences in expression levels between the female pool and male pool. The adjusted P values are shown in table 2 and table 3. The sequence information of these twenty-nine transcripts is presented in table S4. These differentially expressed transcripts may be sexually dimorphic genes in the pectoral fin of the Japanese eel. Subsequently, ten of these transcripts were selected for RT-qPCR to validate the transcriptomic data and assess their feasibility as molecular markers. The expression levels of six candidates in the pectoral fins of female eels were higher than those in the pectoral fins of males. On the other hand, the expression levels of four candidates were higher in the pectoral fins of male eels than in the pectoral fins of females. RT-qPCR of pectoral fin tissues from female and male silver eels Ten candidate genes obtained from pectoral fin tissues of seven female (F1‒F7) and five male (M1‒ Table 2. Information of transcripts whose expression in the libraries of pectoral fins of female silver eels was six-fold higher than that of male silver eels Transcript ID CDS length NR best hit Best hit_E value Hit_seq name AJA12763.p1 453 nt KKF25601.1 4.689e-15 Protein NLRC3 AJA1330.p1 1221 nt XP_018604518.1 1.651e-141 galactosylceramide sulfotransferase-like AJA14095.p1 791 nt XP_018924236.1 1.715e-60 protein fosB-like isoform X4 AJA15563.p1 1080 nt XP_018594957.1 1.641e-155 fatty acid 2-hydroxylase AJA16680.p1 462 nt XP_023684314.1 5.780e-37 CLOCK-interacting pacemaker AJA18689.p1 1398 nt XP_020340178.1 1.494e-128 otopetrin-1 AJA20026.p1 345 nt XP_021180512.1 1.277e-20 src substrate cortactin AJA21056.p1 508 nt XP_023686008.1 1.384e-35 mucin-19-like AJA24273.p1 2355 nt XP_018599663.1 4.190e-210 early growth response protein 3 AJA25664.p1 489 nt XP_023663962.1 2.182e-52 class I histocompatibility antigen, F10 alpha chain-like AJA12577.p1 2210 nt XP_021450334.1 1.637e-64 heat shock protein 30-like AJA27270.p1 267 nt KPP63944.1 9.285e-27 Annexin A11-like, partial AJA4564.p1 396 nt XP_017294726.1 2.513e-08 septin-2-like isoform X2 AJA4608.p1 1263 nt XP_017579547.1 1.058e-106 tyrosine-protein kinase fyna AJA737.p1 1071 nt OWK00640.1 1.469e-23 MBNL2, partial Transcript ID Hit_organism Involved function Adjusted P value from limma R package AJA12763.p1 Larimichthys crocea Ribonuclease inhibitor 3.76e-05 AJA1330.p1 Scleropages formosus Galactosylceramide sulfotransferase activity 5.25e-05 AJA14095.p1 Cyprinus carpio DNA-binding transcription factor 7.70e-06 AJA15563.p1 Scleropages formosus Lipid biosynthetic process 2.64e-05 AJA16680.p1 Paramormyrops kingsleyae Negative regulation of transcription 0.000274 AJA18689.p1 Oncorhynchus kisutch Formation of otoconia and proton channels 0.000258 AJA20026.p1 Fundulus heteroclitus Actin filament polymerization 7.95e-06 AJA21056.p1 Paramormyrops kingsleyae Formation of mucin protein 0.000215 AJA24273.p1 Scleropages formosus DNA binding 2.87e-05 AJA25664.p1 Paramormyrops kingsleyae MHC class I antigen 1.11e-06 AJA12577.p1 Oncorhynchus mykiss Structural constituent of ribosome and phospholipid binding 2.57e-06 AJA27270.p1 Scleropages formosus Calcium-dependent phospholipid binding 0.000653 AJA4564.p1 Kryptolebias marmoratus Transmembrane transporter activity, GTP-binding activity, and innate immunity 0.002781 AJA4608.p1 Pygocentrus nattereri Tyrosine-protein kinase activity 2.56e-06 AJA737.p1 Cervus elaphus hippelaphus RNA splicing protein 4.44e-05 page 7 of 16Zoological Studies 62: 2 (2023) © 2023 Academia Sinica, Taiwan M5) silver eels were subjected to RT-qPCR. The results of RT-qPCR displayed that the expression levels of cipc, LOC111853410, and cttn in the female silver eel group were statistically significantly higher than those in the male silver eel group (Fig. 2D‒F). However, there were no significant differences in the expression of cipc or cttn among some female and male silver eel samples (Fig. S1D and Fig. S1F). In addition, the expression levels of LOC108249696, LOC110515957, and egr3 had no significant differences between the female and male silver eel groups, as shown in figure 2A–C. On the other hand, the RT-qPCR showed that the expression levels of kera, dcn, cited1, and LOC104575574 were significantly higher in the male silver eel group than in the female silver eel group (Fig. 3). However, there were no significant differences in the expression of cited1 or LOC104575574 among some female and male silver eel samples (Fig. S2C and Fig. S2D). Detailed expression levels of ten candidate genes in the pectoral fin tissues from all silver eel samples are displayed in figures S1 and S2. These RT-qPCR results indicated that the expression of LOC111853410, kera, and dcn is consistent with the transcriptomic data and suggested that only LOC111853410, kera, and dcn may be potential sexually dimorphic genes used for sex typing of silver eels. Possible criteria for sex typing of silver eels The best potential candidate sexually dimorphic genes, LOC111853410, kera, and dcn, were selected after analyses of transcriptomic data and RT-qPCR. The ΔCT values of these three candidate genes from the pectoral fins of silver eels were individually analyzed to Table 3. Information of transcripts whose expression in the libraries of pectoral fins of male silver eels was six-fold higher than that of female silver eels Transcript ID CDS length NR best hit Best hit_E-value Hit_seq name AJA11404.p1 495 nt XP_016828489.1 5.355e-35 peptidyl-prolyl cis-trans isomerase FKBP11 isoform X5 AJA13052.p1 822 nt BAP63925.1 1.385e-149 type I collagen proa2 chain AJA14877.p1 1933 nt XP_003445032.2 2.642e-132 keratocan isoform X1 AJA1532.p1 420 nt XP_012692196.1 2.741e-63 collagen alpha-1(X) chain-like AJA1691.p1 702 nt XP_010221213.1 3.972e-55 amine sulfotransferase-like AJA20844.p1 1038 nt KPP76630.1 1.960e-134 decorin precursor-like AJA21947.p1 707 nt XP_006629547.2 4.719e-125 secreted frizzled-related protein 2 AJA22857.p1 582 nt XP_023655487.1 5.314e-81 mimecan AJA26246.p1 297 nt XP_015457250.1 3.073e-43 SH3 domain-binding glutamic acid-rich-like protein AJA2914.p1 327 nt ACN78878.1 4.287e-53 glutathione peroxidase 1 AJA461.p1 372 nt XP_012680047.1 8.861e-22 2-oxoisovalerate dehydrogenase beta subunit, mitochondrial AJA5036.p1 360 nt ELW64932.1 3.881e-41 DCN1-like protein 5 AJA7588.p1 921 nt XP_007242076.1 2.743e-54 cbp/p300-interacting transactivator 1 AJA8438.p1 416 nt XP_020796713.1 9.756e-13 Cytochrome c1, heme protein, mitochondrial Transcript ID Hit_organism Involved function Adjusted P value from limma R package AJA11404.p1 Cricetulus griseus Prolyl isomerase activity 0.000213 AJA13052.p1 Anguilla japonica Formation of type I collagen 4.52e-08 AJA14877.p1 Oreochromis niloticus Formation of keratocan 2.26e-06 AJA1532.p1 Clupea harengus Formation of type X collagen 4.20e-07 AJA1691.p1 Tinamus guttatus Amine sulfotransferase activity 0.000215 AJA20844.p1 Scleropages formosus Formation of an extracellular matrix protein 5.01e-08 AJA21947.p1 Lepisosteus oculatus Inhibition of osteoblast differentiation 2.33e-05 AJA22857.p1 Paramormyrops kingsleyae As a osteoglycin involved in the regulation of energy homeostasis 3.25e-06 AJA26246.p1 Astyanax mexicanus Electron transfer activity 0.000792 AJA2914.p1 Anguilla japonica Response to oxidative stress and glutathione peroxidase activity 2.02e-06 AJA461.p1 Clupea harengus Formation of beta subunit of 2-oxoisovalerate dehydrogenase 0.000747 AJA5036.p1 Tupaia chinensis Formation of an extracellular matrix protein 0.000464 AJA7588.p1 Astyanax mexicanus Regulation of transcription, DNA-templated 1.00e-05 AJA8438.p1 Boleophthalmus pectinirostris Electron transfer activity 8.31e-05 page 8 of 16Zoological Studies 62: 2 (2023) © 2023 Academia Sinica, Taiwan establish the criteria for sex typing. The Vertical Point Plot of every candidate gene shows that the ΔCT values of the S-Female and S-Male groups are clearly separated (Fig. 4). Moreover, there was a significant difference (P < 0.05) in the ΔCT values between the S-Female and S-Male groups for every candidate gene. For the gene of LOC111853410, the lowest ΔCT of the S-Male group and the highest ΔCT of the S-Female group were summed and averaged to infer the threshold of ΔCT for sex typing (Fig. 4A). In addition, for the genes of kera and dcn, the lowest ΔCT of the S-Female group and the highest ΔCT of the S-Male group were summed and averaged to infer the threshold of ΔCT (Fig. 4B and 4C). The inferred thresholds of ΔCT of LOC111853410, kera, and dcn were 11.3, 11.4, and 6.5, respectively. The analytic results indicated that the pectoral fins of female silver eels had the ΔCT of LOC111853410 < 11.3, the ΔCT of kera > 11.4, and the ΔCT of dcn > 6.5. Fig. 2. Comparisons of relative expression levels of six genes possibly higher expressed in the pectoral fins of female silver eels between the “S-Female” (n = 7, four repeats) and “S-Male” (n = 5, four repeats) groups. (A) LOC108249696, (B) LOC110515957, (C) egr3, (D) cipc, (E) LOC111853410, (F) cttn. The arp was a housekeeping gene, and two female samples (F1 and F2) were used as the internal control for calculating relative expression levels. The y-axis was in log2 scale. The “S-Female” meant the female silver eels, and the “S-Male” meant the male silver eels. The solid black line represented the median value, and the solid black circle represented the outliers. Asterisks indicated significant differences (**P < 0.01, ***P < 0.001) in relative expression levels between two groups. page 9 of 16Zoological Studies 62: 2 (2023) © 2023 Academia Sinica, Taiwan Sharba S, Sundh H, Sundell K, Benktander J, Santos L, Birchenough G et al. 2022. Rainbow trout gastrointestinal mucus, mucin production, mucin glycosylation and response to lipopolysaccharide. Fish Shellfish Immun 122:181–190. doi:10.1016/j.fsi.2022.01.031. Sudo R, Suetake H, Suzuki Y, Aoyama J, Tsukamoto K. 2013. Profiles of mRNA expression for prolactin, growth hormone, and somatolactin in Japanese eels, Anguilla japonica: The effect of salinity, silvering and seasonal change. Comp Biochem Physiol A Mol Integr Physiol 164:10–16. doi:10.1016/j.cbpa.2012.09.019. Suzuki S, Kakuta M, Ishida T, Akiyama Y. 2014. GHOSTX: An improved sequence homology search algorithm using a query suffix array and a database suffix array. PLoS ONE 9:e0103833. doi:10.1371/journal.pone.0103833. Takagi Y, Ura K. 2007. Teleost fish scales: A unique biological model for the fabrication of materials for corneal stroma regeneration. J Nanosci Nanotechno 7:757–762. doi:10.1166/jnn.2007.503. Tanaka H. 2015. Progression in artificial seedling production of Japanese eel Anguilla japonica. Fisheries Sci 81:11–19. doi:10.1007/s12562-014-0821-z. Tesch FW. 2003. Developmental stages and distribution of the eel species. In: J. E. Thorpe (ed.) The eel No. 2, pp. 73–118. Tsukamoto K. 1992. Discovery of the spawning area for Japanese eel. Nature 356:789–791. doi:10.1038/356789a0. Tsukamoto K. 2006. Spawning of eels near a seamount. Nature 439:929–929. doi:10.1038/439929a. Weis SM, Zimmerman SD, Shah M, Covell JW, Omens JH, Ross J et al. 2005. A role for decorin in the remodeling of myocardial infarction. Matrix Biol 24:313–324. doi:10.1016/ j.matbio.2005.05.003. Wiberg UH. 1983. Sex determination in the European eel (Anguilla anguilla, L.). A hypothesis based on cytogenetic results, correlated with the findings of skewed sex ratios in eel culture ponds. Cytogenet Cell Genet 36:589–598. doi:10.1159/000131981. Winnemoller M, Schon P, Vischer P, Kresse H. 1992. Interactions between thrombospondin and the small proteoglycan decorin - interference with cell attachment. Eur J Cell Biol 59:47–55. Yamada Y, Okamura A, Mikawa N, Horie N, Tsukamoto K. 2019. A new liquid-type diet for leptocephali in mass production of artificial glass eels. Fisheries Sci 85:545–551. doi:10.1007/ s12562-019-01295-2. Yang VWC, LaBrenz SR, Rosenberg LC, McQuillan D, Hook M. 1999. Decorin is a Zn2+ metalloprotein. J Biol Chem 274:12454– 12460. doi:10.1074/jbc.274.18.12454. Yeh LK, Liu CY, Chien CL, Converse RL, Kao WWY, Chen MS et al. 2008. Molecular analysis and characterization of zebrafish keratocan (zKera) gene. J Biol Chem 283:506–517. doi:10.1074/ jbc.M707656200. Zoeller JJ, Pimtong W, Corby H, Goldoni S, Iozzo AE, Owens RT et al. 2009. A central role for decorin during vertebrate convergent extension. J Biol Chem 284:11728–11737. doi:10.1074/jbc. M808991200. Supplementary materials Table S1. Detailed information of Japanese silver eels used in this study. (download) Table S2. Detailed information of Japanese yellow eels used in this study. (download) Table S3. Data of NGS statistics and mapping rate to the genome of Japanese eel for each library. (download) Table S4. Coding DNA sequences of the 29 transcripts which are candidate sexually dimorphic genes in the pectoral fin of Japanese eel. (download) Fig. S1. Expression profiles of six genes, possibly having higher expression in the pectoral fin tissues of female silver eels, were obtained from silver and yellow eel samples by RT-qPCR. (A) LOC108249696, (B) LOC110515957, (C) egr3, (D) cipc, (E) LOC111853410, (F) cttn. The arp was a housekeeping gene, and two female samples (F1 and F2) were used as the internal control for calculating relative expression levels. The data were presented as means ± SD, and the y-axis was in log2 scale. Different letters indicate significant differences (P < 0.05) in gene expression among the pectoral fin tissues of Japanese eel samples. For the silver eels, F meant the female, and M meant the male. Additionally, Y-M meant the male yellow eels. (download) Fig. S2. Expression profiles of four genes, possibly having higher expression in the pectoral fin tissues of male silver eels, were obtained from silver and yellow eel samples by RT-qPCR. (A) kera, (B) dcn, (C) cited1, (D) LOC104575574. The arp was a housekeeping gene, and two female samples (F1 and F2) were used as the internal control for calculating relative expression levels. The data were presented as means ± SD, and the y-axis was in log2 scale. Different letters indicate significant differences (P < 0.05) in gene expression among the pectoral fin tissues of Japanese eel samples. For the silver eels, F meant the female, and M meant the male. Additionally, Y-M meant the male yellow eels. (download) Fig. S3. Comparisons of ΔCT values of LOC111853410, kera, and dcn between the “S-Male” (n = 5, four repeats) and “Y-Male” (n = 10, four repeats) groups for assessing whether the inferred threshold from silver eels is applicable to juvenile eels (yellow eels) for sex typing. (A) LOC111853410, (B) kera, (C) dcn. A solid black circle represented one ΔCT value. The “S-Male” meant the male silver eels, and the “Y-Male” meant the male yellow eels. For each gene, there were 20 and 40 solid black circles in the “S-Male” group and the “Y-Male” group, respectively. The solid red line meant the inferred threshold of ΔCT from silver eels for sex typing. (download) page 16 of 16Zoological Studies 62: 2 (2023)