Effects of gonadotropin-releasing hormone analog (GnRHa) immunization on the gonadal transcriptome and proteome of tilapia (Oreochromis niloticus)

Background: Gonadotropin releasing hormone (GnRH) plays an important role in the regulation of vertebrate reproduction. Studies have shown that immunization against GnRHa can induce sexually sterile tilapia. To explore the mechanism behind this, in this study, RNA-seq and data-independent acquisition (DIA) techniques were used to study the transcriptome and proteome of the gonad of tilapia immunized with GnRHa. Results : 644 differentially expressed genes (80 upregulated and 564 downregulated) and 1150 differentially expressed proteins (351 upregulated and 799 downregulated) were identied. There were 209 genes with consistent differential expression patterns in the transcriptomic and proteomic analyses, of which 9 were upregulated and 200 downregulated, indicating that the gonad gene expression was inhibited by GnRHa immunization. The downregulated genes were particularly involved in the functions of single-organism process, binding, cellular process, metabolic process and catalytic activity, and associated with the pathways including ECM–receptor interaction, focal adhesion, cardiac muscle contraction and oxidative phosphorylation. The expression of six differentially expressed genes involved in the GnRH signaling pathway was all downregulated. In addition, several important functional genes related to gonadal development after GnRHa immunization were screened. Conclusions: This study conrmed the expression of corresponding genes was affected by GnRHa on the gonad development in tilapia at the molecular level, and laid a foundation for elucidating the mechanism of GnRHa immunization. hormone; LH: luteinizing hormone; FDR: false discovery rate; FC: fold change; T1: immunized group; CK: control group; DEGs: differentially expressed genes; DEPs: differentially expressed proteins; DDA: data-dependent acquisition; GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; qRT-PCR: quantitative real-time reverse transcription-PCR; FPKM: fragments per kilobase of transcript per million


Abstract
Background: Gonadotropin releasing hormone (GnRH) plays an important role in the regulation of vertebrate reproduction. Studies have shown that immunization against GnRHa can induce sexually sterile tilapia. To explore the mechanism behind this, in this study, RNA-seq and data-independent acquisition (DIA) techniques were used to study the transcriptome and proteome of the gonad of tilapia immunized with GnRHa.
Results : 644 differentially expressed genes (80 upregulated and 564 downregulated) and 1150 differentially expressed proteins (351 upregulated and 799 downregulated) were identified. There were 209 genes with consistent differential expression patterns in the transcriptomic and proteomic analyses, of which 9 were upregulated and 200 downregulated, indicating that the gonad gene expression was inhibited by GnRHa immunization. The downregulated genes were particularly involved in the functions of single-organism process, binding, cellular process, metabolic process and catalytic activity, and associated with the pathways including ECM-receptor interaction, focal adhesion, cardiac muscle contraction and oxidative phosphorylation. The expression of six differentially expressed genes involved in the GnRH signaling pathway was all downregulated. In addition, several important functional genes related to gonadal development after GnRHa immunization were screened.
Conclusions: This study confirmed the expression of corresponding genes was affected by GnRHa on the gonad development in tilapia at the molecular level, and laid a foundation for elucidating the mechanism of GnRHa immunization.

