Brain and Pituitary Transcriptome Analyses Reveal the Differential Regulation of Reproduction-Related LncRNAs and mRNAs in Cynoglossus semilaevis

Long noncoding RNAs (lncRNAs) have been identified to be involved in half-smooth tongue sole (Cynoglossus semilaevis) reproduction. However, studies of their roles in reproduction have focused mainly on the ovary, and their expression patterns and potential roles in the brain and pituitary are unclear. Thus, to explore the mRNAs and lncRNAs that are closely associated with reproduction in the brain and pituitary, we collected tongue sole brain and pituitary tissues at three stages for RNA sequencing (RNA-seq), the 5,135 and 5,630 differentially expressed (DE) mRNAs and 378 and 532 DE lncRNAs were identified in the brain and pituitary, respectively. The RNA-seq results were verified by RT-qPCR. Moreover, enrichment analyses were performed to analyze the functions of DE mRNAs and lncRNAs. Interestingly, their involvement in pathways related to metabolism, signal transduction and endocrine signaling was revealed. LncRNA-target gene interaction networks were constructed based on antisense, cis and trans regulatory mechanisms. Moreover, we constructed competing endogenous RNA (ceRNA) networks. In summary, this study provides mRNA and lncRNA expression profiles in the brain and pituitary to understand the molecular mechanisms regulating tongue sole reproduction.


INTRODUCTION
Half-smooth tongue sole is one of the most economically important marine flatfish that is widely cultured in China. However, its reproductive dysfunctions and 3 years required for sexual maturation (Shi et al., 2016;Song et al., 2020) restrict the development of tongue sole aquaculture industry, which influences the economic benefit. Therefore, understanding of the molecular mechanism of reproduction in this species is scientifically and commercially necessary.
Reproduction is a crucial function of the organism and is controlled by complex interactions between the hypothalamus, pituitary and gonads, i.e., the hypothalamic-pituitary-gonadal (HPG) axis (Abreu and Kaiser, 2016). Within this axis, the hypothalamus releases the gonadotrophin releasing hormone (GnRH), which can act on pituitary gonadotropes, triggering the release of follicle-stimulating hormones (FSH) and luteinizing hormones (LH). FSH and LH bind to their respective receptors in the ovary, leading to secretion of the sex steroid hormones, estrogen (E 2 ), which stimulates vitellogenesis, and progesterone (P4), which stimulates oocyte meiosis, follicular maturation and ovulation (Yaron and Levavi-Sivan, 2011). Although the importance of the HPG axis for reproduction has been suggested, in contrast to the wealth of knowledge in mammals, the integrated molecular mechanism of HPG axis regulation in nonmammalian species has not been elucidated.
With the rapid development of RNA sequencing (RNA-Seq), the function of noncoding RNAs (ncRNAs), which were originally considered "junk RNAs" of the mammalian genome, has received increasing attention (Kopp and Mendell, 2018). Long noncoding RNAs (lncRNAs) are a group of transcripts that are longer than 200 nucleotides, without apparent proteincoding role (Quinn and Chang, 2016) and, similar to mRNAs, often contain a poly(A) tail and can be spliced (Hombach and Kretz, 2016). LncRNAs have been reported to play important roles in diverse biological processes, such as signal transduction (Arun et al., 2012), cell differentiation and development (Fatica and Bozzoni, 2014), ontogeny (Abdelmohsen et al., 2013) and immune responses (Chen et al., 2017), through transcriptional and posttranscriptional mechanisms (Dykes and Emanueli, 2017) that regulate gene expression, such as epigenetic modification (Mercer and Mattick, 2013), gene imprinting (Williamson et al., 2011), miRNA sponging (Faghihi et al., 2010) and chromatin remodeling (Tsai et al., 2010). Some studies have shown that lncRNAs are involved in the nervous system function and neurological diseases in humans (Qureshi et al., 2010;Oe et al., 2019). In addition, emerging evidence indicates that lncRNAs are involved in reproduction-related processes, including sex hormone responses (Yang et al., 2013), oocyte meiosis , and ovulation (Lian et al., 2020), although these studies focused mainly on humans and other mammals.
Recently, an increasing number of studies on lncRNAs have been reported in teleosts, focusing on processes such as growth and development (Ali et al., 2018;Wu et al., 2020), stress responses (Dettleff et al., 2020;Quan et al., 2020), immune responses (Liu et al., 2019;Zheng et al., 2021) and sex differentiation (Cai et al., 2019;Feng et al., 2020). However, the related knowledge still lags far behind that in mammals. Our previous research has provided lncRNA expression profiles of the tongue sole ovary, and we investigated the potential role of lncRNAs in the ovary in reproduction (Dong et al., 2021). However, in the HPG axis, the brain and pituitary are also crucial tissues in reproductive regulation. The hypothalamus, an area of the brain, can released a series of hormones that act on the pituitary, which is not only regulated by the hypothalamus but also influences the ovarian function through hormones and other regulatory factors, thus playing a connecting role in the HPG axis. However, lncRNAs in the HPG axis of tongue sole, particularly in brain and pituitary, that regulate reproduction in this fish have not been systematically analyzed.
Thus, it is necessary to identify and characterize mRNAs and lncRNAs in the reproductive neuroendocrine system to investigate their potential function in ovarian development, maturation and ovulation in tongue sole. Therefore, in this study, the brain and pituitary of tongue sole in three stages of ovarian development [stage IV: late vitellogenesis, stage V: maturation and stage VI: after ovulation (Shi et al., 2016)] were collected for sequencing to get their transcriptome profiles and explore the effect of the brain and pituitary on transcriptomic mechanisms controlling reproduction in tongue sole. This research expands the mRNA and lncRNA catalog in the tongue sole brain and pituitary and preliminarily elucidates the molecular mechanisms in which key lncRNAs and mRNAs are involved. Furthermore, our findings reveal candidate regulators for improving the efficiency of tongue sole reproduction.

