Comparative transcriptome analysis reveals ectopic delta-5 and delta-6 desaturases enhance protective gene expression upon Vibrio vulnificus challenge in Tilapia (Oreochromis niloticus)

Tilapia (Oreochromis niloticus) cultures are frequently infected by Vibrio vulnificus, causing major economic losses to production units. Previously, tilapia expressing recombinant delta-5 desaturase and delta-6 desaturase (D56) were found to be resistant to V. vulnificus infection. In this report, we profile the D56-mediated molecular changes underlying this resistance in tilapia. A comparative transcriptome analysis was performed on V. vulnificus-infected wild-type and D56-transgenic tilapia using Illumina’s sequencing-by-synthesis approach. Gene enrichment analysis on differentially expressed unigenes was performed, and the expression patterns were validated by real-time PCR. Comparative transcriptome analysis was performed on RNA-sequence profiles obtained from wild-type and D56-transgenic tilapia at 0, 6 and 24 h post-infection with V. vulnificaus. GO and KEGG gene enrichment analyses showed that D56 regulates several pathways and genes, including fatty acid (FA) metabolism associated, and inflammatory and immune response. Expression of selected FA metabolism-associated, inflammatory and immune responsive genes was validated by qPCR. The inflammatory and immune responsive genes that are modulated by FA-associated D56 likely contribute to the enhanced resistance against V. vulnificus infection in Tilapia. Transcriptome profiling and filtering for two-fold change variation showed that 3795 genes were upregulated and 1839 genes were downregulated in D56-transgenic tilapia. These genes were grouped into pathways, such as FA metabolism, FA elongation, FA biosynthesis, biosynthesis of unsaturated FA, FA degradation, inflammation, immune response, and chemokines. FA-associated genes and immune-related genes were modulated by D56 at 6 h and 24 h post infection with V. vulnificus. The expression patterns of FA-related genes, inflammatory genes, antimicrobial peptide genes and immune responsive genes at 0, 3, 6, 12, 24 and 48 h post-infection suggests these genes are involved in the enhanced resistance of D56 transgenic tilapia to V. vulnificus.


Background
Tilapia (Oreochrombis niloticus) is an important commercial aquaculture species throughout the world, and its production is severely affected by the pathogenic bacteria Vibrio vulnificus, which causes septicemia in fish and humans [1][2][3][4]. Omega-3 polyunsaturated fatty acids (n-3 PUFAs) are known to exert beneficial effects, such as protection of liver, reduction of cholesterol, lower blood pressure and protect from cardiovascular diseases [5,6]. Furthermore, n-3 PUFAs show positive ionotropic effects and minimize tachyarrhythmia in animal models [7]. Many of these effects may be mediated by alterations in the proinflammatory cytokines, TNF-α, IL-1β, IL-6, prostaglandin (PG) E2, and PGF1α, which modulate the immune response in model organisms [8][9][10]. Dietary supplementation with eicosanoids and n-3 PUFAs is well documented to affect immune cell function and B-cell activity [11,12], and a recent report showed that PUFA-rich food limit pathogen infection in the aquatic organisms [13]. Similarly, transgenic expression of n-3 PUFA biosynthesis genes from Atlantic salmon, i.e., Fatty acyl desaturase synthase delta (Fadsd)5 and Fadsd6, in zebrafish limits infection with Vibrio alginolyticus and V. vulnificus [5,14].
Previously we reported the dual expression of SsFadsd5 and SsFadsd6 (D56) in tilapia [15]. The dual expression of these genes is under the control of a TRE-regulated CMV minimal promoter, which drives expression of D56 in liver and muscle [15]. Expression of D56 in tilapia enhances resistance to V. vulnificus infection [15]. In addition, the D56 transgenic tilapia exhibit altered gut microbial profiles [15]. However, the underlying molecular mechanism involved in the resistance to V. vulnificus has not been studied using a transcriptomic approach.
We compared the liver transcriptomes between V. vulnificus-susceptible wild-type tilapia and D56 transgenic tilapia with enhanced resistance to the pathogen to reveal the particular genes responsible for the resistance [15,16]. The alterations in expression of key genes were identified by gene enrichment analysis with KEGG pathway and GO tools. We showed the involvement of fatty acid (FA)-associated genes and immunomodulatory genes in the development of resistance against V. vulnificus infection in tilapia.