Background
Gonadotropin releasing hormone (GnRH) is a neurohormone secreted by hypothalamic neurons, which plays important roles in the regulation of vertebrate reproduction. It is transported to the portal vein system through axons, flows into the anterior pituitary through blood vessels, binds to the GnRH receptor (GnRHR) on gonadotropin cells in the pituitary, and stimulates the synthesis and secretion of follicle-stimulating hormone (FSH) and luteinizing hormone (LH) [1]. GnRH is a small-molecule polypeptide with weak immunogenicity, however, coupling with macromolecular protein carriers can improve its immunogenicity and induce immunized animals to produce GnRH antibodies [2]. GnRH antibody can neutralize endogenous GnRH, leading to infertility in male and female mammals [3]. The efficiency of GnRH purification is low due to the low content of endogenous GnRH, and synthetic gonadotropin releasing hormone analog (GnRHa) is often used as a clinical alternative. GnRHa not only retains the physiological functions of GnRH, but also exhibits a prolonged half-life and improved capacity to bind to GnRH receptors [4]. GnRHa antigen with good immunogenicity can be prepared by coupling GnRHa with bovine serum albumin (BSA) and emulsifying with Freund's incomplete adjuvant [5]. GnRH immunization is considered to be an effective method to improve the growth performance of cultured animals, improve meat quality, control reproduction, and prevent undesired sexual behavior [6,7]. The reproduction of fish is regulated by the release of gonadotropin from the pituitary induced by GnRH released from the hypothalamus [8]. Studies have shown that fish have two or three forms of GnRH [9], of which one or two forms can be found in the pituitary [10]. However, few reports have been published on the regulatory effect of GnRH immunization on the reproductive function of fish.
Tilapia is one of the most important freshwater culture varieties, exhibiting the advantages of rapid growth and strong stress resistance, which is an international fish that the Food and Agriculture Organization of the United Nations focuses on promoting aquaculture [11]. However, owing to its early sexual maturity and rapid breeding cycle, the reproductive behavior of tilapia in mixed-sex culture consumes a lot of energy, which interferes with growth, and at the same time leads to very different qualities of harvested tilapia, a low commodity rate, and a high feed conversion ratio, greatly reducing the product quality and associated economic benefits [12]. As such, it is important to obtain sterile individuals to eliminate the adverse effects of sexual maturity and reproductive behavior of tilapia. Our previous studies showed that the level of GnRH antibody in the GnRHa immunized group was significantly higher than that in the negative control group [13]. To further investigate the mechanism behind the effects of GnRHa immunization in tilapia, the transcriptome and proteome of gonad of tilapia after GnRHa immunization were analyzed, and the genes with functions that are particularly important for gonad development of tilapia were screened. The results confirmed the expression of corresponding genes was affected by GnRHa on gonad development at the molecular level, and laid a foundation for elucidating the mechanism of GnRHa immunization.

RNA-seq
The quality control results of the original data of RNA-seq showed that a total of 375.2 million clean reads were obtained from the eight samples after preliminary filtration, and the average number of clean reads per sample was 46.9 million (ranging from 35.9 to 59.5 million). Overall, 98.45-98.98% high-quality clean reads were obtained after further filtration of the clean reads, and the filtered reads included 0.31-0.36% reads containing an adapter, 0.71-1.21% low-quality reads (containing more than 50% bases with Qvalue ≤ 20), and 53-148 reads containing all A, no reads containing more than 10% N. After quality control, the data were considered to be of sufficient quality to be used for subsequent analysis. Owing to the sample quality and species, the efficiency of removing ribosomal RNA might not be stable, and the contamination of ribosomal RNA would affect the subsequent analysis. Therefore, the read alignment tool was applied to align high-quality clean reads to the ribosome of the species (mismatch number: 0), and the aligned reads were removed; the remaining data were used for subsequent analysis. The alignment results showed that the proportion of unmapped reads obtained from the samples was 98. 35

Identification of DEGs
In this study, edgeR software was used to analyze the difference of gene expression between the groups, and FDR and log 2 FC were used to screen the DEGs. A total of 644 DEGs were obtained. As shown in Fig. 1, 80 upregulated genes and 564 downregulated genes were identified in the immunized group (T1), compared with the control group (CK).

GO and KEGG analyses of DEGs
The DEGs were annotated by GO, and then the upregulated and downregulated genes was conducted by the GO term classification. The results (Fig. 2) showed that the upregulated genes were mainly involved in catalytic activity (13 genes) and metabolic process (11 genes), while the downregulated genes were mainly related to single-organism process (116 genes), binding (113 genes), cellular process (103 genes), metabolic process (83 genes), and catalytic activity (73 genes).
Functional enrichment analysis of the DEGs revealed the 20 most enriched pathways, as shown in

Identification of DEPs
Quality control of the original mass spectrometry data was carried out using QuiC software (Biognosys), and the similarity of quality control indexes between samples was determined. Then, a library was established based on the data obtained by DDA mode using Pulsar software, and was used for the quantification of subsequent DIA data. The identification standard of proteins was as follows:

GO and KEGG analyses of DEPs
The DEPs were annotated to GO terms, which were then classified according to the up-or downregulation (Fig. 5). The results showed that the DEPs were mainly enriched in metabolic process

