Identification and Comparison of microRNAs in the Gonad of the Yellowfin Seabream (Acanthopagrus Latus)

Yellowfin seabream (Acanthopagrus latus) is a commercially important fish in Asian coastal waters. Although natural sex reversal has been described in yellowfin seabream, the mechanisms underlying sexual differentiation and gonadal development in this species remain unclear. MicroRNAs (miRNAs) have been shown to play crucial roles in gametogenesis and gonadal development. Here, two libraries of small RNAs, constructed from the testes and ovaries of yellowfin seabream, were sequenced. Across both gonads, we identified 324 conserved miRNAs and 92 novel miRNAs: 67 ovary-biased miRNAs, including the miR-200 families, the miR-29 families, miR-21, and miR-725; and 88 testis-biased miRNAs, including the let-7 families, the miR-10 families, miR-7, miR-9, and miR-202-3p. GO (Gene Ontology) annotations and KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment analyses of putative target genes indicated that many target genes were significantly enriched in the steroid biosynthesis pathway and in the reproductive process. Our integrated miRNA-mRNA analysis demonstrated a putative negatively correlated expression pattern in yellowfin seabream gonads. This study profiled the expression patterns of sex-biased miRNAs in yellowfin seabream gonads, and provided important molecular resources that will help to clarify the miRNA-mediated post-transcriptional regulation of sexual differentiation and gonadal development in this species.


Introduction
MicroRNAs (miRNAs) are small, endogenous, non-coding RNAs that are approximately 22 nucleotides long; miRNAs regulate gene silencing and translational repression by binding to target messenger RNAs (mRNAs) [1]. The relationships between miRNAs and mRNAs are intricate, in that one miRNA may regulate several genes, and each gene may also be targeted by multiple miRNAs [2]. Since miRNAs were first discovered in 1993, thousands have been identified in diverse animal groups, including mammals, birds, amphibians, and fish, and it has been shown that most mature miRNA sequences are conserved across these groups [3][4][5][6]. MiRNAs play important roles in the regulation of multiple biological processes, such as cell proliferation and differentiation, embryonic development, apoptosis, energy metabolism, immunity, metamorphosis, and gametogenesis [7][8][9][10][11][12][13].
Sexual development is a complex process that is tightly regulated by a series of molecules and signaling networks at transcriptional and post-transcriptional levels [14]. As high-throughput

Identification of Conserved and Novel miRNAs in the Gonads of the Yellowfin Seabream
By searching miRBase 22.0 and predicting novel miRNAs with miRDeep2, 324 known miRNAs and 92 novel miRNAs were identified in the yellowfin seabream gonads. Of these, 37 known miRNAs and 26 novel miRNAs were specific to the testis, whereas 51 known miRNAs and 5 novel miRNAs were specific to the ovary. The remaining 236 known miRNAs and 61 novel miRNAs were expressed in both the testis and the ovary (Figure 2A,B). In total, 360 mature miRNAs were identified in the testis, and 353 mature miRNAs were identified in the ovary. Other non-coding RNAs identified in the yellowfin seabream gonads, including rRNAs, tRNAs, snRNAs, snoRNAs, and repeat sequences, are shown in Table A3.

Identification of Conserved and Novel miRNAs in the Gonads of the Yellowfin Seabream
By searching miRBase 22.0 and predicting novel miRNAs with miRDeep2, 324 known miRNAs and 92 novel miRNAs were identified in the yellowfin seabream gonads. Of these, 37 known miRNAs and 26 novel miRNAs were specific to the testis, whereas 51 known miRNAs and 5 novel miRNAs were specific to the ovary. The remaining 236 known miRNAs and 61 novel miRNAs were expressed in both the testis and the ovary (Figure 2A,B). In total, 360 mature miRNAs were identified in the testis, and 353 mature miRNAs were identified in the ovary. Other non-coding RNAs identified in the yellowfin seabream gonads, including rRNAs, tRNAs, snRNAs, snoRNAs, and repeat sequences, are shown in Table A3.

