Transcriptome changes in response to temperature in the fish pathogen Photobacterium damselae subsp. damselae: Clues to understand the emergence of disease outbreaks at increased seawater temperatures

The marine bacterium Photobacterium damselae subsp. damselae (Pdd) is a generalist and facultative pathogen that causes disease in a wide range of marine animals including fish species of importance in aquaculture. Disease outbreaks in fish farms have been correlated with an increased water temperature during summer months. In this study, we have used RNA sequencing to analyze the transcriptome of Pdd RM-71 cultured at two different temperatures, which simulated temperature conditions experienced during free swimming lifestyle at mid latitudes in winter months (15°C) and during outbreaks in aquaculture in warm summer months (25°C). The enhanced bacterial growth of Pdd observed at 25°C in comparison to 15°C suggests that an elevated seawater temperature contributes to the build-up of a sufficient bacterial population to cause disease. In comparison to growth at 15°C, growth at 25°C resulted in the upregulation of genes involved in DNA synthesis, nutrient uptake, chemotaxis, flagellar motility, secretion systems and antimicrobial resistance. Plasmid-encoded virulence factors, which include a putative adhesin/invasin OmpU, a transferrin receptor and a serum resistance protein, were also upregulated. Transcription factor RpoS, genes involved in cold shock response, modulation of cell envelope and amino acid metabolism, as well as genes of yet unknown function were downregulated at 25°C. Notably, the gene encoding damselysin cytotoxin (Dly) was among the most highly transcribed genes at the two assayed temperatures, at levels comparable to the most highly expressed housekeeping genes. This study contributes to our understanding of the regulatory networks and biology of a generalist marine bacterial pathogen, and provides evidence that temperature regulates multiple physiological and virulence-related functions in Pdd.


Introduction
The Vibrionaceae family includes aquatic bacteria found in ocean environments and temperature and salinity are the main abiotic factors that shape their abundance and distribution [1,2]. Rising ocean temperatures contribute to the increase in prevalence and severity of a wide range of Vibrio-related diseases in marine organisms and also in humans [3,4]. Vibrio and Photobacterium infections have a marked seasonal distribution, and most cases occur during the warmer summer months [5].
Photobacterium damselae subsp. damselae (hereafter Pdd), a marine bacterium of the Vibrionaceae family, is a pathogen of a wide range of marine animals, including fish, molluscs, crustaceans and cetaceans. It is also an emerging pathogen in marine aquaculture systems that causes wound infections and septicaemia in economically important fish species [6][7][8]. In addition, it is an opportunistic human pathogen which causes severe wound infections that can have a fatal outcome [9].
The major reported virulence factors of Pdd are cytotoxins with hemolytic activity [8]. Highly virulent strains harbor the virulence plasmid pPHDD1 that carry the cytotoxin genes dly and hlyA pl [21]. dly encodes damselysin toxin (Dly), a phospholipase-D active against sphingomyelin [22], and hlyA pl encodes the pore-forming toxin phobalysin P (PhlyP) [23]. In addition, the pore-forming toxin phobalysin C (PhlyC), encoded by the hlyA ch gene located on chromosome I, and the phospholipase PlpV also contribute to virulence in fish and to cell toxicity [24,25]. These four cytotoxins are secreted via a type II secretion system [26]. The two-component regulatory system RstAB positively regulates transcription of the hemolysin genes dly, hlyA pl and hlyA ch , and its inactivation severely impairs virulence [27].
Despite the increasing evidence suggesting that outbreaks caused by Pdd in fish farms are triggered by rises in seawater temperature during summer months, the role of temperature in the physiology and gene regulation of this pathogen has not been studied so far. The aim of this study was to investigate the transcriptome of Pdd at two different temperatures, 15˚C and 25˚C, using an RNA-seq approach. We identified a number of genes upregulated at 25˚C-the temperature at which most outbreaks occur in aquaculture farms-that likely contribute to an increased bacterial growth rate and to an enhanced ability to colonize and survive in fish hosts. Additionally, we found genes with higher expression at 15˚C which might aid Pdd to adapt to life in colder waters during winter months. The global transcriptome data also shed light on the relative expression values of genes encoding virulence factors compared to housekeeping genes and demonstrated that the damselysin toxin gene is one of the most highly expressed genes in the cell. Finally, this study has brought to the forefront many previously overlooked genetic networks and gene clusters of this pathogen. Overall, the information generated in this study is expected to provide novel approaches for the prevention and control of vibriosis caused by Pdd in marine fish aquaculture.