Results
Expression of recombinant delta-6 desaturase and delta-5 desaturase alters the transcriptome in tilapia Wild-type and D56-transgenic tilapia were infected with V. vulnificus, and RNA was extracted from liver at 0, 6 and 24 h post-infection (hpi). Transcriptome sequencing of six groups of samples produced a total of 275,304,348 raw reads for wild-type and D56-transgenic tilapia. After filtering the data 48,315,226, 38,578,158 and 35,079,100 clean reads were obtained for wild-type tilapia fish at 0, 6 and 24 h infected samples, respectively, representing 92.12, 89.24 and 92.19% of raw reads (Table 1). Similarly, 37,449,898,49,652,212 and 50,987,302 clean reads (90.62, 88.48 and 91.06% of raw reads) were obtained for D56-transgenic tilapia for 0, 6 and 24 hpi samples, respectively (Table 1). A total of 42,622 unigenes were identified from the RNA-sequencing and filtered for two-fold change in expression between V. vulnificus challenged wild-type and D56-transgenic tilapia (Supplementary Figure S1). At the 0 h time-point, 3795 genes were upregulated and 1839 genes were downregulated in D56-transgenic tilapia (Fig. 1a). At 6 hpi, 4365 genes were upregulated and 1976 genes were downregulated (Fig. 1a). At 24 hpi, 4665 were upregulated and 2202 genes were downregulated (Fig. 1a). We could recognize the relevance of DEG between wild-type and D56transgenic tilapia at different time-points. We found that 1112 DEG existed at three time-points (Fig. 1b).
Fatty acid-associated genes are altered in D56-transgenic tilapia according to KEGG pathway analysis Gene enrichment analysis using the KEGG pathway database showed a total of 24, 30 and 33 pathways were affected in D56-transgenic tilapia at 0, 6, and 24 hpi, respectively ( Table 2, 3 4). Immediately after infection, altered expression of various FA-associated pathways, such as FA metabolism, FA elongation, FA biosynthesis, biosynthesis of unsaturated FA, and FA degradation were observed (Table 2). Differences in the FA degradation pathway were also observed between wild-type and D56-transgenic fish at 6 h post infection (Table 3), and the FA metabolism and FA degradation pathways were identified as differentially expressed at 24 hpi (Table 4).

D56-transgenic tilapia exhibit altered immune-related gene expression in GO analysis
The GO enrichment analysis showed a total of 28, 23 and 35 gene sets that were differentially expressed in . Biological function-associated GO terms, such as defense response to bacterium, angiogenesis, immune response, antigen presenting and presentation and inflammatory response genes were also altered in D56transgenic tilapia. Altered molecular function-related GO terms included iron ion binding protein, cytokine and chemokine activity (Table 5, 6, 7). Taken together, the GO analysis revealed that inflammatory genes, chemokine-associated genes, cytokine-associated genes, immune-related genes and iron binding protein genes are differentially regulated in D56-transgenic and wildtype tilapia after V. vulnificus infection (Table 5, 6, 7).
We selected target genes with significant fold change, immune-related annotation and higher FPKM value for follow up research (Supplementary File 1).