Identification of Conserved and Novel miRNAs in the Gonads of the Yellowfin Seabream
By searching miRBase 22.0 and predicting novel miRNAs with miRDeep2, 324 known miRNAs and 92 novel miRNAs were identified in the yellowfin seabream gonads. Of these, 37 known miRNAs and 26 novel miRNAs were specific to the testis, whereas 51 known miRNAs and 5 novel miRNAs were specific to the ovary. The remaining 236 known miRNAs and 61 novel miRNAs were expressed in both the testis and the ovary (Figure 2A,B). In total, 360 mature miRNAs were identified in the testis, and 353 mature miRNAs were identified in the ovary. Other non-coding RNAs identified in the yellowfin seabream gonads, including rRNAs, tRNAs, snRNAs, snoRNAs, and repeat sequences, are shown in Table A3.
The 10 most abundant miRNAs in the testis were miR-143, miR-30d, miR-15b, and seven novel miRNAs, whereas the 10 most abundant miRNAs in the ovary were miR-143, miR-30d, miR-181c, and 7 novel miRNAs. Of these, miR-143, miR-30d, miR-novel-46, miR-novel-47, and miR-novel-70 were highly abundant in both the testis and the ovary of the yellowfin seabream (Table 1). The highly abundant miRNAs in both testis and ovary are highlighted in red.

Identification of Sex-Biased miRNAs in the Testis and the Ovary
We identified 155 miRNAs that were differentially expressed between the testes and the ovaries of yellowfin seabream ( Figure 3): 88 testis-biased miRNAs (41 known miRNAs and 47 novel miRNAs) and 67 ovary-biased miRNAs (55 known miRNAs and 12 novel miRNAs). Of these, 24 miRNAs were exclusively expressed in the testes and 6 miRNAs were specifically expressed in the ovaries. Excluding these gonad-specific miRNAs, the greatest fold-change in expression level between the testis and the ovary was a 92.5-fold change in testis-biased miRNA expression and a 38.7-fold change in ovary-biased miRNA expression. In addition, the expression levels of 20 testis-biased miRNAs (31.25%) were over 8-fold greater than those in the ovary, whereas the expression levels of 10 ovary-biased miRNAs (16.39%) were over 8-fold greater than those in the testis. Thus, more testis-biased miRNAs were differentially expressed than ovary-biased miRNAs, and the fold-changes in expression levels were greater. Overall, the let-7 families, the miR-10 families, miR-7, miR-9, and miR-202-3p were highly expressed in the testis, whereas the miR-200 families, the miR-29 families, miR-21, and miR-725 were upregulated in the ovary. The 30 most DE miRNAs between the testis and the ovary are shown in

Target Gene Prediction and Functional Analysis
To investigate the potential functions of the DE miRNAs, we predicted the corresponding target genes. Our analyses predicted 16,217 target genes of the DE miRNAs, which were classified into 57 GO (Gene Ontology) subcategories: 26 biological process terms, 16 cellular component terms, and 15 molecular function terms. Of the biological process terms, the target genes were mainly associated with cellular process (GO: 0009987), single-organism process (GO: 0044699), and metabolic process (GO: 0008152). Of the cellular component terms, the target genes were mainly associated with membrane (GO: 0016020), membrane part (GO: 0044425), and cell (GO: 0005623). Of the molecular function terms, the target genes were mainly associated with binding (GO: 0005488), catalytic activity (GO: 0003824), and transporter activity (GO: 0005215) ( Figure 6). KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analysis indicated that the target genes were annotated in 131 signaling pathways. The 20 most enriched KEGG pathways are shown in Figure 7. Notably, eight pathways were significantly enriched (p < 0.05), including the phosphatidylinositol signaling system (ko04070), endocytosis (ko04144), and steroid biosynthesis (ko00100).
Next, we analyzed several important GO terms and KEGG pathways associated with sexual development and reproduction enriched in the putative target genes, including the GO term "reproductive process" and the steroid biosynthesis signaling pathway. We identified 61 potential target genes associated with reproductive processes, including zona pellucida glycoprotein 3 (zp3), sperm acrosome membrane-associated protein 4 (spaca4), annexin A1 (anxa1), mitotic-specific cyclin-B1 (ccnb1), and meiotic recombination protein dmc1. In addition, 17 putative target genes were enriched in the steroid biosynthesis signaling pathway, including 17β-hydroxysteroid