Correlation analysis of transcriptome and proteome in tilapia
The correlations of the transcriptome and proteome of tilapia in the immunized group and the control group were analyzed, and functional enrichment analysis of the related gene sets was carried out.
The results are shown in Fig. 7. The findings for a total of 209 mRNAs were consistent with the differential expression pattern of the corresponding proteins, of which 9 were upregulated and 200 were downregulated. Significance analysis showed that the findings of 48 mRNAs were consistent with the differential expression patterns of the corresponding proteins, of which 1 was upregulated and 47 were downregulated. The detailed information of the genes with a consistent pattern are shown in Table 3. Table 3 Genes with a consistent pattern of differential expression between transcriptome and proteome To further validate the results from the RNA sequencing data, six genes were selected for qRT-PCR analysis. Of the six selected genes, all showed similar expression patterns in the qRT-PCR analysis as observed from RNA-seq data (Fig. 8). The statistical analysis also showed very good correlation (r = 0.997) between the two types of analysis.

Discussion
In addition to genetic factors, external environmental factors also play important roles in the sex determination and differentiation of fish. Therefore, fish can be useful subjects to study sexual differentiation. With the continuous progress and development of transcriptomic and proteomic technologies, transcriptomics and proteomics have been widely used in the study of reproductive and sexual differentiation-related mechanisms in fish, including tilapia, especially in the identification of DEGs, the screening of functional genes, and the discovery and verification of important regulatory pathways [14][15][16][17]. After treatment with progesterone receptor inhibitor (RU486), 7148 genes were differentially expressed in gonad of female tilapia; the results revealed that fshr and lhr were significantly downregulated and ars was significantly upregulated after RU486 treatment, which might account for the masculinization and infertility of female fish [18]. Transcriptomic analysis of the gonads of female and male tilapia at different developmental stages revealed that estrogen may play an important role in female sex determination and maintenance of phenotypic sex, which lays the foundation for future studies into the molecular mechanisms of sex determination and maintenance of phenotypic sex in non-model teleosts [19]. Analysis of the transcriptome of the gonads of control female, high-temperature-treated female, and high-temperature-induced neomale tilapia identified a number of genes that may be involved in GSD + TE (genotypic sex determination + temperature effects), which should be useful for investigating the molecular mechanisms of GSD + TE in fish [20].
Study of the gonadal proteome in fish during sex reversal or gonadal differentiation, and the screening of important functional proteins related to reverse acquisition and gonadal differentiation are of great significance for the study of vertebrate sex determination and differentiation [21,22].
GnRH immunization of mammals, such as mice, sheep, cattle, and pigs, can reduce the levels of FSH and LH, and inhibit the development of the gonads, achieving artificial intervention in animal sex differentiation [12]. Intervention in fish sex differentiation by GnRH immunization, investigation of the mechanism of action of GnRH in regulating sex differentiation, and screening of sex determination factors are of great significance for uncovering the sex-differentiation and sex-controlling mechanisms. However, no transcriptomic or proteomic studies of tilapia gonad after GnRH immunization have been reported.
In this study, the gonadal transcripts of tilapia immunized with GnRHa were sequenced, and an average of 88.28% clean reads were mapped to the reference genome. Sun   Forty-seven downregulated genes were obtained by combined analysis of the transcriptomic and proteomic data, among which the differential expression of zona pellucida sperm-binding protein 3 (ZPBP3), apolipoprotein A-I (apo A-I), and cytoplasmic 1 was most significant. Zona pellucida sperm binding protein (ZPBP) is closely related to sperm motility, capacitation, acrosome reaction, and sperm-egg binding. Studies have shown that ZPBP1 and ZPBP2 were mainly expressed in testicular tissue of mammals, while ZPBP3 was found in tilapia [23,24]. In this study, it was found that ZPBP3 could be expressed in the gonad of female tilapia, which may be related to the specific sex differentiation mechanism of fish. Meanwhile, the expression of ZPBP3 in gonad of tilapia was significantly downregulated after GnRHa immunization, indicating that ZPBP3 plays an important role in the gonad development of tilapia. Apo A-I is the major component of high-density lipoprotein, which plays an important role in reverse cholesterol transport and lipid metabolism, and is also a very important innate immune molecule, providing a platform for the assembly of several immune complexes [25,26]. Previous studies showed that apo A-I was present in multiple tissues of Branchiostoma belcheri, and the expression of apo A-I mRNA was highest in the gonads; the expression of apo A-I mRNA increased significantly in the process of infection, suggesting that apo A-I may be involved in immune stress response and play an important role in the innate defense immune system of Amphioxus [27].