Ectopic D56 alters FA metabolism-related genes
Since the KEGG analysis showed FA-associated pathways are altered in D56-transgenic tilapia, the expression of selected FA pathway-associated genes was analyzed by real-time PCR at 0, 3, 6, 12, 24 and 48 hpi with V. vulnificus. Since the D56 transgenes (delta-6 desaturase and delta-5 desaturase) are associated with FA biosynthesis, we analyzed the expression pattern by qPCR at many time-points. Significant alterations in the expression of FA-associated genes were observed from the qPCR data ( Fig. 2). Notably, ApoA4b was downregulated by D56transgenic tilapia at 24 hpi (Fig. 2a). CPT1 was upregulated at 24 and 48 hpi (Fig. 2b), and PCK1 was PPARα was upregulated at 6 and 12 hpi, but it was downregulated at 24 hpi (Fig. 2e). These results showed that FA metabolism-related genes are altered in transgenic tilapia upon V. vulnificus infection.

Ectopic D56 alters pro-inflammatory cytokines and CFD in whole blood sample
We also measured gene expression in whole blood samples after challenge. Expression of cytokines, inflammatory factors and complement-related genes was analyzed by real-time qPCR. For inflammatory factors, cytokines and complement-related genes, we found that were several significant differences between wild-type and D56transgenic tilapia whole blood. TLR-5 was altered at 24 hpi; NF-κB2 was altered at 6, 12, 24 and 48 hpi; NF-κBI was altered at 0, 6, 12, and 24 hpi; IL-1β was altered at 24 and 48 hpi; C1qb was altered at 0, 3 and 24 hpi; CFD was altered at 48 hpi ( Fig. 7a-f). In whole blood samples, the expression level of CPT1 was also different at 3 and 48 hpi (Fig. 7g). According to these results, protective and immune-related genes are induced in transgenic tilapia upon V. vulnificus infection.