Target Gene Prediction and Functional Analysis
To investigate the potential functions of the DE miRNAs, we predicted the corresponding target genes. Our analyses predicted 16,217 target genes of the DE miRNAs, which were classified into 57 GO (Gene Ontology) subcategories: 26 biological process terms, 16 cellular component terms, and 15 molecular function terms. Of the biological process terms, the target genes were mainly associated with cellular process (GO: 0009987), single-organism process (GO: 0044699), and metabolic process (GO: 0008152). Of the cellular component terms, the target genes were mainly associated with membrane (GO: 0016020), membrane part (GO: 0044425), and cell (GO: 0005623). Of the molecular function terms, the target genes were mainly associated with binding (GO: 0005488), catalytic activity (GO: 0003824), and transporter activity (GO: 0005215) ( Figure 6). KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analysis indicated that the target genes were annotated in 131 signaling pathways. The 20 most enriched KEGG pathways are shown in Figure 7. Notably, eight pathways were significantly enriched (p < 0.05), including the phosphatidylinositol signaling system (ko04070), endocytosis (ko04144), and steroid biosynthesis (ko00100).

Integrated Analysis of the DE miRNAs and Corresponding Target Genes
To understand the relationships between DE miRNAs and the corresponding sex-related target genes, we constructed regulatory networks (Figure 8). MiR-202-3p was predicted to target spaca4, anxa1, ccnb1, sqle, and sperm flagellar protein 1; these genes are associated with gametogenesis, cellular signaling, and hormonal secretion [30][31][32][33][34]. The testis-biased miRNA let-7j was predicted to target cyp3c and lhb, two ovary-biased genes that may be associated with steroids synthesis and the sex hormone responsible for follicular development [35]. The ovary-biased miR-29b might target cyp11b1, which encodes a rate-limiting enzyme of steroid hormones synthesis that is required for testis development [36]. MiR-29b may also target meiotic recombination protein dmc1, which is essential for spermatogenesis [37]. cellular signaling, and hormonal secretion [30][31][32][33][34]. The testis-biased miRNA let-7j was predicted to target cyp3c and lhb, two ovary-biased genes that may be associated with steroids synthesis and the sex hormone responsible for follicular development [35]. The ovary-biased miR-29b might target cyp11b1, which encodes a rate-limiting enzyme of steroid hormones synthesis that is required for testis development [36]. MiR-29b may also target meiotic recombination protein dmc1, which is essential for spermatogenesis [37].
The overall distribution of small RNAs differed between the testis and the ovary. We observed both miRNA and piRNA peaks in the ovary, but only one piRNA peak in the testis. PiRNAs were more plentiful in the testis, whereas miRNAs were more enriched in the ovary. The most abundant small RNAs in the gonads were the piRNAs, which play critical roles in gametogenesis, especially testis development [43,44]. During the long term of oogenesis, many materials, including mRNAs, proteins, and miRNAs, are stored in the egg, as preparation for early embryonic development.
To better compare the similarities and differences of the miRNA expression profiles in Acanthopagrus latus and six reported fishes. Except for the species-specific miRNAs, we summarized nine testisor ovary-biased miRNAs, including let-7, miR-10, miR-7, miR-9, miR-202, miR-200, miR-29, miR-21, and miR-725, which exhibited significant sex-biased expression in at least four fish species (Table 2) [18][19][20][21][22][25][26][27]. MiR-7 exhibited a significant testis-biased expression in Acanthopagrus latus and three other fish species, indicating that it might play an important role in testis development. However, the other eight miRNAs exhibited different gender biases in different fish species. For example, let-7 was highly expressed in the testis of Acanthopagrus latus, Pelteobagrus fulvidraco, Acipenser schrenckii, and Odontobutis potamophila, whereas it was identified as an ovary-biased miRNA in Oplegnathus punctatus; miR-9 was a testis-biased miRNA in Acanthopagrus latus, Pelteobagrus fulvidraco, Acipenser schrenckii, and Trachinotus ovatus, but was mainly expressed in Oreochromis niloticus. This result indicates that they may play different roles in gonadal development in different species. Table 2. Nine important tests-or ovary-biased miRNAs in yellowfin seabream and six reported fishes.
It has been reported that miR-143 is highly expressed during gonadal development, and plays an important role in this process. For example, miR-143 inhibits steroid hormone synthesis and granulosa cell apoptosis by targeting FSHR (follicle stimulating hormone receptor) in the bovine ovary [45]. In mice, miR-143 regulates granulosa cell proliferation by targeting the FSH signaling pathway [46]. The miR-30 family was highly expressed in the testes of Trachinotus ovatus [26], yellow catfish [19], and Amur sturgeon [19,25]. The miR-181 family is highly expressed in the gonads of tilapia [21], mice [47], and humans [48]. In yellowfin seabream gonads, miR-143, miR-30d, and miR-181c were the most highly expressed miRNAs, indicating that these might participate in gonadal development in the yellowfin seabream, similar to what has previously been shown in other animals.
In mice, the let-7 family was expressed in the spermatogonia and spermatocytes, and was associated with spermatogonial differentiation through the retinoic acid signaling pathway [49].
The let-7 miRNA families were also identified as testis-biased in Amur sturgeon [25], Atlantic halibut [50], and yellow catfish [19]. However, the expression patterns of the let-7 families differed in spotted knifejaw, where the let-7 families were more highly expressed in the ovary than in the testis [22]. In mice, miR-7a2 affected sexual maturation and reproductive processes by regulating the synthesis and secretion of gonadotropic and sex-steroid hormones; notably, the deletion of miR-7a2 caused infertility in both sexes [51]. MiR-7 was also highly expressed in yellow catfish [18], suggesting that this miRNA might participate in testis development in fish. It has been suggested that MiR-9 might regulate sexual development in several species. For example, miR-9 was upregulated in the tilapia ovary and was predicted to affect sexual determination and differentiation by regulating the expression of dmrt1 [20]. In contrast, in Monopterus albus, miR-9 might affect the physiological processes that promote oocyte degeneration and might stimulate spermatogenesis via foxl3 expression during natural sex change [52]. Here, the let-7 families, miR-7, and miR-9 exhibited testis-biased expression patterns, suggesting that these miRNAs might participate in spermatogenesis and testis development in yellowfin seabream.
Some miRNAs were also strongly expressed in the yellowfin seabream ovary, including the miR-200 family (miR-200a, miR-200b, and miR-429a), the miR-29 family (miR-29a, miR-29b, and miR-29c), and miR-21. In mice, miR-200b and miR-429 inhibited luteinizing hormone synthesis by negatively regulating zeb1, which is required for ovulation and female fertility [53]. In zebrafish, miR-200 clusters regulate sperm motility by targeting amh, wt1a, and srd5a2b [41]. The miR-200 family is also highly expressed in the ovaries of several other fish species [18,22]. In mice, the miR-29 family might participate in male meiosis by downregulating Dnmt3a, which plays an important role in the response to DNA damage [54], while miR-29-3p and its target Tbx21 help to regulate the onset of puberty and reproduction in female mice by modulating Gnrh1 expression [55]. In addition, miR-29 is an ovary-biased miRNA in tilapia [21]. MiR-21 is required to block the apoptosis of granulosa cells during mouse periovulation [56]. Moreover, miR-21 regulates ovarian development in rainbow trout [57]. Thus, the ovary-biased expression patterns of the miR-200 family, the miR-29 family, and miR-21 suggested that these miRNAs might be involved in ovarian development in yellowfin seabream.
As miRNAs themselves cannot encode functional protein products, these small RNAs affect biological function by the post-transcriptional regulation of target genes. The predicted target genes of the DE miRNAs were classified into 57 GO subcategories. Of these, 61 putative target genes were enriched in the GO term "reproductive process", suggesting that the corresponding DE miRNAs were involved in yellowfin seabream reproduction. In addition, the putative target genes were significantly enriched in eight KEGG pathways, including "steroid biosynthesis", "endocytosis", and "phosphatidylinositol signaling system". Sex steroids are tightly associated with sexual differentiation, gonadal development, and sexual reversal [58][59][60]. Endocytosis is required for the cellular uptake of sex steroids and yolk during oogenesis [61,62]. The phosphatidylinositol signaling system is primarily associated with the phosphatidylinositol-3 kinase (PI3K) pathway, which is critical for gonadal development [63]. Therefore, the DE miRNAs might critically affect gonadal development by targeting their gonad-development-associated target genes.
Each miRNA may interact with many different target genes, and each coding gene may be targeted by many miRNAs. MiR-202-3p was predicted to target spaca4, anxa1, ccnb1, sqle, and sperm flagellar protein 1. Conversely, cyp11b1 was predicted to be the target of miR-29b, miR-181a, miR-199-3p, and miR-novel-64. Therefore, miRNAs and mRNAs may form complicated interaction networks. Here, we constructed integrated networks of the DE miRNAs and known genes associated with sexual differentiation. These networks demonstrated negatively correlated expression patterns: ovary-biased DE miRNAs always targeted genes associated with testis development, while the testis-biased DE miRNAs always targeted genes associated with ovarian development. The testis-determining gene sox9 was targeted by three ovary-biased miRNAs (miR-124, miR-novel-86, and miR-novel-87), which suggests that inactivating sox9 expression might potentially maintain ovarian development in yellowfin seabream. Indeed, it has been demonstrated that miR-124 prevents sox9 expression during ovarian development [64]. Meanwhile, the oocyte-polarity marker gene foxh1 was targeted by miR-103, miR-107a-3p, and miR-novel-49; three miRNAs that were significantly downregulated in the yellowfin seabream ovary. Therefore, repression of miR-103, miR-107a-3p, and miR-novel-49 might be required for oogenesis in yellowfin seabream. The miRNA-mRNA interaction network provides significant information for the identification of miRNAs associated with gonadal development. The direct interactions between miRNAs and their predicted targets, and its role in gonadal development will be analyzed in our future studies.