The apo A-I expression profile in Monopterus albus immunized with
Aeromonas hydrophila showed that infection could downregulate the expression of apo A-I in the small intestine, spleen, and liver, indicating that apo A-I might be involved in the innate immune system of this species [28]. In this study, the expression of apo A-I in gonad of tilapia immunized with GnRHa was significantly downregulated, which also indicated that apo A-I plays an important role in the immune response of tilapia. Actin is one of the most highly conserved proteins in eukaryotic cells.
It is mainly present in the cytoplasm, is the main component of cytoskeleton microfilaments, and is necessary for a variety of cell functions, such as the division, movement, and growth of cells [29].
Actin (cytoplasmic 1) can be expressed in several tissues of tilapia, but its expression level in gonad is much higher than that in other tissues, indicating that this gene may play an important role in gonad development and gonad functions of tilapia [24]. In this study, the expression of actin (cytoplasmic 1) in gonad of tilapia immunized with GnRHa was significantly downregulated, which indicated that GnRHa immunization could inhibit the expression of important genes involved in gonad development and provide a reference for revealing the mechanism of GnRH immunization.

Conclusions
To elucidate the mechanism behind the immunosuppressive effect of GnRHa on gonad development in tilapia, this study analyzed the transcriptome and proteome of tilapia gonad after GnRHa immunization. The results showed that, compared with that in the control group, the gonad gene expression after GnRHa immunization was significantly different, and the number of downregulated genes was significantly higher than that of the upregulated genes. Six differentially expressed genes involved in the GnRH signaling pathway were screened, all of which were downregulated, indicating that GnRHa immunization inhibited the expression of genes related to the gonadotropin releasing hormone signal transduction pathway in tilapia. In addition, several important functional genes related to GnRHa immunity were screened, which laid a foundation for revealing the mechanism behind the effect of GnRH on gonad development in tilapia. Further studies on the expression patterns of screened GnRHa immune-related genes and the verification of gene functions are currently underway.

Fish and GnRHa immunization
All female tilapia (Oreochromis niloticus) with average weight of 30.22 ± 1.80 g were provided by the National Tilapia Seed Farm (Nanning, Guangxi, China). According to the amino acid sequence of tilapia GnRH in Uniprot database (entry name: Q76FQ2), we had designed and synthesized tilapia GnRHa (Pyr-His-Trp-Ser-Tyr--Leu-Arg-Pro-NHEt). The GnRHa antigen was prepared according to the method we reported before [13]. The tilapia were randomly divided into two groups. One group (the GnRHa-immunized group; T1) was immunized by intraperitoneal injection with 0.5 mL/fish of GnRHa Four replications were included in each group with 60 tilapia per replication, which were named T1-1, T1-2, T1-3, and T1-4 and CK-1, CK-2, CK-3, and CK-4. During the study, the water temperature of tilapia culture was kept at 27.5 ± 1.5 °C, and the water was aerated continuously for 24 h. The water was refreshed every 2 days and the fish were provided with feed equivalent to 3% of their body weight daily. The basic fish feed comprising soybean meal, fish meal, and flour was used as feed (purchased from Nanning Tongwei Feed Company; floating fish feed).

Sampling and gender identification
Feeding was stopped 24 h prior to sampling. The sex of tilapia was judged according to the external genitalia. The fishes were anesthetized with MS222. Then, the complete gonad was isolated and washed with PBS. For each repetition, the gonads of three tilapia were mixed, put into an RNase-free cryotube, and stored at − 80 °C after freezing in liquid nitrogen.
RNA extraction, RNA-seq library construction, and Illumina deep sequencing The total RNA of tissues was extracted using Trizol reagent (Invitrogen), in accordance with the manufacturer's instructions. Then, genomic DNA was removed from the RNA sample using DNase.
Two microliters of total RNA was separated by 1% agarose gel electrophoresis to preliminarily determine the integrity and purity of RNA samples. Then, the concentration, quality, and integrity of RNA were determined using a NanoDrop 2000 spectrophotometer (Thermo Scientific) and an Agilent