Growth analysis
Cells were routinely grown at 15 or 25˚C on tryptic soy agar (TSA) or in tryptic soy broth (TSB) supplemented with NaCl up to 1% (TSA-1 and TSB-1, respectively). For growth curves, three replicates for each temperature of the assay (15 and 25˚C) were grown in TSB-1 until obtaining exponentially growing precultures (OD 600 : 0.3). Then, 1:100 dilutions of each preculture were grown in 100 μl of TSB-1 in 96 well plates and the optical density (OD 600 ) was measured during 48h using the spectrophotometer Epoch2 microplate reader (Biotek). Three replicates were measured for each temperature condition.

RNA-seq
Experimental design, RNA extraction and purification. As for RNA-seq approach, 3 biological replicates were performed for each condition. 15˚C was chosen as the control condition and is close to the temperature that this bacterium finds during free swimming lifestyle in mid latitudes, whilst 25˚C represents the higher temperature condition that usually precedes aquaculture outbreaks. For each temperature, three independent precultures were started and grown until an OD 600 : 0.3. Then, 1:100 dilutions of each preculture were grown in 10 ml of TSB-1 in 100 ml flasks until they reached a sharp OD 600 of 0.55. Cells were immediately treated with RNAprotect Bacteria Reagent (Qiagen) for stabilization of RNA following manufacturer's instructions. Pelleted cells were then carefully resuspended in TE buffer (30mM Tris�Cl, 1mM EDTA, pH 8.0) containing 15mg/ml lysozyme (Sigma Aldrich) and the appropriate volume of Proteinase K (Qiagen). RNA extraction was subsequently carried out using RNeasy Mini Kit (Qiagen) following manufacturers' instructions. An extra DNase I treatment was carried out using the on-column kit RNase-free DNase (Qiagen) to eliminate genomic DNA contamination. RNA was eluted using nuclease-free water. The quality and the quantity of the total RNA was determined using a Bioanalyzer 2100 (RNA 6000 Nano chip assay) and a Qubit 3.0 (Quant-It dsRNA BR Assay).
Libraries preparation and sequencing. Total RNA was rRNA-depleted using the Ribo-Zero rRNA Removal Kit (Gram Negative Bacteria) (Illumina) and cDNA libraries were obtained using the TruSeq RNA kit following Illumina´s recommendations. Briefly, rRNAdepleted RNA was chemically fragmented prior to reverse transcription and cDNA generation. The cDNA fragments then went through an end repair process, the addition of a single 'A' base to the 3' end and then ligation of the adapters. Finally, the products were purified and enriched by PCR to create the indexed final double stranded cDNA library. The pool of libraries was sequenced on an Illumina HiSeq 2500 sequencer.
Mapping and quantification of transcripts. The quality control of the raw data (raw reads) was performed using the FastQC [http://www.bioinformatics.babraham.ac.uk/projects/ fastqc/] program. The raw pair-end reads were first mapped against the reference genome of the Photobacterium damselae subsp. damselae type strain CIP102761 (GenBank Acc. No. NZ_ADBS00000000.1). Reads that did not map to the reference genome (corresponding to genes of RM-71 strain not present in CIP102761) were subsequently mapped to the draft genome sequence of strain RM-71 (GenBank Acc. No. NZ_LYBT00000000.1). The two processes were completed using the Bowtie2 [28] v2.2.6 algorithm. Several quality control steps were performed. Reads displaying a very low quality were removed by using Samtools [29] and Picard Tools software [30]. Furthermore, one of the key factors that can condition the sequencing process is the GC content of samples which was checked as normal (distribution between 40-60%) in our experiment. Likewise, distribution of duplicates was evaluated to confirm the normal small proportion. The process of genetic quantification was carried out by the HTSeq [31] software (0.6.1 version).
Comparison between samples. Concordance between samples of the same condition (replicates of each of the two assayed temperatures) was carried out by a study of correlation and distance considering the whole transcriptome normalized by the size of the library. This process was made using the statistics program R. Differential expression analysis was assessed using DESeq2 [32] method (1.18.1 version). The analysis of Differentially Expressed Genes (DEG) was done by using statistical packages designed by Python and R, using the DESseq2 [32] algorithm applying a differential negative binomial distribution for the statistics significance. Comparison between the two different conditions (25˚C vs. 15˚C) was set as fixed effect in DEseq2. A Python script developed at Sistemas Genómicos (Valencia, Spain) was employed to generate a data matrix for each group condition with the counts obtained from HTSeq count for each sample (each of the three replicates at each of the two temperatures). We considered differentially expressed genes those with Fold Change (FC) value lower than -1.5 or higher than 1. 5  were used to represent the expression of each individual gene. FPKM is used for normalization of the data since it indicates the number of lectures of a given gene per kilobase (independently of the length of the gene), and per million reads (independently of the size of the library).

Growth dynamics of Pdd at 15 and 25˚C
Pdd RM-71 was the strain selected for the present study. It was isolated during a disease outbreak in a turbot (Scophthalmus maximus) farm in Galicia (NW Spain), when the water temperature increased suddenly from 18˚C to 22-24˚C in the summer of 1988 [10]. It is a highly virulent, strongly hemolytic and cytotoxic strain and contains the pPHDD1 virulence plasmid [21]. Growth of RM-71 was analyzed in 48-h continuous cultures at 15 and 25˚C in TSB-1 medium (S1 Table). The two assayed temperatures simulated an a priori non-risky condition (15˚C) and warm water episodes that trigger aquaculture outbreaks (25˚C). The beginning of the exponential phase was delayed at 15˚C compared to growth at 25˚C and there was a great difference between optical density values after 15 h of cultivation at 15˚C (OD 600 : 0.129) and 25˚C (OD 600 : 0.527) (Fig 1). These observations suggest that 25˚C is closer to the optimal growth temperature of Pdd than 15˚C.
This substantial increase in bacterial proliferation at 25˚C with respect to 15˚C during the first 15 h of cultivation might contribute to the rapid progression of Pdd outbreaks when the sea water temperature increases. Proliferation of species of the Vibrionaceae family is favoured by warm (>15˚C) sea waters [35], and recent studies have demonstrated that following an increase in water temperature, Vibrios can go from barely detectable to being the predominant bacteria in a very short time [36,37]. Pdd isolation was first reported from ulcers in damselfish (Chromis punctipinnis) during the summer and fall seasons in southern California and was shown to follow a seasonal pattern of infectivity. It was proposed that elevated water temperatures might allow the build-up of sufficient bacterial populations to cause disease in damselfish, hence the seasonal infectivity of Pdd [38]. Previous studies have provided sound evidence that seawater transmits the disease caused by Pdd and that the spread of this bacterium largely depends on water temperature [39]. Skin is suggested to be a potential route of penetration for this pathogen, which is able to specifically adhere to fish mucus [39]. Hence, even if fish are colonized by a small number of bacteria, fast proliferation enhanced by warm temperature will facilitate the infecting bacterial population to evade host immune responses by a variety of mechanisms [40]. High numbers of bacterial cells at the infection site might cause exhaustion of complement components as well as of phagocytes [41] leading to systemic infection and fish death. Overall, the results of the growth dynamics analysis at 15 and 25˚C contribute to understand why increased water temperatures precede most outbreaks caused by this pathogen.

RNA sequencing results
Strain RM-71 was grown at 15 and 25˚C, and cDNA prepared from mRNA isolated from cultures at the two different temperatures was subjected to Illumina sequencing. Around 55 to 71 million reads were generated for each biological replicate (S2 Table). Growth at 15˚C was defined as the control condition. The comparative analyses of the transcriptomes at 15 and 25˚C resulted in a total of 1195 differentially expressed genes (DEGs): 641 genes with lower expression at 25˚C (FC lower than -1.5) and 554 genes with higher expression at 25˚C (FC higher than 1.5) (Fig 2, S3 and S4 Tables).
Similar to other members of the Vibrionaceae family, Pdd RM-71 contains two chromosomes. In addition, this strain harbors the virulence plasmid pPHDD1 [21]. Using the complete sequences of chromosome I and II of the type strain CIP102761, and the complete pPHDD1 plasmid sequence of strain RM-71 (GenBank Acc. No. NC_014653) as references, we distributed the DEGs into each replicon. Notably, we observed an imbalance in the number of DEGs between the two chromosomes (Fig 3). In chromosome I similar numbers of DEGs are upregulated and downregulated. However, chromosome II contains 204 downregulated genes and only 77 upregulated genes at 25˚C. The chromosome I of Vibrionaceae species contains most of the essential genes, whereas chromosome II has a more flexible gene content and is responsible for adaptation to environmental changes [42,43]. Interestingly, among the 31 DEGs in the virulence plasmid pPHDD1, 24 corresponded to genes whose expression is upregulated at 25˚C, an observation of special interest since this plasmid constitutes a hallmark of highly virulent isolates [8]. Indeed, the present study has unveiled potential virulence factors

Genes involved in growth and virulence are upregulated at 25˚C
Growth at 25˚C resulted in the upregulation of 533 genes that mapped to the genome of the type strain (S3 Table) and of 21 additional genes unique to strain RM-71 (S4 Table). A list of the 50 top DEGs upregulated at 25˚C plus additional selected genes is shown in Table 1.

Fig 2. Smear plot of differentially expressed genes (DEGs) in
Pdd RM-71 exposed to two temperatures, 15˚C and 25˚C. The smear plot shows the relationship between the log FC and mean of normalized counts. Grey points represent genes with non-significant changes in expression, whereas red points represent genes that are significantly differentially expressed.
https://doi.org/10.1371/journal.pone.0210118.g002 Nutrient acquisition and metabolism. Genes encoding membrane proteins, nutrient transporters and porins were upregulated at 25˚C ( Table 1). The nucleoside permease NupC (VDA_001789) was the most upregulated gene at 25˚C. A V. cholerae nupC deletion mutant was impaired for nucleoside acquisition leading to diminished fitness in nutrient-limited environments [44]. Growth at 25˚C upregulated two ribonucleotide reductases belonging to class Ia and II, one of them, VDA_001894, is coenzyme B12-dependent. In accordance, the vitamin B12 receptor BtuB was also upregulated ( Table 1). Upregulation of the uracyl phosphoribosyltransferase VDA_001405, an enzyme necessary for the synthesis of precursors of all pyrimidine nucleotides, was also observed. Among upregulated amino acid biosynthesis enzymes and aminotransferases was the glutamine synthetase type I, which has a central role in amino acid biosynthesis. The putative trypsin superfamily protein encoded by A0J47_09785, which has a possible function in peptide degradation, as well as a number of ribosomal and translation-related proteins were upregulated suggesting that growth at 25˚C enhances protein synthesis and renovation.
Pdd degrades extracellular lipids, which may serve as carbon and energy sources [10,11,25]. An operon encoding an extracellular lipase (VDA_001610) and a fatty acid transporter https://doi.org/10.1371/journal.pone.0210118.g003 Table 1. List of selected DEGs with enhanced expression at 25˚C including the 50 top upregulated genes at 25˚C. Enhanced expression at 25˚C is denoted by positive FC values. Genes with VDA codes correspond to the annotation in the CIP102761 genome, and genes with A0J47 codes correspond to the annotation in the RM-71 genome. FadL (VDA_001611) were 2-fold upregulated at 25˚C, and so was the long-chain fatty acid transporter VDA_002568 ( Table 1), suggesting that exogenous lipid degradation and uptake of fatty acids might constitute an advantage for Pdd fast replication at 25˚C. Iron acquisition plays a role in Pdd virulence for fish [45][46][47]. Our analysis unveiled upregulation at 25˚C of a pPHDD1 plasmid-borne gene that encodes a putative transferrin binding protein (VDA_000111) (see below), and of VDA_000794 encoding a TonB-dependent siderophore receptor. The ligand(s) transported through VDA_000794 are unknown, but a recent study demonstrated that expression of this gene is enhanced in Pdd under iron-limitation conditions [48].

Gene ID Product/Function Fold Change p-value Location
Motility and chemotaxis. Motility and tissue colonization constitute important factors in Pdd pathogenicity [39,49]. Four flagellum-related genes were found among the 50 most upregulated genes (Table 1). Additional upregulated genes included flagellar hook protein FlgE (VDA_002616), flagellar motor rotation protein MotB (VDA_002670), flagellar hook-associated protein FlgK (VDA_002609) and flagellar basal-body rod modification protein FlgD (VDA_002617) (S3 Table). Previous studies have reported that Pdd can infect new hosts through seawater and that increased water temperature boosted infection by this route [49]. The upregulation of motility-related genes at 25˚C supports these previous observations and surely calls for further studies along these lines.
Chemotaxis is initiated by membrane chemoreceptors dubbed methyl-accepting chemotaxis proteins, which bind ligands and transduce a signal cascade that modulates flagellum activity. Notably, two chemotaxis-related genes were found among the 50 most upregulated genes at 25˚C and correspond to the methyl-accepting chemotaxis proteins VDA_003044 and VDA_001198 (Table 1). A recent study showed that mutants in chemotaxis genes are not only impaired in swimming motility in Pdd but also exhibit diminished production of the major virulence factor PhlyP and impaired adhesion to eukaryotic cells [50]. Vibrio fischeri and V. anguillarum mutants in chemotaxis functions also display a high reduction of virulence in fish and the inability to colonize host tissues [51,52]. Collectively, these results suggest that growth of Pdd at 25˚C could enhance chemotaxis and flagellum-dependent motility, contributing to access and adhesion to fish hosts and increasing chances of outbreaks in aquaculture farms during summer months. Stress response. A number of chaperones, heat shock and stress-related proteins were listed among the 50-top DEGs ( Table 1), suggesting that growth at 25˚C constitutes a heat stress condition. Upregulation was also found for DnaJ chaperonin (VDA_002770), heat shock protein GrpE (VDA_002772) and heat-shock chaperonin VDA_003125 (S3 Table). Upregulated proteases included the outer membrane stress sensor protease DegQ and the ATP-dependent protease HslV, two peroxidases and the iron-sulfur cluster assembly protein SufB. SufB synthesizes Fe-S clusters that act as cofactors in cellular processes under conditions of iron starvation or oxidative stress in E. coli [53]. Altogether, it is conceivable that an increase in temperature might constitute a signal for Pdd to activate its molecular machinery against reactive oxygen species formation by host cells, therefore linking temperature rise with virulence for fish hosts.
Transcriptional regulators and intracellular signalling. A DeoR family transcriptional regulator and a XRE family regulator were listed within the top-50 DEGs. Although studies about these transcriptional regulators are scarce, in Shigella flexneri and Salmonella typhi DeoR regulators are important for virulence and for intracellular growth [54, 55]. Thus, these two genes may regulate processes related to Pdd pathogenicity and further investigation into these regulators is required. Cyclic-di-GMP is an intracellular second messenger involved in environmental signalling and regulates a number of phenotypes in bacteria. VDA_002825, encoding a cyclic-di-GMP phosphodiesterase A, was 3.67-fold up-regulated at 25˚C. In V. cholerae cyclic-di-GMP regulates motility and biofilm formation [56,57] and DNA repair [58], among other functions.
Secretion systems. The Pdd type II secretion system (T2SS) plays a major role in the secretion of the four cytotoxins Dly, PhlyP, PhlyC and PlpV [25,26]. The complete cluster of eps (extracellular protein secretion) genes (VDA_003114-VDA_003123), encoding part of the T2SS machinery, was moderately upregulated at 25˚C relative to 15˚C (S3 Table), suggesting that the T2SS secretome serves important functions when this organism is growing at the temperature conditions that enhance outbreaks in aquaculture farms.
Potential virulence factors encoded within pPHDD1 plasmid are upregulated at 25˚C. A recent study reported that a double mutant of RM-71, with deletion of Dly and PhlyP-encoding genes, was still more virulent in a sea bass fish model than a naturally plasmidless strain [25], suggesting that pPHDD1 encodes additional yet uncharacterized virulence genes. The analysis of the differential gene expression profiles along pPHDD1 at the two temperatures brought to the forefront a collection of plasmid-encoded genes which have putative roles in virulence and that surely will deserve special attention in future studies (Fig 4, S5  Table).
A cluster of 5 upregulated genes (VDA_000154 to VDA_000158) encode the TolC protein and AcrAB, plus two additional proteins. TolC and AcrAB form a tripartite multidrug and toxin secretion efflux pump. Expression of tolC and acrAB in the fish pathogen Yersinia ruckery is increased at 28˚C with respect to 18˚C with a concomitant increased resistance to antibiotics and toxic substances such as acriflavine [59]. The role of this system in the biology of Pdd remains unknown. In addition to a possible role in efflux of toxic substances, we cannot rule out the possibility that this system participates in the secretion of virulence factors, as TolC is part of the type I secretion system that exports hemolysins and other virulence factors in several Gram negative pathogens [60].
VDA_000110 is 36% identical to V. vulnificus Vep07, an outer membrane lipoprotein that confers resistance to eel serum. Mutation of vep07 caused the loss of virulence in eels in V. vulnificus biotype 2 [61]. Of note, V. vulnificus Vep07 is encoded within pVvbt2, a transferable virulence plasmid [62]. VDA_000113 is homologous to Vibrio tasmaniensis OmpU, a protein with an essential role in adhesion to and invasion of mollusc host cells [63]. OmpU homologues in other Vibrio species also play major roles in host-cell recognition and pathogenesis [64,65]. Finally, the protein encoded by VDA_000111 shares 60% identity with the V. vulnificus plasmid-encoded Vep20 protein, a transferrin receptor with a major role for virulence in eels [66]. A previous study has demonstrated that Pdd strains can multiply in the presence of transferrin as the sole iron source [46]. Surely, the role of this gene deserves further investigation as it may play a main role in iron acquisition from host transferrin.
The three genes A0J47_18120, A0J47_18115 and A0J47_18110 are exclusive of RM-71 and have no known homologues in other pPHDD1-containing strains of Pdd studied so far. These genes are predicted to participate in defence against other bacteria which occupy the same ecological niches and compete for the same resources. A0J47_18120 is a proline-alanine-alaninearginine (PAAR) domain-containing protein. Homologues of these proteins in V. cholerae are essential for secretion through the type VI secretion system (T6SS) and killing of target cells [67]. A0J47_18115 encodes a putative RNase toxin dubbed Ntox44 which is found in many bacterial groups and is predicted to be exported by the type II, type VI or type VII secretion systems [68].

Growth at 25˚C causes downregulation of functions related to cell envelope, metabolism and stress response
Albeit 15˚C being far from the optimal growth temperature (25˚C) of Pdd in laboratory conditions, it is conceivable that this bacterium lives in marine ecosystems in vast areas of the globe at temperatures below its optimum. Indeed, seawater temperatures during the summer of 1988 previous to the heat wave that caused a Pdd outbreak in a turbot farm in Galicia were around 18˚C [10], and seawater temperatures in the same geographical area during the winter months are known to fluctuate between 13 and 15˚C [69]. Growth at 25˚C resulted in the downregulation of 614 genes that mapped to the genome of the type strain (S3 Table) and of 27 additional genes unique to strain RM-71 (S4 Table), whose changes in expression are denoted by a negative FC value. A list of the top-50 downregulated genes at 25˚C, plus additional genes, is shown in Table 2 and reveals an important number of loci organized in operons.
The most important change is experienced by a putative operon constituted of VDA_001578, VDA_001579 and VDA_001580 (Table 2). VDA_001579 is the most differentially expressed gene in the whole transcriptome of this pathogen in this study, and has no known homologues with a demonstrated function so far. VDA_001580 is a predicted serine protein kinase PrkA, a family of proteins involved in cell wall homeostasis [70], saline stress, motility [71] and virulence [72]. VDA_001578 is predicted to be a member of the Stage V sporulation protein SpoVR family. SpoVR confers resistance to Bacillus subtilis spores and it has been hypothesized that homologues in other species might play a role in peptidoglycan synthesis regulation [73]. Also related to the cell wall is VDA_002031, which encodes a LrgB-family protein, a group of enzymes responsible for modulation of murein hydrolase activity [74]. VDA_001762 encodes lysophosphatidic acid acyltransferase PlsC, an integral membrane protein involved in phospholipid biosynthesis [75]. Adjustment of membrane composition is a conserved strategy that bacteria use to face variations in environmental parameters [76,77]. The importance of the cell envelope in Pdd acclimatization to temperature changes surely will deserve further investigation.
Growth at 25˚C downregulated the expression of peptidases, membrane transporters and metabolic enzymes ( Table 2). Of note is the downregulation of the oligopeptide permease system oppABCDF whose main function is predicted to be nutritional [78]. Interestingly, the expression of oppABCDF in Vibrio alginolyticus was found to be sensitive to temperature [79]. The role of the opp operon in Pdd physiology has not been studied so far and further evaluation of these genes is necessary to assess their possible contribution to bacterial fitness at low temperatures. VDA_001251 contains a NifB/NifX domain for synthesis of iron-molybdenum cofactors. These cofactors bind the active site of dinitrogenase enzyme which participates in nitrogen fixation [80]. Although nitrogen fixation has not been studied in Pdd, some members of the Vibrionaceae family have this ability [81]. Enzymes of amino acid metabolism were also downregulated at 25˚C. The upregulation of specific amino acid biosynthetic pathways at both temperatures might be in part attributable to the need to provide optimal amino acid ratios for the different types of proteins produced at each temperature.
Genes of the Krebs cycle were downregulated at 25˚C: citrate synthase, malate synthase and isocitrate lyase. In addition to its role in the cell bioenergetics, the Krebs cycle is connected to iron acquisition in Pdd where endogenous citrate is used as iron scavenger [45,82]. Thus, the use of ferric citrate as an iron source might be favoured in Pdd at low temperatures.
Biosynthesis and uptake of betaine glycine constitutes an adaptation for growth at low temperatures and is part of the cold-stress response of the marine fish pathogen Vibrio anguillarum [83]. Of note, a glycine betaine transporter showed a strong downregulation at 25˚C in our study (Table 2). Only one gene among the main 50 downregulated genes at 25˚C encoded a cold shock protein (VDA_003169). In contrast, some mesophilic bacteria such as V. cholerae [84] and V. parahaemolyticus [85] overexpress cold shock proteins following shifts to low temperatures close to 15˚C. The alternative sigma factor RpoS (VDA_000629) is downregulated at 25˚C. RpoS is a major regulator of the general stress response pathway in bacteria [86]. Important traits regulated by RpoS include virulence and colonization in V. cholerae and V. parahaemolyticus [87,88]. A three-gene operon encoding a glutamate decarboxylase, a glutaminase and a glutamate/GABA antiporter is potentially involved in acid resistance [89]. The biological roles of this system in Pdd are unknown, but amino acid decarboxylation has been reported in other species of the Vibrionaceae family as a strategy for acid tolerance [90]. An operon which encode the methionine sulfoxide reductase system MsrPQ was downregulated in Pdd at 25˚C. The MsrPQ systems participate in the repair of oxidative damage [91].

Growth at 25˚C does not upregulate expression of the major cytotoxins of Pdd: Dly cytotoxin is within the 10 most expressed genes at 15˚C
Considering their importance in the pathogenicity of Pdd for fish, expression of the cytotoxins Dly, PhlyP, PhlyC and PlpV would be expected to be upregulated at 25˚C. Unexpectedly, expression of damselysin (Dly) (VDA_000159), PhlyP (VDA_000160) and PhlyC (VDA_002420) did not experience significant expression changes in growth at 25˚C compared to 15˚C (S3 Table), and PlpV (VDA_002242) was slightly downregulated at 25˚C (Table 2). This clearly indicates that levels of cytotoxin expression at 15˚C do not represent a limiting step that would prevent disease outbreaks in fish from occurring at low temperatures. These observations support that, as suggested above, the major influence of increased seawater temperatures in the onset of Pdd outbreaks in fish farms may be connected to the upregulation of other cellular processes (higher division rate, motility, chemotaxis, other plasmid-encoded putative virulence factors, etc) and not to a differential production of the four major cytotoxins.
This observation prompted us to analyze the RNAseq data to identify which are the most expressed genes in the genome of Pdd at each temperature of the study. Transcript abundance was quantified as Fragments Per Kilobase of transcript per Million mapped reads (FPKM), a method that allows the comparison of transcripts abundance among samples and conditions (S6 Table). Notably, the dly gene was the ninth most expressed gene at 15˚C (Table 3) with transcript abundance levels similar to genes of ribosomal proteins, which are among the most actively transcribed genes in fast-growing prokaryotic cells [92]. Two cold shock proteins and the NAD-dependent glyceraldehyde-3-phosphate dehydrogenase were included within the 10 most expressed genes at 15˚C. The 10 most expressed genes at 25˚C all corresponded to ribosomal protein genes (Table 4). LSU ribosomal protein L24p (L26e) (VDA_003450) was the most expressed gene under both conditions.
To illustrate the dominance of damselysin toxin transcripts, the FPKM values of dly and other virulence-related genes, as well as a selection of genes related with secretion systems and housekeeping cellular functions were compared (Fig 5). Although far from the top 10 most expressed genes, the FPKM values of the mRNA levels for the two pore-forming toxins PhlyP and PhlyC were also higher than those of housekeeping genes as gyrB, recA, mreB and ftsZ. Dly was the most highly expressed virulence factor at the two temperatures of the study, being particularly the case at 15˚C, and its transcript levels were almost 3 orders of magnitude higher than those of the PlpV phospholipase, which is considered to only have a minor contribution to virulence and cell toxicity [25]. The plasmid-encoded putative virulence factors OmpU, Vep07 and Vep20 showed transcript abundance levels largely inferior to Dly cytotoxin, again reinforcing the dominance of Dly as the top-expressed virulence factor in Pdd. This observation is in agreement with early studies which described highly virulent Pdd strains are producers of "large amounts of a cytolytic toxin in vitro" [93] and reinforces the major role of   Table) were obtained for a selection of virulence and regulatory genes, secretion system genes and housekeeping genes including the top-expressed gene encoding ribosomal protein L24p, and compared using a logarithmic scale. Vertical error bars represent standard deviation of biological triplicates.
https://doi.org/10.1371/journal.pone.0210118.g005 pPHDD1 as a virulence plasmid that contributed to the evolution of highly hemolytic and highly virulent lineages of Pdd.