Ethical Procedures
All procedures with Acanthopagrus latus were approved by the Ethics Committee of Sun Yat-Sen University (protocol no. 20200110008, approval date: 11 October 2019) and the methods were carried out following the approved guidelines.

Sample Collection
Yellowfin seabreams were obtained from Hailv aquaculture station (Zhuhai, China). Six yellowfin seabreams were selected: three 1-year-old males, with an average body weight (ABW) of 113.3 ± 8.5 g and an average body length (ABL) of 15.3 ± 0.3 cm, and three 3-year-old females, with an ABW of 519.3 ± 55.3 g and an ABL of 26.1 ± 0.8 cm. Each fish was anesthetized using 100 mg/L tricaine methane sulfonate (MS-222). Fish gonads were extracted, snap frozen in liquid nitrogen, and stored at −80 • C for RNA extraction.

RNA Extraction, Library Construction, and Small RNA Sequencing
Samples were removed from the −80 • C freezer and homogenized for 40 seconds with Bioprep-24 (Tomos, Shanghai, China) in 1 mL TRIzol reagent (Invitrogen, Carlsbad, CA, USA). Total RNAs were extracted according to the manufacturer's instructions. The extracted RNA was further treated with DNaseI (Vazyme, Nanjing, China) to eliminate genomic DNA contamination. RNA degradation and contamination were further monitored by 1.5% agarose gel. The quantity and quality of the isolated RNA were measured using a BioSpec-nano Spectrophotometer (Shimadzu, Japan), and the RNA integrities were assessed by Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Only RNAs with an OD260/280 ≥ 1.8 and an RIN (RNA integrity number) ≥ 7 were retained. The total RNA of each sample was standardized to 200 ng/µL, then an equal amount of high-quality RNA (10 µL) from each sample was pooled (testes and ovaries were pooled separately), and sent to Majorbio Bio-pharm Technology Co., Ltd. (Shanghai, China) for small RNA library construction and sequencing. Briefly, RNA was purified by PAGE to enrich for 18-32 nt molecules, ligated with 5 -and 3 -adaptors, and was then sequenced using the Illumina HiSeq platform.