Sample Preparation and Ethics Statement
We previously explained the feeding and sampling process of experimental fish (Dong et al., 2021). In brief, we collected triplicate brain and pituitary samples from tongue sole in three stages of ovarian development (stages IV, V and VI). The identification of ovarian developmental stages has been described in detailed in our previous studies (Dong et al., 2021). Tissues were immediately frozen in liquid nitrogen and stored at −80°C to prevent degradation before total RNA isolation.
The acquisition of the fish was approved by the Animal Care and Use Committee of the Chinese Academy of Fishery Sciences, and all experiments were conducted at the Chinese Academy of Fishery Sciences in accordance with the Guidelines for the Care and Use of Laboratory Animals at the Chinese Academy of Fishery Sciences, China.

RNA Isolation, Library Preparation and RNA Sequencing
According to the manufacturer's procedure, total RNA was isolated using a TRIzol reagent kit (Invitrogen, Carlsbad, CA, United States). The quality of the RNA was measured via an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, United States) and RNase free agarose gel electrophoresis. The rRNAs were detached from the total RNA with a Ribo-Zero TM Magnetic Kit (Epicentre, Madison, WI, United States) to retain mRNAs and ncRNAs. After rRNA detached, the enriched mRNA and ncRNAs were broken into short pieces by using fragmentation buffer and were then reverse transcribed into cDNA with random primers from an RNA Library Prep Kit (NEB, United States). Before the end repair, poly (A) addition, and Illumina sequencing adapter ligation, QIAquick PCR Extraction Kit (Qiagen, Venlo, Netherlands) was used to purify the cDNA fragments after second-strand cDNA synthesis. After the digestion of second-strand cDNA with uracil-N-glycosylase (UNG), the products were filtered by agarose gel electrophoresis, PCR amplification and sequenced on the Illumina HiSeq TM 4000 platform by Gene Denovo Biotechnology Co., Ltd. (Guangzhou, China).