Discussion
Tilapia (Oreochromis niloticus) is a staple product of the aquaculture industry, with annual consumption exceeding 3.7 million metric tons, as of 2014 [17]. Presently tilapia are grown in fresh-water pond culture systems in approximately 125 countries [18]. Breeding programs were used to develop improved versions of tilapia with high biomass [18]. Availability of whole genome sequence and RNA sequence data in recent years has allowed a greater understanding of the genetic makeup and expression profiles of different strains or groups of tilapia [19]. In fresh-water and brackish-water cultures, tilapia is prone to infection with the aquatic bacterial pathogen, V. vulnificus, which severely threatens the tilapia production [4,18,20,21]. V. vulnificus is a halophytic Gram-negative bacillus-type bacterium that can cause skin lesions, soft tissue dysfunction, and sepsisinduced mortality in tilapia or people who consume raw fish containing this pathogen [22,23]. The V. vulnificus strain 93 U204 has been isolated from an infected tilapia fish and its genome was previously sequenced [3]. Regulation of gene expression plays a major role in an organisms defense against pathogens [16]. Sequencingby-synthesis on an Illumina RNA-sequence platform has become a widely applied method for comparative transcriptome analysis [24,25]. The primary sequence-bysynthesis data is contained in a multiplexed and mixed form with sequence information of all the groups in a Bcl file [26]. The complex Bcl form has to be converted into a readable form, such as FASTQ format prior to further analysis [14]. The comprehensive computationbased resource, Gene ontology (GO), is extensively used to analyze the large amounts of FASTQ converted reads obtained from transcriptome analysis [27,28]. Hence, we analyzed our transcriptome data with the KEGG pathway tool to search for major pathways altered in D56 transgenic tilapia compared to wild type. The major enriched genes in identified KEGG pathways and major genes or gene groups from GO analysis were measured by real-time PCR to confirm the regulation of gene expression by D56 transgenic tilapia fish compared to wild types. Several reports have shown that resistance to infection can be enhanced in tilapia. AMPs, such as TP3 and TP4 have been reported to decrease the bacterial counts of V. vulnificus [29,30]. In addition, the multifunctional growth factor, PGRN, has been reported to modulate the immune response and improve survival of zebrafish infected with V. vulnificus [20]. Similarly, the granulin peptide, GRN-41, has been reported to exert antibacterial function against V. vulnificus [31]. When tilapia are fed with Epinecidin (Epi)-1-expressing transgenic Artemia, the mortality rate caused by V. vulnificus infection is decreased [32]. These studies demonstrate the utility of AMPs in controlling V. vulnificus infection. D56-transgenic zebrafish and tilapia exhibit resistance to V. vulnificus infection [5,15]. The genes expressed in D56-transgenic fish, Atlantic salmon Fadsd5 and Fadsd6, play an important role in n-3 PUFA biosynthesis. Interestingly, exogenous FAs assimilated into the vibrio species affect the swimming motility, bacterial membrane structure, permeability, and virulence [33].
A comparative analysis of liver transcriptomes of wildtype and D56-transgenic tilapia infected with V. vulnificus was performed in this study. The multiplexed data of the six groups in the Bcl file were quality controlled, and sequence data for each group was multiplexed into FASTQ format. From the short reads, 88.48-92.19% were mapped to the RNA genome in all samples. Thus, nearly 90% of the transcriptome was mapped for the comparative analysis. The KEGG pathway analysis showed several FA-associated pathways in wild-type and D56-transgenic tilapia groups at 0 hpi with V. vulnificus ( Table 2). KEGG analysis also revealed that the drug metabolism pathway, amino acid pathways and Cytochrome 450 pathway were altered by D56 (Table 2). These results showed that the expression of D56 enhances FA biosynthesis, especially biosynthesis of unsaturated fatty acids, and it may also lead to the modulation of additional pathways. Several such pathways were also differentially activated between wild-type and D56-transgenic fish group 6 h and 24 h post infected with V. vulnificus (Table 3, 4). These data demonstrated the D56-mediated FA pathway may control several amino acid biosynthesis mechanisms in addition to stress response mechanisms related to pathogen resistance (Table 2, 3, 4).
The GO analysis showed that genes associated with inflammation, chemokine synthesis, iron homeostasis and immune response are altered by D56 expression in   Table 5). The iron binding genes are related to the AMP hepcidin [34][35][36]. Hence, it is possible that tilapia hepcidins (TH), such as TH1-5 and TH2-3, may be functionally affected by D56 expression. Previously, the anti-V. vulnificus activity of TH2-3 has been reported [37]. Thus, further investigations may explore the possibility that D56 regulates tilapia hepcidins. In addition, several AMPs that are secreted by tilapia fish are associated with innate immunity and immune modulatory functions. Hence, the FA-associated D56transgene expression may also directly alter the level of AMPs in tilapia fish. To address this possibility, the expression AMP-associated genes in tilapia was studied by qPCR (Fig. 4). TH and BP1 were upregulated at 24 h and 48 hpi with V. vulnificus in D56-transgenic tilapia, suggesting the activation of these peptides was promoted by D56 to combat V. vulnificus infection (Fig. 4). The complement genes, C1qb, CFHR1 and CFD were also upregulated in D56-transgenic tilapia after V. vulnificus infection (Fig. 3). It is also possible that the plasma antibacterial function is stimulated by D56 expression in transgenic tilapia [38,39]. The inflammation-associated gene expression included early upregulation of TLR-2 and TNF-α in D56-transgenic fish (Fig. 5). The expression level of NF-κB did not vary, however, it is possible that the translocation of NF-κB from cytoplasm to nucleus triggered the activation of downstream genes, such as TNF-α [40]. The CD8-related Peroxiredoxin (PRDX)1, G-protein receptor associated Atypical Chemokine Receptor (ACKR)4, Tissue inhibitor of Metalloprotease (TIMP)2 associated with extracellular matrix were also altered by the presence of D56 (Fig. 6). Probing the immunomodulatory gene expression in whole blood revealed that genes such as NF-κB2, TLR-5, IL-1β, CFD and C1qb were upregulated by D56 after V. vulnificus infection (Fig. 7). Altogether, the FA-associated pathways triggered in the D56-transgenic tilapia appear to regulate a variety of immune-related genes that may serve to enhance resistance to V. vulnificus infection.