Conclusions
Regulation of expression of virulence factors in response to temperature has been widely studied in pathogenic bacteria infecting homeotherms [94][95][96]. Meanwhile, much less is known about the temperature-dependent regulation of virulence factors in fish bacterial pathogens [59,97]. The aim of the current study was to investigate which genes are differentially regulated in Pdd in a heat stress condition (25˚C), relative to a colder condition (15˚C), considering that the majority of fish farm outbreaks occur during warm summer months at temperatures close to 25˚C.  within pPHDD1 related to iron acquisition (transferrin receptor), adhesion (OmpU), serum resistance (Vep07-like) and defence against competitors, were also upregulated at 25˚C (Fig 6). Overall, this study unveils a large set of previously overlooked genetic networks in Pdd and points at a number of cell functions related to virulence and acclimatization to changes in the environment.
Notably, damselysin toxin was one of the most highly expressed genes in the cell at the two assayed temperatures, with expression levels comparable to the most expressed genes encoding ribosomal proteins. This finding highlights the importance of this phospholipase D for the bacterium and suggests that this toxin might fulfil other biological roles in addition to virulence. To the best of our knowledge, this study represents the first transcriptome-based analysis of Pdd, and it has allowed us to identify a large set of gene functions that surely will constitute a foundation for future studies. Currently we are in the process of investigating the role of the top-differentially expressed genes identified in this study in the pathobiology of Pdd.
Supporting information S1