Sequencing Data Processing and Analysis
Raw reads were produced by sequencing and filtered with fastp (version 0.18.0) (Chen et al., 2018) removing low-quality reads. Short reads were mapped to the rRNA database with Bowtie2 (version 2.2.8) (Langmead and Salzberg, 2012). After removing mapped rRNA reads, the remaining clean reads were further aligned to the tongue sole reference genome (Chen et al., 2014) by using TopHat2 (version 2.1.1). Rebuilding and identification of transcripts was conducted with Cufflinks (v2.1.1), TopHat2 (Kim et al., 2013), Cuffmerge (Trapnell et al., 2012), and Cuffcompare. We used the Coding-Non-Coding-Index (CNCI, version 2) (Sun et al., 2013), Coding Potential Calculator (CPC, http://cpc.cbi. pku.edu.cn/) (Kong et al., 2007) and the SwissProt protein database to distinguish the protein-coding capability of the new transcripts. The intersections of the identified transcripts without protein-coding potential were described as the collection of lncRNAs, which could be classified to 5 categories: intergenic lncRNAs, bidirectional lncRNAs, intronic lncRNAs, antisense lncRNAs, and sense-overlapping lncRNAs.

Differentially Expressed Transcripts and Functional Enrichment Analysis
We used per kilobase of transcript per million mapped reads (FPKM) to measure and normalize the expression levels of both mRNAs and lncRNAs. RSEM software is used to calculate the FPKM (Li and Dewey, 2011). The edgeR package (http://www.rproject.org/) was used to identify significant differentially expressed (DE) transcripts among the groups. We identified DE mRNA and DE lncRNA with absolute value of log 2 (Fold change) > 1 and a false discovery rate (FDR) < 0.05, and DE miRNA with absolute value of log 2 (Fold change) > 1 and p-value < 0.05. Gene Ontology (GO, https://geneontology.org) and Kyoto Encyclopedia of Genes and Genomes (KEGG, https://www. genome.jp/kegg/) analyses were carried out to have an insight to the biological function of the DE mRNAs. Both the GO and KEGG analyses were guided by Gene Denovo Biotechnology Co., Ltd (Guangzhou, China). GO terms and KEGG pathways with a p-value < 0.05 were considered significantly enriched.

Gene Set Enrichment Analysis
We conducted Gene Set Enrichment Analysis using GSEA software from the Broad Institute (http://www.broadinstitute. org/gsea/index.jsp) and MSigDB to identify whether a set of genes in specific GO terms\KEGG pathways\Disease Ontology (DO) terms showed a significant difference between the two groups. In brief, we input a gene expression matrix and classed genes by their signal-to-noise ratios. Enrichment scores and p-value were calculated with default parameters, p-value, p < 0.05; q-value, q < 0.25; and normalized enrichment score (| NES|) > 1.

LncRNA-mRNA Association Analysis
We performed antisense lncRNA, cis regulation and trans regulation analyses to predict the target genes of lncRNAs. Some lncRNAs may have capacities to set gene silencing, modulating transcription and impact mRNA stability. We recognize these lncRNAs as antisense lncRNAs. The relationships of antisense lncRNA and its target genes were predicted by RNAplex software (http://www.tbi.univie.ac.at/ RNA/RNAplex.1.html). Some lncRNAs can regulate their nearby genes in same allele, which situated less than 10 kilobases upstream or downstream of genes. We recognize these lncRNAs as cis lncRNAs. Besides, some lncRNAs can act on genes located distantly. Trans role is lncRNA's regulation for expression of trans chromosome genes. The correlation expression of lncRNAs and genes was conducted to identify the target genes of lncRNAs with p ≥ 0.9. Cytoscape software (v3.6.0) was used to establish and visualize the interaction network.

Competing Endogenous RNA Network Construction
Based on the ceRNA hypothesis, we constructed mRNA-miRNA-lncRNA networks. The negatively co-expressed mRNA-miRNA pairs and lncRNA-miRNA pairs were assessed by the Spearman rank correlation coefficient (SCC), with a SCC < -0.7. The coexpressed lncRNA-mRNA pairs were assessed using the Pearson correlation coefficient (PCC) with a PCC >0.9. The significance correlations between the common miRNA sponges and the two genes were calculated using the method of hypergeometric cumulative distribution. The mRNA-miRNA-lncRNA pairs were filtered with a p-value < 0.05. The interaction networks were visualized using Cytoscape software (v3.6.0).

Real-Time Quantitative PCR Validation
Sixteen transcripts-eight mRNAs and eight lncRNAs-were randomly selected to validate the RNA-Seq data by RT-qPCR. Primers as shown in Supplementary Table S1 were designed with Primer 5 software (Premier Biosoft International) and synthesized by Sangon biotech (Shanghai, China). RNA was reverse transcribed to cDNA using HiScript III RT SuperMix for qPCR (+gDNA wiper) (Vazyme Biotech, China). The RT-qPCR reaction system was mixed with cDNA template, primers, ChamQ SYBR Color qPCR Master Mix (Vazyme Biotech, China) and RNase-free water. All reactions were contained in 96 well plate and then conducted in a StepOne Plus Real-Time PCR system (Applied Biosystems, United States) with triplicate. Melting curves analysis was used to re-confrm the amplifcation of only a single PCR product. Using beta-2microglobulin (β-2-m) as a reference, the relative expression levels of mRNA and lncRNA were quantified with the 2 −ΔΔCt method (Livak and Schmittgen, 2001).

Overview of the RNA-Sequencing
In this study, 18 cDNA libraries were constructed for sequencing to conduct a mechanistic investigation of the brain and pituitary transcriptome mechanism and analyse the potential role of lncRNAs regulating the reproduction in tongue sole. A total of 734.6 million (734,595,534) and 720.1 million (720,070,248) clean Frontiers in Genetics | www.frontiersin.org December 2021 | Volume 12 | Article 802953 reads were generated from the brain and pituitary. After removing low-quality reads, an average of 98.66 and 98.54% of the high-quality reads from the brain and pituitary, respectively, were unmapped to the rRNA database. And the rRNA mapped reads were deleted. Among the unmapped reads, the mapping ratios and unique mapping ratio of tongue sole reference genome were 76.98 and 75.97%, respectively, in the brain and 74.36 and 72.22%, respectively, in the pituitary (Supplementary Table S2).

Functional Annotation and Classification of the Assembled Sequences
The novel transcripts were identified based on the location of the assembled transcripts in the reference genome, a transcript length of ≥200 bp and an exon number of ≥2. The functional annotation information of the novel transcripts is illustrated in Supplementary Table S3. The intersection of the CPC/CNCI/ SwissProt data was taken to screen for lncRNAs. A total of 3,897 lncRNA transcripts were obtained, including 1,941 intergenic lncRNAs (49.8%), 355 bidirectional lncRNAs (9.1%), 763 antisense lncRNAs (19.6%), 384 senseoverlapping lncRNAs (9.9%) and 455 other lncRNAs (11.7%). No intronic lncRNAs were identified ( Figure 1A). In addition, most novel lncRNA transcripts mainly contained between 1 and 4 exons and the number of lncRNA transcripts decreased with the increase of the number of exons ( Figure 1B).

Analysis of DE mRNAs and LncRNAs
In the brain, we identified 324 Figure 2B).
In the pituitary, we identified 603 DE mRNAs, 392 upregulated and 211 downregulated, in the ovarian stages IV vs. V comparison group and identified 5,027 DE mRNAs, 2,484 upregulated and 2,543 downregulated, in the ovarian stages V vs.    Figure 2D).
DE mRNAs and lncRNAs were illustrated in circos plots. The DE mRNAs and lncRNAs in the brain and pituitary in two comparison groups (stages IV vs. V and V vs. VI) are presented in Figure 3 and Supplementary Figure S1.

GO and KEGG Pathway Enrichment Analyses of DE mRNAs
We performed GO and KEGG enrichment analyses to identify the biological functions of the DE mRNAs. As shown in Supplementary Table S6, in the ovarian stages IV vs. V comparison, the most enriched GO terms in the brain were extracellular region part and NAD(P) + transhydrogenase activity, which were in the cellular component and molecular function categories, respectively. In the biological process category, the significantly enriched GO terms included cell migration, cell motility and cellular localization. In the ovarian stages V vs. VI comparison, the most enriched GO terms in the brain were cytosol and activin receptor signaling pathway, which were in the cellular component and biological process categories, respectively. As shown in Supplementary  Table S6, in the ovarian stages IV vs. V comparison, the significantly enriched GO terms in the pituitary included non-membrane-bounded organelle, cytoskeleton, oxidoreductase activity, growth factor binding, carbohydrate metabolic process and so on. In the ovarian stages V vs. VI comparison, the significantly enriched GO terms in the pituitary included cytoskeleton, transport vesicle membrane, nucleosidetriphosphatase activity, oxidoreductase activity and chromatin organization. Furthermore, several GO terms related to reproduction, including female pregnancy sex chromosome and fatty acid synthase activity, were enriched.
In the KEGG pathway analysis, we established a p-value of 0.05 as the threshold. In the brain, 8 and 23 pathways were significantly enriched in the ovarian stages IV vs. V and stages V vs. VI comparison groups, respectively. Among these pathways, the Jak-STAT signaling, calcium signaling and melanogenesis pathways were significantly enriched in both comparison groups ( Figure 4A). In the ovarian stages IV vs. V comparison, the cytokine-cytokine receptor interaction pathway was significantly enriched. In the ovarian stages V vs. VI comparison, several signaling pathways, including the GnRH signaling pathway, MAPK signaling pathway, ErbB signaling pathway, Wnt signaling pathway and insulin signaling pathway, were identified ( Figure 4B). In the pituitary, 7 and 18 pathways were significantly enriched in the ovarian stages IV vs. V and stages V vs. VI comparison groups, respectively. Among these pathways, the adipocytokine signaling pathway and cell cycle pathway were significantly enriched in both comparison groups.
In the ovarian stages IV vs. V comparison, pathways associated with reproduction and the cell cycle, such as oocyte meiosis, ECM-receptor interaction and p53 signaling pathway, were identified ( Figure 4C). In the ovarian stages V vs. VI comparison, several pathways related to metabolism, including fatty acid biosynthesis, degradation and metabolism, were identified ( Figure 4D).

GSEA
We carried out GSEA to identify more potential pathways that are not significantly different but have important biological significances. Generally, the significance of gene sets was established based on the following criteria: absolute normalized enrichment score (|NES|)>1, nominal (NOM) p-value < 0.05, and FDR q-value < 0.25 (Supplementary Table S7). In the brain, GSEA showed that 6 and 4 pathways were significantly enriched in the ovarian stages IV vs. V and stages V vs. VI comparisons, respectively. Among these pathways, some were related to genetic information processing, such as proteasome, protein export and ribosome biogenesis in eukaryotes. In addition, some pathways were related to immunity, such as intestinal immune network for IgA production and Toll-like receptor signaling pathway. In the pituitary, GSEA showed that 27 and 12 pathways were significantly enriched in the ovarian stages IV vs. V and V vs. VI comparisons, respectively. Some pathways, such as the p53 signaling pathway, cell cycle, glycosphingolipid biosynthesis-lacto and neolacto series, NOD-like receptor signaling pathway, DNA replication and Toll-like receptor signaling pathway, were significantly enriched in both comparison groups.

Identification and Enrichment Analysis of LncRNA Target Genes
To further predict the potential roles of the lncRNAs involved in modulating the reproductive process, we constructed interaction networks of lncRNAs and their antisense, cis and trans target genes ( Figure 5). In the brain, 57 lncRNAs had the potential to regulate 40 of the genes listed in Supplementary Table S8. Among these lncRNAs, 11 exhibited antisense complementary correlations with 12 DE mRNAs, 19 showed potential cis target regulatory relationships with 16 DE mRNAs, and 19 lncRNA-mRNA coexpression pairs were identified. In the pituitary, 42 lncRNAs had the potential to regulate 38 of the genes listed in Supplementary Table S8. Among these lncRNAs, 9 exhibited antisense complementary correlations with 12 DE mRNAs, 26 showed potential cis target regulatory relationships with 23 DE mRNAs, and 9 lncRNA-mRNA coexpression pairs were identified. The top 20 KEGG pathways of the three types of regulatory relationships are shown in Supplementary Figure S2 (also see Supplementary Table S9). The target genes of these DE lncRNAs were enriched in numerous KEGG pathways, including some signal transduction pathways, such as TGF-beta signaling pathway, Wnt signaling pathway and ErbB signaling pathway; some signaling molecule and interaction pathways, such as cytokine-cytokine receptor interaction, neuroactive ligandreceptor interaction, ECM-receptor interaction and cell adhesion molecules (CAMs); and some pathways related to cell growth and death, such as apoptosis, cell cycle, fatty acid biosynthesis and steroid hormone biosynthesis pathways. These pathways directly or indirectly influence reproductive processes, indicating that these lncRNAs may play a role in reproduction.

mRNA-MiRNA-LncRNA Interaction Analysis
We used MIREAP, miRanda and TargetScan to predict the targets of each miRNA after DE mRNAs, lncRNAs and miRNAs were recognized. CeRNA networks were built based on the mRNA-miRNA and lncRNA-miRNA pairs. In these ceRNA networks, each pair of mRNAs and lncRNAs were targeted and negatively co-expressed with a joint miRNA (Supplementary Table S10).
According to the result of ceRNA network and pathway enrichment analysis, we selected 4 reproduction-related genes in the brain-cga, lhb, cyp19a1 and collagen alpha-1 (IV) chain (col4a1) -to predict their interactions with lncRNAs. A total of 26 lncRNAs interrelated with these 4 mRNAs through competitive binding of 25 miRNAs. We also selected 5 reproduction-related genes in the pituitary-cyp19a1, mmp2, egfl6, smad4 and gprc5c-to predict their interactions with lncRNAs. A total of 63 lncRNAs interrelated with these 5 mRNAs through competitive binding of 10 miRNAs ( Figure 6).
GO and KEGG enrichment analyses were conducted to analyze the functions of the mRNAs in the ceRNA networks. GO terms and KEGG pathways with p-value < 0.05 were regarded as significant enrichment. In the brain, overall, 99 GO terms were identified. In the pituitary, overall, 132 GO terms were identified (Supplementary Table S11). In addition, brain mRNAs involved in ceRNA network were annotated to 9 pathways, including MAPK signaling pathway and endocytosis. In pituitary, mRNAs involved in ceRNA network were annotated to 13 pathways, including reproduction-related pathways such as MAPK signaling pathway, fatty acid elongation, oocyte meiosis and progesterone-mediated oocyte maturation (Supplementary Figure S3).

Validation of RNA-Seq Data by RT-qPCR
An equal number (8)  levels in the brain and pituitary from the ovarian stages IV vs. V and V vs. VI comparison groups were quantified by RT-qPCR. The expression patterns of selected transcripts from the RNA-Seq analysis were similar with RT-qPCR experimental validation results. These findings showed that our RNA-Seq results were usability and accuracy (Figure 7).

DISCUSSION
Reproduction, a complex biological process involving many organs, is precisely controlled by the HPG axis in vertebrates. The hypothalamus and pituitary can secrete a series of hormones triggering ovarian development, maturation and ovulation to maintain the reproductive cycle. However, more recent studies of the role of lncRNAs in reproduction in tongue sole have focused mainly on the ovary (Dong et al., 2021). In this study, the transcriptome analysis of the brain and pituitary of tongue sole at three ovarian developmental stages was performed to investigate the key gene functions and the relationships of mRNAs and lncRNAs. We conducted KEGG pathway analysis to identify the biological functions of the DE genes. In this study, some pathways were significantly enriched in the brain at ovarian stages V vs. VI, for example, GnRH signaling pathway, calcium signaling pathway and MAPK signaling pathway. Previous studies have shown that the principal molecule controlling the HPG axis in mammals is the decapeptide GnRH, which is secreted by the hypothalamus and is a major regulator of reproduction. GnRH acts on the pituitary where its receptors located to control the synthesis and release of the gonadotropins,  LH and FSH (Blanco, 2020). These hormones are glycoprotein heterodimers, which are composed of a common α submit (Cga, or αGSU) and unique β subunits (Lhβ and Fshβ) (Shi et al., 2015;Coss, 2018).
As shown in Figure 8, the GnRHR, which belongs to the 7 transmembrane domain-containing class A superfamily of G-protein-coupled receptors (GPCRs), binds to G q/11 proteins to activate phospholipase C (PLC), which conveys its signal to diacylglycerol (DAG) and inositol 1,4,5-trisphosphate (IP3) (Chang and Pemberton, 2018;Janjic et al., 2019). In this study, we found that members of the GPCR family, including GPCR family C group 5 member C (Gprc5c) GPR52, and GPR4, were significantly differentially expressed in the pituitary at ovarian stages V vs. VI. All of these genes were expressed mainly in the brain and pituitary. DAG and IP3 activate the intracellular protein kinase C (PKC) pathway and stimulate the release of intracellular calcium, respectively. The activation of PKC transactive the epidermal growth factor (EGF) receptor and active mitogen-activated protein kinases (MAPKs). Active MAPKs translocation to the nucleus leads to activation of transcription factors and rapid induction of early genes, including lhβ, fshβ and cga. In the present study, the cga was highly expressed in the pituitary compared to the ovary and brain. In addition, cga was significantly differentially expressed in the brain in the ovarian stages IV vs. V comparison and exhibited the highest expression in the brain during ovarian stage IV. Similar to cga, lhβ was also expressed mainly in the pituitary and brain, especially in the pituitary. The lhβ was significantly differentially expressed in the brain in the ovarian stages IV vs. V comparison and exhibited the highest expression in the brain during ovarian stage IV.
In the present study, we found that the cytokine-cytokine receptor interaction pathway was significantly enriched in the brain at ovarian stages IV vs. V. Cytokines are soluble extracellular proteins or glycoproteins that are crucial intercellular regulators of cells involved in innate immunity, cell growth, cell differentiation, and cell death. Cytokines are released by diverse cells in the body, and they stimulate responses by binding to specific receptors on the surface of target cells. As shown in the GnRH signaling pathway (Figure 8), EGF, as a kind of cytokine, binds to EGFR on the cell surface to activate downstream pathways. In vertebrates, hypothalamic peptides affect the production of pituitary hormones, which are regulated by several growth factors, including EGF, one of the most extensively studied growth factors (Smith et al., 2015). In mammals, EGF can induce prolactin (PRL) and LH secretion, suggesting that EGF might play a crucial role in growth and reproduction (Przylipiak et al., 1988). A previous study in zebrafish (Danio rerio) showed that EGF did not affect lhβ expression in pituitary cells (Lin and Ge, 2009). However, EGF was found to significantly reduce lhβ expression in grass carp (Ctenopharyngodon idellus) (Hu et al., 2020). In our study, epidermal growth factor-like protein 6 isoform X2 (egfl6) was found to be significantly differentially expressed at ovarian stages V vs. VI and exhibited the highest expression during ovarian stage IV. Moreover, lhβ exhibited the highest expression during ovarian stage IV in the pituitary, with the same expression trend as egfl6, indicating that EGF may also affect pituitary hormones in tongue sole. Furthermore, we found that epidermal growth factor receptor isoform X2 (egfr) was significantly differentially expressed in the brain at ovarian stages V vs. VI. Metabolism has a strictly connection with reproduction, and many studies have focused on elucidating the mechanism by which signals from the periphery are transmitted to the HPG axis in multiple metabolic states (Sliwowska et al., 2014). Maintaining the balance between reproduction and energy metabolism is important. The process of reproduction requires energy, which is stored mainly as fat, glycogen and glucose. Insulin has been recognized as a key hormone for metabolic regulation in vertebrates (Irwin, 2019). The major function of insulin is to maintain peripheral glucose homeostasis through stimulation of glucose uptake, oxidation and storage (Deck et al., 2017). However, insulin also plays an important role in regulating reproduction and may play a vital role in the interaction between metabolism and reproduction in mammals (Sliwowska et al., 2014). Studies have shown that twice daily insulin treatment can reverse the 50% reduction in LH in diabetic male rats (Rattus norvegicus) compared to nondiabetic controls (Dong et al., 1991). However, research showed that in a growthrestricted, hypogonadotropic lamb model, central injection of insulin did not increase LH secretion (Hileman et al., 1993). Therefore, the role of insulin as a regulator in the HPG axis in mammals remains controversial. In present study, we found that the insulin signaling pathway was significantly enriched in the brain at ovarian stages V vs. VI, indicating that insulin may impact reproduction in tongue sole. In the silver pomfret (Pampus argenteus), gonadal development, especially vitellogenesis and spermatogonia proliferation, is closely related to IGFs (Gu et al., 2021). IGF has also been proven to control reproduction in the sapphire devil, Chrysiptera cyanea (Mahardini et al., 2018). In our study, insulin-like growth factor 2 (igf2) was differentially expressed in the pituitary at ovarian stages V vs. VI. Other metabolic pathways were enriched in the pituitary at ovarian stages V vs. VI, for example, fatty acid metabolism, biosynthesis and degradation and galactose metabolism. Fatty acids have already been proven to play an important role in reproduction (Sorbera et al., 2001).
The extracellular matrix (ECM) is a noncellular component of tissues and organs that maintains cell and tissue structure and functions. The ECM exerts direct or indirect control over cellular activities such as cell adhesion, migration, differentiation, proliferation and apoptosis (Walker et al., 2018). The ovary has been proven to undergo extensive tissue remodeling throughout each reproductive cycle (Curry Jr and Osteen, 2003). The matrix metalloproteinase (MMP) family is the important enzyme family accountable for ECM remodeling. The MMP system is generally believed to participate in follicle rupture during ovulation in medaka (Oryzias latipes) (Ogiwara et al., 2005;Takahashi et al., 2013). Mmp2 has been observed to possibly be involved in regulating the activity of stem/progenitor cells that give rise to new granule neurons (Sîrbulescu et al., 2015). Studies have shown that the activation of extracellular signalregulated kinases 1 and 2 (ERK1/2), which is included in response to GnRHR, occurs through the sequential activation of PKC, MMPs and EGF receptor transactivation (Shah and Catt, 2004). Further studies showed that MMP2 and MMP9 may participate in the endocrine regulation of pituitary gonadotrophs, suggesting a vital role for MMPs in GnRH signaling (Roelle et al., 2003). In addition, MMP2 induces the hydrolysis of type IV collagen, which constitutes the basement membrane, and membranetype 2-MMP degrades type I collagen present in the theca cell layer in medaka (Ogiwara et al., 2005). Research has indicated that the col4a1 knockdown decreases cell viability and suppresses cell cycle progression in breast cancer cells (Salem et al., 2016). Moreover, the type IV collagen has been detected in the basement membranes of the rat anterior pituitary (Vila-Porcile et al., 1992). Collagen, as a component of the ECM, participates in the tissue remodeling. We found that mmp2 was differentially expressed and exhibited the highest expression at ovarian stage VI in the pituitary. We also found that col4a1 was differentially expressed and exhibited the highest expression at ovarian stage V in the pituitary. In ovarian stages V vs. VI, mmp2 exhibited increased expression, while col4a1 exhibited decreased expression. To date, the detailed function of mmp2 and col4a1 in the reproduction of teleosts is poorly characterized. Considering these observations, it is tempting to speculate that mmp2 and col4a1 in the pituitary are associated with neuroendocrine regulation by a novel mechanism in the reproduction of tongue sole.
Recent studies have identified the roles of lncRNAs in reproduction, focusing mainly on the ovary. However, the expression patterns and potential physiological funtions of lncRNAs in the brain and pituitary are unclear. Therefore, we used transcriptome sequencing to explore lncRNA expression profiles in the brain and pituitary at three ovarian stages. A total of 22,549 lncRNAs (9,816 known and 12,733 novel) were identified. The present study indicated that brain and pituitary DE lncRNAs and their target DE genes might have regulatory roles in reproduction in tongue sole. Thus, we constructed lncRNA-mRNA interaction networks based on their antisense and cis regulatory mechanisms or coexpression relationships. The networks indicate that smad4 is the cis target gene of lncRNA XR_522655.2. The porcine SMAD4 protein has high expression in granulosa cells (GCs) and oocytes from primary, preantral, and antral follicles and lower expression in apoptotic ovarian GCs, indicating that SMAD4 may participate in ovarian development and selection . SMAD4 is a key component of the TGF-β signaling pathway, which participated in the various processes in the mammalian ovary and participates in feedback regulation between the ovary and pituitary Du et al., 2020). Smad4 dysregulation is associated with disorders and delays in nervous system development . In goldfish (Carassius auratus), smad4 has been proven to have a high expression in the pituitary and might affect FSH synthesis (Lau and Ge, 2005). In addition, accumulating evidence shows that several lncRNAs may have a potential role in regulating smad4, for example, lncRNA PCAT7 in prostate cancer cells (Lang et al., 2020), lncRNA 9884 in renal inflammatory cells (Zhang et al., 2019), LINP1 in lung cancer cells (Zhang et al., 2018) and LIN-LET in bladder cancer cells (Zhuang et al., 2017). In our study, we also found that lncRNA XR_522655.2 might have a cis regulatory relationship with smad4 in the brain and pituitary of tongue sole for reproduction.
In the lncRNA-mRNA networks, the lncRNA TCONS_00069197 is antisense-acting toward its target gene prl, a gene encoding PRL. PRL is an anterior pituitary Frontiers in Genetics | www.frontiersin.org December 2021 | Volume 12 | Article 802953 hormone with one of the broadest ranges of function among vertebrate hormones (Whittington and Wilson, 2013). PRL is a cytokine based on its structure and receptor type, which possesses 1 transmembrane domain, together with growth hormone (Dobolyi et al., 2020). In mammals, PRL has been proven to take part in reproductive development and cycling, including oocyte meiosis and maturation (Sorokowski et al., 2019). PRL also plays an important role in reproduction of nonmammalian species. In chickens, PRL participates in the mechanisms orchestrating ECM turnover during ovarian follicular development in the hen (Gallus domesticus) ovary by regulating the transcription, translation, and/or activity of some components of the MMP system (Hrabia et al., 2021). Similarly, in fish, PRL has been shown to play a role in reproductive development and sexual maturity. The prl has a significantly lower expression level in juvenile male blue gouramis (Trichogaster trichopterus), than in reproductively mature males, suggesting that PRL may involve in endocrine control of gonadal development and reproduction (Degani et al., 2010). Furthermore, prl has been detected in all the brain areas of orange-spotted grouper (Epinephelus coioides) (Zhang et al., 2004), yellow perch (Perca flavescens) (Lynn and Shepherd, 2007), goldfish (Imaoka et al., 2000), and rainbow trout (Oncorhynchus mykiss) (Yang et al., 1999). In this study, prl displayed high expression in the pituitary and brain, especially in stage IV. Similar to prl, the lncRNA TCONS_00069197 also had relatively high expression in the pituitary and brain and especially in stage IV, indicating that this lncRNA may play a role in regulating the nero-and endocrine reproductive axes, such as steroidogenesis and gonadogenesis, in tongue sole. We also constructed mRNA-miRNA-lncRNA networks based on the ceRNA hypothesis. Interestingly, in the brain, the genes closely related to reproductive processes, including cga, lhb, cyp19a1 and col4a1, were potentially modulated by 26 lncRNAs via competition for 25 miRNAs. In the pituitary, the genes closely associated with reproductive processes, including cyp19a1, mmp2, egfl6, smad4 and gprc5c, were potentially modulated by 63 lncRNAs via competition for 10 miRNAs. Aromatase is encoded by the cyp19a1s gene. It is a key enzyme that converts androgens into estrogens and plays a crucial role in regulating numerous physiological functions in vertebrates, including sex differentiation, gonadal development and feedback effects of steroid hormones on gonadotropin secretion (Lin et al., 2020). In teleosts, there are two types of cyp19a1s, encoding aromatase a and b, which are expressed mainly in the gonads and brain, respectively. Teleost cyp19a1a expressed in the gonad is involved in gonadal sex differentiation to trigger and maintain ovarian differentiation (Guiguen et al., 2010). The expression of teleost cyp19a1b in the brain is connected with brain cell proliferation, neurogenesis and brain development and repair (Ulhaq and Kishida, 2018). In this study, cyp19a1 was expressed mainly in the brain and pituitary. Sixteen lncRNAs were shown to modulate cyp19a1 via 8 miRNAs in the brain, and 4 lncRNAs were shown to modulate cyp19a1 via 1 miRNA in the pituitary. Moreover, the KEGG pathways progesterone-mediated oocyte maturation, oocyte meiosis and MAPK signaling pathway were enriched in DE mRNAs in the ceRNA network. These results indicate that lncRNAs may be regulators of gene expression in the brain and pituitary via miRNAs, which could play a potential role in tongue sole reproduction regulated by the HPG axis.

CONCLUSION
This study revealed the key lncRNAs and mRNAs involved in the HPG axis regulating the reproduction of tongue sole. The enrichment analysis indicated that the target genes of the DE lncRNAs were enriched in a series of pathways related to reproduction. We constructed lncRNA-mRNA regulatory networks based on antisense, cis and trans regulatory mechanisms. In addition, based on the ceRNA hypothesis, we constructed mRNA-miRNA-lncRNA (i.e., ceRNA) networks to predict the target genes of lncRNAs via miRNAs. In summary, this study provides a valuable resource for lncRNAs modulating the expression of mRNAs that play a crucial role in tongue sole reproduction and offers deeper insight into the molecular mechanisms fundamental to improving the efficiency of tongue sole reproduction.

ETHICS STATEMENT
The animal study was reviewed and approved by the Animal Care and Use Committee at the Chinese Academy of Fishery Sciences.

AUTHOR CONTRIBUTIONS
BS conceived the project and designed the experiments. YD and LL analyzed the data. YD performed the experiments and wrote the first draft of manuscript. BS, HW, and YD critically revised the manuscript. All authors contributed to the article and approved the submitted version.