Conclusions
We compare the results between NGS and qPCR (Supplementary Table S2). We found that TLR-2, TLR-5 and ACKR4 were induced in the D56-transgenic line. TLR-5 and TLR-2 may stimulate a pro-inflammatory response, and ACKR4 blocks chemokine signaling. Through the regulation of CFHR1 and CFD, n-3 PUFAs can modulate the complement system. Moreover, AMP-associated genes, including TP3, TP4, TP5, TH, PGRN, BPI and LEAP2, may possibly be regulated by n-3 PUFAs (Fig. 8a).
In terms of downstream effects, HNF4A and PPARs may induce CPT1 and PCK1. This action would enhance ATP production from FA oxidation and help the host overcome infection by pathogenic bacteria (Fig. 8b). Furthermore, NF-κB downstream genes, such as TNF-α and IL-1β, may suppress pathogens (Fig. 8c). In this study, we measured the transcriptome-wide responses of wildtype and D56-transgenic tilapia to V. vulnificus infection. Several pathways associated with inflammation, immune response, chemokine and cytochrome were specifically altered in D56-transgenic tilapia, as revealed by KEGG pathway and GO analysis of the wild-type and D56transgenic tilapia at the baseline (0 hpi) and after infection (6 and 24 hpi). These results suggest that D56 may modulate pro-inflammatory cytokines, AMPs and immune response genes to enhance resistance to V. vulnificus (Fig. 8).

Tilapia fish and culture
Tilapia fish Oreochromis niloticus were acclimatized in FRP tanks with 2000 l capacity under the controlled conditions of 28°C with a 13 h light/ 11 h dark cycle with Alanine rich food at 4% body weight per day for 45 days. Wild-type and transgenic tilapia fish with dual expression of delta-6 desaturase in liver and delta-5 desaturase in muscle (D56) were gift from Dr. Jen-Leih Wu's lab from Institute of Cellular and Organismic Biology, Academia Sinica, Taipei, Taiwan [15]. After 45 days the wild-type tilapias reached an average of 8.93 ± 0.98 cm length and 11.52 ± 3.55 g weight. The transgenic tilapias were 10.06 ± 1.02 cm in length and 16.25 ± 3.61 g weight. The experimental fish were handled after getting approval from the Institutional Animal Care and Use Committee (IACUC) at Academia Sinica, Taiwan, and the IACUC guidelines were followed. According to the formula: n = (Z_1-beta * sqrt(p1q1 + p2q2) + Z_1-alpha/2 * sqrt(2 * p_avg * q_avg))^2 / delta^2, p2 = p1 + delta, p_ avg = (p1 + p2) / 2, q_avg = 1 -p_avg. The minimum number of samples in each group was 5 (K = 2, α = 0.05, 1-β = 1-0.2 = 0.8, p1 = 37% and Δ = 74%). A total of 66 animals were used in this study (5 samples * 2 groups * 6 time-points = 60 for qPCR analysis; 6 animals were used for transcriptome analysis). The method for euthanasia was rapid chilling by soaking the animal in ice for 10 min. The ratio of ice to water was roughly 5:1. Ice and water were both measured volumetrically for a total volume of roughly 10 L [41].

Expression of D56 in transgenic Tilapia
The vector construct and expression of D56 has been explained previously [15]. In brief the delta-6 desaturase transgenic line driven by muscle-specific CKMb promoter express the Fadsd6 with TcFP11 reporter was created by co-injecting the plasmid into the embryo. Similarly delta-5 desaturase transgenic line driven by liver-specific Fabp10 promoter to express the Fadsd5 with TcFP13 reporter was created separately. The delta 5 and delta 6 transgenic tilapia have been crossed and the offspring carrying D56 has been selected for the experiment. Atlantic salmon (Salmon salar) delta-6 desaturase gene Fadsd6 and delta-5 desaturase gene Fadsd5 were used in the construct.