Bioinformatic analysis of RNA-seq data
To acquire clean reads, the raw reads were filtered by removing reads containing adapters, reads containing more than 10% unknown nucleotides (N), and low-quality sequences [reads containing more than 50% low-quality (Qvalue ≤ 10) bases]. To avoid interference of ribosomal RNA, we mapped the clean reads to the ribosome database using bowtie [30] and removed the matched reads from the ribosome database. We mapped the remaining reads to the genome sequence of tilapia (GCF_001858045.2) using TopHat2 [31]. The sequencing data were deposited in the NCBI Sequence Read Archive under accession numbers SRR9937067-SRR9937074. The aligned reads were assembled into transcripts using cufflinks, and the assembled transcripts were merged using cuffmerge [32]. The abundance of gene transcripts was calculated via FPKM (fragments per kilobase of transcript per million mapped reads) [33], and genes with mean abundance of > 0 FPKM in any one of these samples were regarded as being expressed.
The raw read counts were used for differential expression analysis using edgeR [34]. Genes with fold change of FDR < 0.05, |log 2 FC| > 1, and p < 0.05 were considered to be differentially expressed. In this study, the differentially expressed genes (DEGs) between the immunized group (T1) and the control group (CK) were analyzed. After obtaining the DEGs, Gene Ontology (GO) classification was performed on them using the GOseq R packages based on a Wallenius noncentral hypergeometric distribution [35], which can adjust for gene length bias in DEGs. The statistically significant enrichment of DEGs in KEGG pathways was tested using KOBAS software [36].
Protein extraction, protein digestion, and separation high pH reverse phase  Mass spectrometric data analysis Raw data of DDA were processed and analyzed by Spectronaut Pulsar 11.0 (Biognosys AG). Pulsar was set up to search the database of O. niloticus assuming trypsin as the digestion enzyme.
Carbamidomethylation was specified as the fixed modification. Oxidation was specified as the variable modification. Raw data of DIA were processed and analyzed by Spectronaut Pulsar 11.0 (Biognosys AG) with the default parameters. After Student's t-test, differentially expressed proteins (DEPs) were filtered if their Qvalue was < 0.05 and absolute AVG log 2 ratio was > 0.58.

Protein functional annotation and enrichment analysis
Proteins were annotated against the GO, KEGG, and COG/KOG databases to obtain their functions.
Significant GO functions and pathways were examined among the DEPs with Qvalue ≤ 0.05.

Validation of DEGs by qRT-PCR
To validate the DEGs obtained by RNA sequencing, qRT-PCR was carried out with the total RNA used for RNA sequencing to examine the expression of 6 DEGs (2 upregulated unigenes and 4 downregulated unigenes). In the qRT-PCR analysis, the β-actin gene was amplified as an endogenous control. All samples were tested in triplicate. Specific primers were designed by the Primer Express 3.0 Software (Applied Biosystems), and are listed in Table 4. The relative changes of gene expression were calculated by 2 −ΔΔCt methods [38]. The value of log 2-fold change was used for graphing. Table 4 Primer sequences used for quantitative real-time reverse transcription-PCR. Transcript name  Forward Primer  Reverse Primer  ncbi_109194882  uncharacterized  GCCATCACTGTTCCCATCC  A  CACCCTCCGTAACCTCGT  CT  ncbi_109195722  uncharacterized  GCCATCACTGTTCCCATCC  A  CACCCTCCGTAACCTCGT  CT  ncbi_100695709  EGFR  ATGCTGTGGACGCTGATG  AA  GATGCTGCTGTTGAGGCT  TG  ncbi_100695964  MT1  GGATAAGTGGTTCTGGCG  GG  GGTGTAGTAGAGAGCAGC  GT  ncbi_100700606  MMP2  GCTACCCGAAAAGGTTGT  CC  TGATTTCACCCAACTTCAC  GAT  ncbi_100703685 c The genes and proteins differentially expressed between the immunized group and the control group.   The location of the six genes in the GnRH signaling pathway.   Correlation analysis of the protein and mRNA expression between the immunized group and the control group. The x-axis represents the log2 value of the fold differential gene expression, the y-axis represents the log2 value of the fold differential protein expression, and the dotted line indicates the threshold of significant differential expression.

Supplementary Files
This is a list of supplementary files associated with this preprint. Click to download. NC3RsARRIVEGuidelinesChecklist2014.pdf