Differential miRNA Expression Analysis
MiRNA expression levels were standardized based on the number of transcripts per million reads (TPM). Differentially expressed miRNAs (DE miRNAs) between the testes and the ovaries were identified using DEGSeq Version 1.30.0 (http://bioconductor.org/packages/stats/bioc/DEGSeq/), correcting the false discovery rate with the Benjamini-Hochberg method. MiRNAs with an absolute fold-change value ≥ 2 and an adjusted p-value < 0.001 were identified as significant DE miRNAs.

Real-Time Quantitative PCR Verification (RT-qPCR)
To validate the transcriptome data, a total of 10 miRNAs were screened for RT-qPCR, including 5 testis-biased miRNAs and 5 ovary-biased miRNAs. The RNAs used for RT-qPCR verification were the same as that for transcriptome sequencing. Briefly, testis and ovary tissues were dissected from three one-year-old and three-year-old individuals, respectively. Total RNA was isolated from each sample using TRIzol reagent (Invitrogen, USA). The extracted RNA from three individuals in each group (testes and ovaries) were mixed in equal mass ratios, and then reverse transcribed using Mir-X miRNA First-Strand Synthesis Kits (Takara, Dalian, China), following the manufacturer's instructions. Roche Light Cycler 480 (Roche, Switzerland) was applied to perform RT-qPCR with TB Green Premix Ex Taq (Takara, Dalian, China). The thermal cycling conditions were as follows: 94 • C for 5 min; 40 cycles of 94 • C for 12 s, 58 • C for 12 s, and 72 • C for 10 s; and an additional 72 • C for 2 min. After thermal cycling, a final disassociation curve analysis was performed. Each sample was analyzed in triplicate, and miRNA relative expression levels were normalized to those of U6 snRNA using the 2 −∆∆CT method. Statistical analysis was performed by SPSS 20.0. Differences with p values < 0.05 were considered statistically significant. The primers used in this study are listed in Table A1.

Prediction of DE miRNA Targets
The potential target genes of the DE miRNAs were predicted using miRanda (http://www.miranda. org/), Targetscan (http://www.targetscan.org/), and RNAhybrid (http://bibiserv.techfak.uni-bielefeld. de/rnahybrid/). The 3 UTRs (untranslated region) of all of the genes were extracted from the yellowfin seabream gonad transcriptome (SRA accession number: PRJNA622226), and matched to the miRNAs. Gene Ontology (GO) annotations and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed to identify the functions and pathways associated with these candidate target genes.

Conclusions
In summary, the small-RNA profiles of the yellowfin seabream testis and ovary were obtained. We identified 324 conserved known miRNAs and 92 novel miRNAs. Of these, 155 sex-biased DE miRNAs were significantly enriched in sex-associated pathways or GO terms, such as the steroid biosynthesis pathway and reproductive processes. Our integrated miRNA-mRNA analysis demonstrated that sex-biased DE miRNAs and known sexual differentiation-associated genes formed an intricate signaling pathway that regulated sexual differentiation in the yellowfin seabream. Our data provide important molecular resources that help to clarify the roles and mechanisms by which miRNAs affect gonadal development in the yellowfin seabream.

Acknowledgments:
We thank LetPub (www.letpub.com) for its linguistic assistance during the preparation of this manuscript.

Conflicts of Interest:
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Appendix A