V. vulnificus culture and infection in to tilapia
The pathogen V. vulnificus strain 93 U204 was initially cultured in a TCSB (Thiosulfate Citrate Bile Salts Sucrose, BD Difco™) agar plate for 16 h at 28°C [3]. From overnight cultured plate a single colony was picked and cultured in 3 ml TSB (Tryptone Soy Broth, BD Difco™) with 1.5% NaCl in an incubated shaker at 28°C with 200 rpm up to four hours until the OD 600mm reached 0.7 to 0.9. At the end of acclimatization the wild-type and transgenic fish were subjected to V. vulnificus challenge by i.p. injecting 50 μl of diluted TSB containing 1.9 × 10 2 CFU of V. vulnificus for each fish. Thirty fishes in each group was infected with V. vulnificus and five fish per each group was sacrificed at 0, 3, 6, 12, 24 and 48 h post infection for molecular analysis.

Extraction of total RNA
The fish was sacrificed at the end of each time point of the experiment and liver tissue were collected for total RNA extraction. The liver tissue was homogenized and total RNA was extracted using Trizol® Reagent by following the company protocol (Invitrogen, USA). Purified RNA was quantified at OD260nm using a ND-1000 spectrophotometer (NanodropTechnology, USA) and   Figure S1). Functional enrichment assay in differentially expressed genes of each experiment design was performed using clusterProfiler v3.6. Genes with low expression level (< 0.3 FPKM value) in either or both of the treated and control samples were excluded (Supplementary File S1). Genes with p value ≤0.05 and ≥ 2-fold changes were considered significantly differentially expressed.
Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) analysis GO enrichment analysis of the differentially expressed genes (DEGs) was implemented using the clusterProfiler package. Pathway enrichment was used for KEGG analysis (https://www.genome.jp/kegg/), which uses public databases to examine Oreochromis niloticus biological pathways. The datasets generated and analyzed during (See figure on previous page.) Fig. 7 Delta 5 and Delta 6 transgenic (D56) tilapia fish modulates the pro-inflammatory cytokines and immune-related genes in whole blood. Wild-type and D56 transgenic tilapia fish challenged with V. vulnificus and whole blood samples from 0, 3, 6, 12, 24 and 48 h post infected conditions were collected for qPCR. Expression of a TLR-2, b IL-1β, c NF-κB2, d NF-κBI, e CFD, f C1qb, g CPT1 genes relative to EF-1α was estimated. Values represented as Mean ± SEM (n = 5). Significance was determined by T-TEST (*P < 0.05,**P < 0.01,***P < 0.001) Fig. 8 Proposed mechanism of Delta 5 and Delta 6 transgenic (D56) mediated ω-3 fatty acid metabolism to enhance resistance against V. vulnificus in Tilapia. a The ω-3 fatty acids synthesized in D56 transgenic fish regulates or activates the extra cellular acting agents such as Antimicrobial peptide (AMP), Complement system and ACKR4 by altering their expression. b In cytosol D56 mediated ω-3 fatty acid regulates HNF4A, TLR-5 and TLR-2, PPARs, TIMP2 and PRDX1. c In nucleus transcription factors such as NF-κB is altered to regulate the proinflammatory cytokines. (Image from ourselves) the current study are available in the NCBI repository, https://www.ncbi.nlm.nih.gov/bioproject/705417 (Accession: PRJNA705417 ID: 705417) (Accession: PRJNA705417 ID: 705417).
Gene expression analysis with qPCR cDNA was obtained from 0.5 μg of total RNA using ReverTra Ace® qPCR RT Master Mix with gDNA Remover kit (TOYOBO, Japan). qPCR was performed using Applied BiosystemsStepOnePlus™ System (ABI, USA) machine. The qPCR reaction mix comprise of 5 μl SYBR® Green Realtime PCR Master Mix (TOYOBO, Japan), 0.5 μl (10 μM) each of gene specific primers (Supplementary Table S1), and 4 μl of twenty times diluted cDNA. Amplification was performed with an initial 95°C for 1 min and 40 cycles of 95°C for 15 s, 60°C for 15 s and 72°C for 45 s (Supplementary Figure S2). Gene expression relative to EF-1α was estimated by 2 T -△△C method as described [42].