Dual RNAseq highlights the kinetics of skin microbiome and fish host responsiveness to bacterial infection

doi:10.21203/rs.3.rs-83693/v1

Download PDF

Research

Dual RNAseq highlights the kinetics of skin microbiome and fish host responsiveness to bacterial infection

https://doi.org/10.21203/rs.3.rs-83693/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 07 May, 2021

Read the published version in Animal Microbiome →

Version 1

posted 05 Oct, 2020

You are reading this latest preprint version

Background: Tenacibaculum maritimum is a fish pathogen known for causing serious damage to a broad range of wild and farmed marine fish populations worldwide. The recently sequenced genome of T. maritimum strain NCIMB 2154 T provided unprecedented information on the possible molecular mechanisms involved in the virulence of this species. However, little is known about the dynamic of infection in vivo, and information is lacking on both the intrinsic host response (gene expression) and its associated microbiome community. Here, we applied complementary omic approaches, including dual RNAseq and 16S rRNA gene metabarcoding sequencing using Nanopore and short-read Illumina technologies to unravel the host–pathogen interplay in an experimental infection system using the tropical fish Platax orbicularis as model.

Results: We showed that the T. maritimum transcriptomic landscape during infection is characterised by an enhancement of antibiotic catalytic and glucan catalytic functions but a decrease in the sulfate assimilation process, compared with in vitro cultures. The fish host concurrently displays a large panel of immune effectors, notably involving innate response and triggering acute inflammatory response. In addition, our results suggest that fish activate an adaptive immune response visible through the stimulation of T-helper cells, Th17, with congruent reduction of Th2 and T-regulatory cells. Fish were, however, largely sensitive to infection, and less than 25% survived after 96 hpi. These surviving fish showed no evidence of stress (cortisol levels) or significant difference in microbiome diversity compared with controls at the same sampling time. The presence of Tenacibaculum in resistant fish skin and the total absence of any skin lesions suggest that these fish did not escape contact with the pathogen, but rather that they prevented pathogen entry. In these individuals, we detected up-regulation of specific immune-related genes differentiating resistant individuals from controls at 96 hpi, which suggests a possible genomic basis of resistance, although no genetic variation in coding regions were found.

Conclusions: Here we focus in detail on the interplay between common fish pathogens and host immune response during experimental infection. We further highlight key actors of defence response, pathogenicity and possible genomic bases of fish resistance to T. maritimum .

General Microbiology

Microbiome

Gene expression

16S rRNA

Nanopore

Tenacibaculum maritimum

Co-infection

Pathogens remain a significant threat to biodiversity, livestock farming and human health [1]. Host–pathogen interactions rely on a complex balance between host defences and pathogen virulence. Through constant selective pressure, pathogens evolve mechanisms to overcome the host's immune system. Reciprocally, the host adapts to counteract and limit pathogen virulence. Although changes in gene expression as a result of host–pathogen interactions appear to be common [2–4], the mechanisms involved often remain poorly understood. A more in-depth understanding of host–pathogen interactions could potentially improve our mechanistic understanding of pathogenicity and virulence, thereby defining novel preventive, therapeutic and vaccine targets [5].

Dual RNAseq sequencing fulfils the need to simultaneously assess the expression of both host and pathogen genes [6–8]. Studies applying this approach to fish bacterial infection systems have flourished recently and show promise for deciphering the complexity of host–pathogen interplay [9–11]. Yet, none of these studies have simultaneously explored dysbiosis and associated changes in microbiome communities. In nature, co-occurrence of multiple pathogen species (co-infection) is frequent. Species interactions can be neutral, antagonistic or facilitative and most often shape strain virulence plasticity, resulting in increased disease virulence [12–14]. Despite their commonness, remarkably few studies have explored such models, i.e. when a host interacts simultaneously with multiple pathogens during co-infection [15]. During tenacibaculosis outbreaks in Platax, T. maritimum burden is also commonly associated with other pathogen co-occurrences, namely Vibrio spp [16]. Nevertheless, such an approach is seriously impaired by the unbalanced representation of the sequences from each compartment, most often favouring the host compartment [8,17]. This bias can be minimised by specific library preparation (i.e. mRNA depletion), in silico normalisation procedures, and/or by investigating models in which the pathogen burden is high.

Tenacibaculum maritimum is a fish pathogen with a worldwide distribution, known for its lethal effects on a broad range of wild and farmed marine fish populations. Major efforts have been undertaken to lessen the impact of this pathogen and/or increase fish immune resistance [18]. The mechanisms of infection and fish response remain largely unknown, however, which has significantly held back aquaculture development. Nevertheless, recent sequencing of T. maritimum strain NCIMB 2154^T genome provided unprecedented information on the putative molecular mechanisms involved in virulence [19]. These authors note, for instance, that T. maritimum displays a large array of evolutionarily conserved stress resistance-related effectors, as well as an expanded capacity for iron mobilisation [19].

Tenacibaculum maritimum adheres and rapidly colonises mucosal surfaces [16,20]. Infected fish show multiplication of T. maritimum on their external tissues leading to severe skin lesions followed by rapid fish death [16]. Similarly to other infection systems, therefore, mucosal surfaces, here mainly skin mucus, are considered as the first barrier against pathogens [21]. This physical and chemical barrier constituted by mucus also includes host immune effectors (innate and adaptive) that orchestrate a complex interaction network with the commensal bacterial community [22,23]. The identification of these multi-specific interactions within the mucus brings to the fore the microbiome as the cornerstone of host–pathogen interactions [21,24]. Indeed, dysbiosis (i.e., the imbalance or alteration of the microbial ecosystem leading to a diseased status) is directly involved in the severity of a disease [25–27]. Recent studies on zebrafish (Danio rerio) raised in axenic conditions or in the presence of probiotic bacteria underlined the crucial role of the microbiota in the development of the immune system, mucosal homeostasis and resistance to stress and pathogens ([28], for review see [29]).

In French Polynesia, recurrent tenacibaculosis infections have been the main obstacle to sustainable local fish aquaculture. Indeed, T. maritimum affects the only the locally-farmed orbicular batfish (Platax orbicularis) leading to very high mortality rates shortly after transferring hatchery fingerlings to off-shore marine cages.

Here we combined dual RNAseq and 16S rRNA metabarcoding sequencing approaches to investigate the molecular responses of the host and microbiome community (gene expression and microbiome composition) simultaneously during T. maritimum infection and recovery phases, using orbicular batfish as a model. We also integrated comparison with in vitro liquid cultures of T. maritimum to refine our knowledge of virulence-related genes and explore genomic and genetic bases of resistance in P. orbicularis.

Animal husbandry

Platax fingerlings used in the bacterial challenge were 58 days old (days post-hatching; dph). Fish were obtained from a mass spawning of six females and eight males induced by desalinisation. Broodstock included wild individuals caught in French Polynesia that had been maintained at the Centre Ifremer du Pacifique (CIP) hatchery facility for seven years, under the supervision of the Direction des Ressources Marines. Eggs were randomly distributed into six black circular 210-L fiberglass tanks at 50 eggs.L^-1 to give an average density of 30 larvae.L^-1. Half of the larvae were then reared in conditions following the standard procedures used in the CIP facility (open water system, normal salinity around 36 psu), referred to as ‘Standard’ rearing conditions [18]. The remaining half of the larvae were reared in conditions supposedly optimal according to previous experiments and referred to as ‘Recirculated’ rearing conditions. Here, they were put in a recirculating system that was desalinated to 24 psu until day 34 whereupon salinity was progressively raised to normal (36 psu). Moreover, commercial clay (Clay Bacter ®) was added daily at a rate of 1g/d/m³ per percent of hourly water renewal from day 1 to day 19, the beginning of the live-prey weaning period. Probiotic Pseudoalteromonas piscicida B1 (local strain), produced in the CIP facilities (see Supplementary Methods), was also added daily to the fish feed (0.5mL/day/tank) and water (0.5 mL/day/tank) at a concentration of 10⁹cfu/mL of bacterial suspension from day 1 to day 57. At the end of the larval phase, at day 20, 700 fingerlings/tank were randomly selected. When these reached an average weight of 1 g, they were sorted to remove the largest and smallest extremes from the batches and then kept at a density of 200 ind./ tank (1g/L). It should be mentioned that around day 40 some of the fish, mainly those in the recirculating system, started to show a decrease in appetite and heavy mucus losses although there was no mortality. After an analysis of the water flora, it was shown that Vibrio harveyi was present but not T. maritimum.

Platax larvae were fed four times a day on live prey (Brachonius sp. and Artemia spp.) before being weaned from day 16 to 23. Fingerlings were then fed on commercial micro pellets ranging from 0.3 to 1 mm for the Micro-Gemma and Gemma ranges (Skretting, Stavanger, Norway) and from 1 to 1.3 mm for Ridley (Le Gouessant, Lamballe, France) according to the standard previously established [18]. Seawater supplied to both systems was pumped from the lagoon, filtered through a 300-µm sand filter and two 25- and 10-µm mesh filters and UV treated (300 mJ/cm²). The recirculating system included a 500-L biological filter to regulate levels of ammonia and nitrite. All tanks were supplied with saltwater held at 28.4 ± 0.3°C at a constant photoperiod (12L:12D) and oxygen saturation was maintained above 60% in the tanks with air distributed via air stones. Water renewal ranged from 36 to 360 L/h and new water input into the recirculating system was of 11 ± 1%. Levels of ammonia and nitrite were monitored once a week by spectrophotometry (HANNA Instruments ®) to assess biofilter performance. Temperature, salinity and dissolved oxygen were measured daily (YSI ®) and uneaten food and faecal material was removed once a day.

a)Bacterial challenge

We used strain TF4 for the experimental infection of 58-dph fingerlings. TFA4 was isolated from the skin of an infected Platax orbicularis in French Polynesia in 2013 and was shown to belong to Tenacibaculum maritimum by whole-genome sequencing, displaying an average nucleotide identity of 99.6% with the reference strain NCIMB 2154^T [30]. Strain TFA4 was cultivated in nutrient Zobell medium (4 g L^-1 peptone and 1 g L^-1 yeast extract Becton, Dickinson and Company, Sparks, MD in filtered and UV-treated seawater) under constant agitation (200 rpm) at 27°C for 48 h. On the infection day, 50 juvenile Platax orbicularis were transferred into 40-L tanks supplied with air and infected by the addition of 10 mL of a bacterial suspension of strain TFA4 to the tank water. Final bacterial concentration in the 40-L tanks, determined by the ‘plate-counting’ method, reached 4.10⁴ CFU.mL^-1. Moreover, for each rearing treatment, 16 to 17 fish were randomly selected from each tank to be transferred into one 40-L tank to form a mock-treated group, which we call the control (N = 50 fish) and 10 mL Zobell medium were added. After 2 hours of bathing, the fish were caught with a net, rinsed successively in two 40-L buckets filled with clean filtered UV-treated seawater, before being transferred back into their respective tanks. Twice a day, one third of the water in the tanks was replaced with filtered UV-treated seawater to maintain good water quality. This also made it possible to inactivate the T. maritimum in the sewage by bleach treatment. Dead animals were also collected and recorded at these times. After 115 hours post-infection (hpi), all living infected animals were considered as survivors and the challenge was ended. All the remaining fish were euthanised.

b) Animal sampling

We randomly sampled five individuals per tank at 58 dph (average weight 7.21g ± 0.28 se). Samplings in the experiment then consisted of five individuals per tank at 24 hpi and 96 hpi (Figure 1A). To increase our sampling dataset for controls we included the ‘standard’ control group (see previous section, one tank) together with the ‘recirculated’ (one tank) individuals. Consequently, our design consisted of four groups, namely control_24h(N = 10 individuals, replicate tanks)_,control_96h(N = 10 individuals, replicate tanks)_,infected_24h(N = 15 individuals, triplicate tanks) and resistant_96h(N = 15 individuals, triplicate tanks)_. For each sampling, at 24 hpi and 96 hpi, individuals were lethally anaesthetised using a benzocaine bath (150 mg.L^-1) and a lateral photograph was taken using a digital fixed camera (Leica Microsystems; Figure 1B and C). Microbiome and host sampling consisted in making gentle fish skin smears with sterile swabs that were directly placed in TRIZOL Reagent (Life Technologies) on ice to prevent RNA degradation. Swabs were disrupted using a mixer mill MM200 (Retsch) for 5 min at a frequency of 30 Hz and stocked at -80°C for later analysis. In parallel, water was also sampled from each tank but was not included in the analysis due to its very low DNA yield.

c) T. maritimum in vitro liquid culture sampling

The TFA4 strain was cultivated in 6 mL Zobell medium under constant agitation (200 rpm) at 27°C for 48 h, following exactly the same procedure and time of incubation as the culture used for bacterial challenge. Five culture replicates were performed. For each replicate, 4 mL at 10⁸ CFU/mL were centrifuged 5 min at 10,000 g at room temperature. Three inox beads and 2 mL TRIZOL (Life technologies) were quickly added to each bacterial pellet and cells were immediately disrupted using a mixer mill MM200 (Retsch) for 5 min at 30 Hz to prevent RNA degradation. RNA was extracted following manufacturer’s instructions, using a high salt precipitation procedure (0.8 M sodium citrate and 1.2 M NaCl per 1 mL of TRIZOL reagent used for the homogenisation) to reduce proteoglycan and polysaccharide contamination. Quantity, integrity and purity of total RNA were validated by both NanoDrop readings (NanoDrop Technologies Inc.) and on a Bioanalyzer 2100 system (Agilent Technologies). DNA contaminants were removed using a DNAse RNase-free kit (Ambion). A total of five RNA samples (1.048 ± 0.019 µg) were further dried in RNA-stable solution (Thermo Fisher Scientific) following manufacturer’s recommendations and shipped at room temperature to McGill sequencing platform services (Montreal, Canada). One library was removed prior to sequencing because it did not meet the minimum quality requirements.

d) Fish mortality and cortisol measurements

Mortality was recorded at 0, 19, 24, 43, 48, 67, 72, 91, 96 and 115 hpi. To estimate log-rank values, we used the non-parametric Kaplan-Meier approach implemented in the survival R package [31]. Differences in survival probability were considered significant when P < 0.05. We assessed stress levels in fish by measuring scale cortisol content [32]. The scales were collected from both sides of each individual, washed and then vortexed three times (2.5 min; 96% isopropanol) to remove external cortisol origination from the mucus. Residual solvent traces were evaporated under nitrogen flux and samples frozen at -80°C. To ensure the scales were dry, they were lyophilised for 12 hours and then ground to a powder using a ball mill (MM400, Retsch GmbH, Germany). Cortisol content was extracted from ~50 mg of dry scale powder by incubation in 1.5 mL methanol (MeOH) on a 30°C rocking shaker for 18 h. After centrifugation at 9500 g for 10 min, the supernatant was evaporated using a rotary evaporator and reconstituted with 0.2 mL of EIA buffer from a Cortisol assay kit (Neogen® Corporation Europe, Ayr, UK). Cortisol concentrations were determined in 50 µL of extracted cortisol using a competitive EIA kit (Neogen® Corporation Europe, Ayr, UK) according to a previously published protocol [33]. After validation of normality and homoscedasticity, differences among groups were tested using a two-way ANOVA followed by Tukey’s HSD post-hoc tests. Differences were considered significant when P < 0.05.

e) RNA and DNA extraction and sequencing

Dual RNAseq. Total RNA was extracted using the same procedure as described above. RNA was then dried in RNA-stable solution (Thermo Fisher Scientific) following manufacturer’s recommendations and shipped at room temperature to McGill sequencing platform services (Montreal, Canada). Ribo-Zero rRNA removal kit (Illumina, San 260 Diego, Ca, USA) was used to prepare rRNA-depleted mRNA libraries that were multiplexed (13–14 samples per lane) and sequenced on a HiSeq4000 100-bp paired-end (PE) sequencing device. Infected individuals 24 hpi were sequenced twice to insure sufficient coverage (Table S1).

Short-read 16S rRNA MiSeq microbiome sequencing. Total DNA was extracted from the same TRIZOL Reagent (Life Technologies) mix described above. DNA quantity/integrity and purity were validated using both a Nanodrop (NanoDrop Technologies Inc.) and a Bioanalyzer 2100 (Agilent Technologies). The V4 region was amplified by PCR using modified 515F/806rb primer constructs (515F: 5’-GTGYCAGCMGCCGCGGTAA-3’; 806rb: 5’- GGACTACNVGGGTWTCTAAT-3’) recommended for microbial survey [34]. Amplicon libraries were multiplexed and sequenced on a single lane of a MiSeq 250bp PE Illumina machine at Genome Québec McGill, Canada. Details of the sequencing statistics are given in Table S2.

Full 16S rRNA Nanopore sequencing. For a broad range amplification of the 16S rRNA gene, DNA was amplified using the 27F/1492R barcoded primer products (27F: 5’- AGAGTTTGATCMTGGCTCAG-3’; 1492R: 5’- TACGGYTACCTTGTTACGACTT-3’). In the PCR experiment, we included eight randomly selected individuals from the infected_24h group, two negative PCR controls (clean water) and one positive control (Acinetobacter DNA).

The PCR mixtures (25 μL final volume) contained 10 ng of total DNA template or 10 µL of water, with 0.4 μM final concentration of each primer, 3% DMSO and 1X Phusion Master Mix (Thermo Fisher Scientific, Waltham, MA, USA). PCR amplifications (98°C for 2 min; 30 cycles of 30 s at 98°C, 30 s at 55°C, 1 min at 72°C; and 72°C for 10 min) of all samples were carried out in triplicate in order to smooth out intra-sample variance. Triplicates of PCR products were pooled and purified by 1x AMPure XP beads (Beckmann Coulter Genomics) clean-up. Amplicon lengths were measured on an Agilent Bioanalyzer using the DNA High Sensitivity LabChip kit, then quantified with a Qubit Fluorometer.

An equimolar pool of purified PCR products (except for negative controls) was made and one sequencing library was finally prepared from 100 ng of the pool using the 1D Native barcoding genomic DNA protocol (with EXP-NBD103 and SQK-LSK108) for R7.9 flow cells run (FLO-MAP107) then sequenced on a MinION device. Details of the sequencing statistics are provided in the Supplementary Material.

f) Microbiome community analyses

Microbiome dynamics with the MiSeq short-reads dataset. Raw reads were filtered to remove the Illumina adaptors as well as for quality and length using Trimmomatic v.0.36 [35], with minimum length, trailing and leading quality parameters set at 100 bp, 20 and 20, respectively. Remaining reads were analysed with functions implemented on the QIIME2 platform v2019.10. Briefly, we used the DADA2 algorithm [36] to cluster sequences in amplicon sequence variants (ASVs). The resulting ASVs were mapped against GreenGenes v13.9 99% OTUs database [37]. We explored alpha-diversity [Shannon, Fisher and Faith’s phylogenetic diversity (PD) indexes] and beta-diversity (Bray-Curtis, unweighted and weighted Unifrac distances) using the phyloseq R package [38]. Dissimilarity between samples was assessed by principal coordinates analysis (PCoA). Differences in alpha-diversity were tested using pairwise Wilcoxon rank tests and were considered significant when P < 0.01. Differences in beta-diversity were tested using PERMANOVA (999 permutations) as implemented in the 'adonis' function of the vegan R package [39] and were considered significant when P < 0.01. We also searched for a ‘core’ microbiome of fish skin, considering those ASVs present in all the individuals across all treatments (infected, control and resistant) as belonging to this group. We finally searched for significant differences in specific ASV abundance across groups using Wald tests implemented in the DESeq2 R package [40]. We used the apeglm method for log2FC shrinkage to account for dispersion and variation of effect size across individuals and treatments, respectively [41]. Differences were considered significant when FDR < 0.01 and |FC| > 2.

Microbiome diversity analysis with the Nanopore dataset. Sequences were called during the MinION run with MinKnow software (v. 1.7.14). The demultiplexing and adaptor trimming were done with the porechop tool (https://github.com/rrwick/Porechop) using the option discard_middle. For each barcode, all Nanopore reads were mapped on the GreenGenes database (v.13.5, http://greengenes.lbl.gov) with minimap2 (v2.0-r191) with the ‘map-ont’ pre-set options [42]. All reference sequences of the GreenGenes database covered by more than 0.01% of all reads were kept for the next step. A second round of mapping (using the same parameters) was done on the selected references in order to aggregate reads potentially mis-assigned during the first round of mapping. SAMtools and BCFtools were used to reconstruct consensus sequences for each reference sequence covered by more than 10 Nanopore reads with the following programs and options: mpileup -B -a -Q 0 –u; bcftools call -c --ploidy 1;vcfutils.pl vcf2fastq. Individuals and consensus sequences were blasted (e-value < 10^-5) against the NCBI nt database(https://ftp.ncbi.nlm.nih.gov/blast/db; accessed October 2019).

g) Compartment-specific differential expression analyses

Read pre-processing. For each individual, raw reads were filtered using Trimmomatic v0.36 [43], with minimum length 60 bp, trailing 20 and leading 20. Filtered PE reads were mapped against a combined reference including the host’s transcriptome (See Supplementary Material Tables S1 and S3 for details of the transcriptome assembly) and the genomes of Alteromonas mediterranea strain: AltDE1 (Genbank accession: GCA_000310085.1), and Pseudoalteromonas phenolica strain: KCTC 12086 (Genbank accession: GCA_001444405.1), Tenacibaculum maritimum strain: NCIMB 2154T (Genbank accession: GCA_900119795.1), Sphingobium yanoikuyae strain ATCC 51230 (Genbank accession: GCA_000315525.1), Vibrio alginolyticus strain: ATCC 17749 (Genbank accession: GCA_000354175.2,) and Vibrio harveyi strain: ATCC 43516 (Genbank accession: GCA_001558435.2). To prevent multi-mapping biases, we used GSNAP v2017-03-17 [44] with minimum coverage set at 0.9, a maximum of five mismatches allowed, and removal of improperly paired and non-uniquely mapped reads (option ‘concordant_uniq’). Reads with low mapping quality (MAPQ) were removed using SAMtools v1.4.1 [45] with the minimum MAPQ threshold fixed at five. A matrix of raw counts was built using HTSeq-count v0.9.1 [46]. Transcripts from the host and bacterial species origin were then separated into different contingency tables using homemade scripts.

Host transcriptome analysis. Low coverage transcripts with count per million (CPM) < 1 in at least nine individuals were removed, resulting in a total of 22,390 transcripts. Similarly, transcript over-representation was assessed using ‘majSequences.R’ implemented in the SARTools suite [47]. We used distance-based redundant discriminant analysis (db-RDA) to document genetic variation among groups and correlation with group (infected or control), weight and time (24 and 96 hpi) as the explanatory variables. Briefly, we computed Euclidean distances and PCoA using the ‘daisy’ and ‘pcoa’ functions, respectively, implemented in the ape R package [48]. PCo factors (n = 6) were selected based on a broken-stick approach [49,50] and used to produce a db-RDA. Partial db-RDAs were used to assess the factor effect, checking for the other factor variables. We tested the significance of the models and individual factors using 999 permutations. Effects were considered significant when P < 0.01.

Differential expression was assessed using the DESeq2 R package [40], using pairwise comparisons with Wald tests. Logarithmic fold changes (logFC) were shrunk using the ‘apeglm’ method, implemented in the DESeq2 R package [40], to account for dispersion and effect size across individuals and treatments [41]. Differences were considered significant when FDR < 0.01 and FC > 2. Group comparisons included infected_24h vs control_24h and resistant_96h vs control_96h. Gene ontology (GO) enrichment was tested using GOAtools v0.6.5 [51] and the go-basic.obo database (release 2017-04-14) using Fisher’s test. Our background list included the ensemble of genes in the host transcriptome. Only GO terms with Bonferroni adjusted P < 0.01 and including at least three differentially expressed genes were considered. Significant GO enriched terms were used for semantic similarity-based clustering in REVIGO (http://revigo.irb.hr/).

Tenacibaculum maritimum gene expression in vitro or during infection. A validation step for searching for transcript over-representation was assessed using ‘majSequences.R’ implemented in the SARTools suite [47], similarly to the fish transcriptome. Most represented sequences were attributed to ssrA-coding genes, but represented less than 8% of the total library. We applied a similar shrinkage method and pairwise comparisons (infected vs in vitro), as for the host comparisons, but used more stringent thresholds as commonly observed in similar studies [8], and considered significant differences when FDR < 0.01 and FC > 4. Gene ontology (GO) enrichment was similar to the methods used for the host.

h) Species-specific weighted co-network gene expression analyses in the host

We built a signed weighted co-expression networks for the host compartment to cluster co-expressed genes and identify putative driver genes using the WGCNA R package [52]. Variation in normalised counts were previously controlled using the ‘vst’ method implemented in the DESeq2 R package [40].

We reduced the expression noise in the dataset by keeping only transcripts with minimum overall variance (> 5%). Briefly, we fixed a soft threshold power of 14 using the scale-free topology criterion to reach a model fit (|R|) of 0.80. The modules were defined using the ‘cutreeDynamic’ function (minimum 30 genes by module and default cutting-height = 0.99) based on the topological overlap matrix, a module eigengene distance threshold of 0.25 was used to merge highly similar modules. For each module, we defined the criteria for module membership (kME, correlation between module eigengene value and gene expression). We looked for significant correlation (Pearson’s correlation; P < 0.001) between modules and physiological data, including cortisol levels (pg.mg^-1 in scales), fish weight (g) and treatment (coded ‘1’ for control_24h, control_96h and resistant_96h and ‘2’ for infected₂₄). Gene ontology (GO) enrichment for each module was tested using same protocol and parameters as described above.

i) The genetic bases of fish resistance

We further explored the putative genetic variation between resistant and infected fish by focusing on resistant fish because of their established phenotype, (i.e. survivors with no signs of lesions after bacterial challenge). We followed GATK recommendations for SNP identification based on RNAseq data. Briefly, BAM files were pre-treated using the ‘CleanSam’ function, duplicates were picked out with the ‘MarkDuplicates’ function, and cigar string split with the ‘SplitNCigarReads’ function. All functions were implemented in GATK v4.0.3.0 software [53,54]. Final SNP calling was conducted with Freebayes v1.1.0 (https://github.com/ekg/freebayes) requiring minimum coverage of 15 and minimum mapping quality of 20, forcing ploidy at 2 and removing indels (‘—no-indels’) and complex polymorphisms (‘—no-complex’). The raw VCF file was filtered for minimum allele frequency (‘—min_maf=0.2’), minimum coverage (‘—minDP=20’) and using Vcftools v0.1.14 [55] to allow no missing data. We computed relatedness (‘—relatedness2’) within and among groups with Vcftools v0.1.14 [55]. We further used distance-based redundant discriminant analysis (db-RDA) to document genetic variation among groups and correlation, with cortisol, treatment and weight as the explanatory variables. Briefly, we computed Euclidean distances and PcoA using the ‘daisy’ and ‘pcoa’ functions, respectively, implemented in the ape R package [48]. Pco factors (n = 6) were selected based on a broken-stick approach [49,50] and used to produce a db-RDA. We tested the model significance using 999 permutations, effects were considered significant when P < 0.01.

We compared infected (with T. maritimum TF4 strain) with mock-treated P. orbicularis groups (hereafter referred to as ‘controls’) in experimental infection conditions, and individuals were sampled at 24 and 96 hpi. Surviving individuals in infected groups at 96 hpi are referred to here as ‘resistant’ fish. Sampling consisted in taking individual skin swabs for later DNA and RNA extraction. The DNA was used to assess microbiome communities (using 16S rRNA and Nanopore technologies), while the RNA served for gene expression analysis of the microbiome and the fish simultaneously (dual RNAseq). In addition, to detect putative virulence factors, we compared gene expression of T. maritimum during the in vivo infection with expression in the in vitro liquid culture. Fish condition, survival and cortisol levels were also monitored.

a) Fish weight, cortisol levels and mortality

Mortality rate in challenged fish reached 77.36 ± 18.35 (mean ± standard error; se) while no mortality was observed in the control group (Kaplan-Meier analysis, P < 0.001; Figure 2A). Mortality started at 24 hpi in the infected group and no novel mortality was observed, even after 72 hpi. Cortisol levels in fish scales vary significantly across groups (ANOVA; F = 9.46; P < 0.01; Figure 2B). Overall, cortisol levels were higher in the infected_24h group compared with all the other groups (Tukey’s HSD; P = 0.01). Cortisol levels in the control_24h group were also higher than in the control_96h and resistant_96h groups (Tukey’s HSD; t = -3.28; P = 0.01 and t = -3.42; P < 0.01, respectively). However, no difference was observed between control_96h and resistant_96hgroups (Tukey’s HSD; t = 0.12; P = 0.99).

b) Dynamics of host transcriptomic response to infection and search for genomic bases of resistance

Global mean unique mapping rate for skin smear samples reached 71.64 ± 2.99% relative to a combined reference for host and microbe compartments. Datasets were predominantly composed of host-origin sequences (mean 86.57 ± 13.48%), with the infected_24h group showing a significantly higher proportion of reads of non-host origin [mean 30.70 % ± 0.03 se] than other groups (Dunn’s test; Benjamini-Hochberg adj. P < 0.05; Figure S1). Details of the host transcriptome and individual mapping are provided in Tables S1 and S3.

Fish response to infection. Differential expression analyses revealed strong differences in host gene expression profiles between control_24h and infected_24h, with a total of 3,631 and 2,388 down- and up-regulated genes in infected_24h, respectively, compared with control_24h, (|FC| > 2; FDR < 0.01; Table S4). The infected_24hgroup responded to infection mainly by activating immune system response (Biological Process; BP), sterol biosynthetic process (BP), defence response (BP), inflammatory response (BP), regulation of biological quality (BP), lipid metabolic process (BP), iron ion homeostasis (BP), complement binding (Molecular Function; MF), heme binding (MF), oxidoreductase activity (MF), sulfur compound binding (MF), and (1->3)-beta-3-D-glucan binding (MF). A complete list of GO enrichment for each module is provided in Table S4.

We then used co-expression network analysis (WGCNA) to draw clusters of co-regulated genes associated with discrete (treatment) or continuous variables (weight and cortisol) and to identify putative hub genes. No gene module correlated significantly with fish mass, suggesting that mass had no significant effect on gene expression profiles. A total of three modules showed negative correlations (P < 0.01) with disease status (coded 1 for control_24h_,control_96h and resistant_96h groups and 2 for infected_24h), namely module_{turquoise-host} (r = -0.97, P < 0.001), module_black-host (r = -0.5, P < 0.001) and module_green-host (r = -0.49, P < 0.001; Figure 3). Inversely, two modules showed positive correlation with the treatment, namely module_blue-host(r = 0.97, P < 0.001) and module_pink-host(r = 0.48; P = 0.001). Almost all of these modules (with the exception of module_black-host) also correlated significantly with cortisol levels. The genes found up-regulated in control_24h clustered mostly in module_{turquoise-host} (n = 3,468; 95.5%), module_black-host (n = 72; 2.0%) and module_green-host (n = 39; 1.0%). Nearly all the genes found to be up-regulated in infected_24h clustered in module_blue-host(n = 2,352; 98.5%). Not surprisingly, the main driver genes (‘hub-genes’) include several transcriptional activators, such as, for module_{turquoise-host}, several Zinc finger proteins, Transcription factor GATA-3 (gata-3), Forkhead box protein O3 (foxpo3), activators of the autophagy pathways and the main drivers of naïve specific T-cell differentiation and activation [56], Runt-related transcription factor 2- (runt2) coding genes, involved in osteoblast differentiation, mineral- depositing cells and sialoproteins [57,58]. In module_blue-host, ‘hub-genes’ mainly reported actors of the innate immune system, inflammatory response, wound healing, oxidative and adhesion activity (Figure 3B).

Genomic bases of resistance. We found 38 DEGs between control_96h and resistant_96h (16 and 22 down and up-regulated in resistant_96h, respectively; |FC| > 2; FDR < 0.01). GO analyses show that adaptive immune response (BP) tends to be activated (uncorrected P < 0.001) in resistant_96h that includes genes related to pathogen recognition and immune response, such as C-type lectin domain family 4 member M, low affinity immunoglobulin gamma Fc receptor II-like and T-cell receptor beta variable 7-2-coding genes. Inversely, resistant_96hshow inactivation of ERK1 and ERK2 cascade (BP) regulation and repressors of the response to wounding (BP), regulation of the transforming growth factor-beta secretion (BP), alcohol biosynthetic process (BP) and regulation of interleukin-8 production (BP). However, GO enrichments were not considered significant under our threshold (Bonferroni adj. P > 0.05). A total of 27 (71.1%) out the 38 DEGs identified also showed different expression levels between control_24h and infected_24h. Among the 11 remaining genes (28.9%), we found Arf-GAP with dual PH domain-containing protein 1, C-type lectin domain family 4 member M, protein KIAA1324-like homolog, Ankyrin repeat and fibronectin type-III domain-containing protein 1 and T-cell receptor beta variable 7-2, up-regulated in resistant_96h group. Inversely, we found the Sal-like protein 1, Early growth response protein 1 and the Low affinity immunoglobulin gamma Fc receptor II-like down-regulated in resistant_96h. The complete list of DEGs and GO term enrichment is provided in Table S4.

Finally, we searched for genetic variation (SNPs) across resistant_96h and infected_24h(two established phenotypes) in order to identify putative variants associated with resistance capacities. We identified a subset of 13,448 filtered bi-allelic SNPs. Genetic variation analyses did not suggest any significant difference among groups (relatedness, Fst) and was not correlated with cortisol levels or fish mass of any of the groups based on the 13,448 markers (PERMANOVA; 1000 permutations; P = 0.18; Figure S2).

c) Microbiome flexibility and interactions among pathogen species and host response

Dynamics of microbiota communities on fish skin. The MiSeq sequencing strategy with amplification of the 16S rRNA V4 region of the 16S rRNA resulted a mean number of PE of 231,164.02 ± 36,542.99 sd out of which a mean of 82.47 ± 2.54 remained after filtering (Table S2). Species richness (Shannon) was lower in infected_24h than in other groups (Mann-Whitney-Wilcoxon; MWW; Holm adj. P < 0.001). Similarly, control_24h showed reduced species diversity values compared with control_96h (MWW; Holm adj. P < 0.001), but no difference was observed between resistant_96h and control_96h (MWW; Holm adj. P = 0.71; Figure S3A). A similar observation was made for Fisher’s alpha parameter, but this showed a significant difference between control_24h and resistant_96h (MWW; Holm adj. P = 0.05; Figure S3B). Bacterial communities varied significantly across groups (PERMANOVA; F = 30.93; P < 0.001); however, the beta dispersion also differed significantly (‘betadisper’, ANOVA, F = 30.93; P < 0.001). Specifically, based on Bray-Curtis distance values, communities were less variable in the infected_24h group compared with all other groups (Tukey’s HSD, [control_24h-infected_24h: 95% CI = 0.05–0.19; resistant_96h-infected_24h: 95% CI = 0.17–0.31; control_96h-infected_24h: 95% CI = 0.12–0.28]; P < 0.001, Figure S4). Communities were also more variable in resistant_96h compared with control_24h(Tukey’s HSD, 95% CI = 0.04–0.20; P = 0.001). The ASVs associated with Tenacibaculum (Flavobacteriales) are largely enriched in infected_24h compared with control_24h, but also significantly enriched in resistant_96h compared with control_96h (shrunken |log2FC| > 2; FDR < 0.01; Figure S4).

Long-read refinement of bacterial communities in infected fish. Results from Nanopore sequencing on full 16S rRNA sequences served to refine the taxonomy at the species levels, which might be limited with short-reads approaches. We amplified the full 16S rRNA of eight individuals randomly subsampled from the 24 hpi infected group, resulting in a mean number of SE reads of 60,019.62 ± 33,778.99 sd after pre-processing (the individual with the lowest coverage had a total of 29,520 sequences). Nanopore shows an over-dominance of T. maritimum in infected₂₄ (73.51 ± 4.89%) and confirmed the presence of other genera, including Vibrio, Polibacter, Alteromonas and Pseudoalteromonas (Figure S5). Among these genera, some species, such as V. harveyi [60], are known as potential fish pathogens, while others (Pseudoalteromonas) have been proposed as putative probiotics [61].\

d) Microbial compartment transcriptomic activity

Gene expression of Tenacibaculum maritimum during experimental infection vs in vitro. We examined the gene expression levels for T. maritimum in the fish during the peak of infection compared with in vitro to highlight putative genes associated with virulence (Figure 4). Mean total mapped reads against T. maritimum in vivo reached 5.65M ± 0.82 se (Figure S6), which is sufficient to conduct differential expression analysis [62].

We found a total of 72 and 142 DEGs up-regulated during experimental infection (in vivo) and in vitro, respectively (Shrunken |log2FC| > 2; FDR < 0.01). We only found the sulfate assimilation (BP) function enriched in vitro (Bonferroni adj. P < 0.1). Among the GO positively enriched during experimental infection (Bonferroni adj. P < 0.1) we found the glucan catabolic process (BP), external encapsulating structure part (cellular component, CC), pattern binding (MF) and the antibiotic catabolic process (BP). These processes include genes already highlighted in genome scale comparisons of Tenacibaculum species as possible virulence-related factors, such as the catalase/peroxidase katG or several Ton-B receptor dependent and Suc-C and SucD-like proteins [19] (Figure S7). In parallel, we also detected putative candidates involved in host membrane interactions and integrity such as Ulilysin, Streptopain, Pneumolysin toxin and peptidoglycan-associated factor lipoprotein or adhesines_.A complete list of GO is provided in Table S4.

Tenacibaculosis is a worldwide fish disease responsible for considerable farmed fish mortality events, but knowledge is lacking on the microbiome kinetics during infection and the concomitant host–pathogen interaction. The dual RNAseq method was chosen as provides unparalleled simultaneous data on the molecular features of the infection. It is particularly suitable for systems characterised by a massive pathogen burden with readily accessible material but for which cultures are not available [17,63]. Here, we adapted this approach for tenacibaculosis in P. orbicularis fish skin samples with the goal for comprehensively assessing the genomic basis and kinetics of infection as well as associated resistance mechanisms.

Infection modulates a wide range of innate and adaptive host immune effectors. Bath exposition of T. maritimum was highly efficient at inducing tenacibaculosis in juvenile orbicular batfish. The low survival rate and kinetics of infection support what is usually observed for other fish species [16,20,64]. Infected fish, sampled at the peak of infection (24 hpi), show large skin lesions characteristic of tenacibaculosis together with a high cortisol concentration in their scales [32]. Cortisol mediates changes in individual energy balance (e.g. mobilisation of energy stores, immunity, cognition, visual acuity or behaviour) [65,66]. This initial cascade of physiological and behavioural changes enables the organism to cope with acute stressors by mobilising adequate bodily functions, while concurrently inhibiting non-essential functions (e.g. reproduction, digestion) [67]. Here, increasing cortisol level reflects a local stress response to an unfavourable environment and is most likely involved in triggering the fish rapid immune response [68].

As expected, fish immune response, especially the innate immune system is strongly solicited at 24 hpi. Infected fish show activation of acute inflammatory response, mainly through driver genes, including interleukin-8 (IL-8) [69], but also activation of pathogen recognition receptors (PRRs), chemokines and antimicrobial-related humoral effectors. For instance, infection triggers co-expression of the cascade Toll-like receptor 5 (TLR5) and Myeloid differentiation primary response protein (MyD88), as previously reported in bony fish during bacterial infection [70]. However, the diversity of fish immune actors combined with the relatively limited knowledge we have on specific effector functions significantly hampered the comprehensive understanding of the mechanisms involved in our non-model species. For instance, in parallel to the TLR5, several other TLRs show reduced expression in infected fish, including TLR2 type-1, TLR-8 and non-mammalian (‘fish-specific’) TLR21. Despite previous effort towards assessing diversity of TLR sequences, protein-specific function remains poorly known in teleosts [71]. Similar observations have be made for the complement system, specifically complement C3, a key component of the immune system involved in ‘complementing’ antibodies for bacterial cell killing [72], for which several isoforms are reported in the Platax transcriptome. The different isoforms here have divergent patterns of expression (both up- and down-regulated in the infected_24hgroup), which support previously observed differences in target surface binding specificities [73].

Innate immune response is generally tightly linked to cellular homeostasis regulation and precedes adaptive immune response. The ability of the fish to maintain cellular homeostasis during infection is of primary importance when facing infection, and mechanisms include redox, biological quality control (autophagy) as well as ion level maintenance [74,75]; all of which were found to be affected in Platax. For instance, infected fish largely activate effectors of iron ion homeostasis. Iron, albeit largely present in the environment, is poorly accessible to organisms and iron sequestration and maintenance is a major mechanism developed by the host to limit pathogen growth as well as to regulate macrophage cytokine production [76]. In parallel, infected₂₄ individuals activated the (1->3)-beta-D-glucan binding process and thus contribute to the body of literature revealing receptor capacities in fish and pathway conservation (through PRRs, C-lectin and/or TLRs) across vertebrates and invertebrates [77,78]. Indeed, supplementation of β-glucan stimulates immune response in fish and increases resistance of the host to viruses and other pathogens (probably by reducing bacterial adhesion through lectin binding [79]); it therefore represents a promising immunostimulant for aquaculture [78,80]. Effects of β-glucan vary depending on species, exposure time, source of glucan, organs and markers monitored [77,81] and further studies will be needed to evaluate its potential at the production scale. Nevertheless, β-glucan is also relevant in bridging inflammatory response and activation / differentiation of T-cells in the adaptive immune response [82].

The adaptive immune response was also modulated at 24 hpi and its fine-tuned orchestration offers the opportunity to separate the preferential immune paths that can fight against T. maritimum infection. We identified several hallmarks of differentiated T-cells, indicative of the specialisation of the adaptive immune response to T. maritimum infection. Among the main driver genes of the response to infection in Platax, we noted a reduced expression of foxp3 and gata-3 in infected_24h. Both transcription factors are important regulators of the fate of Naïve CD4+ naïve T-cells, encouraging differentiation to T-regulatory (Treg) [83] and T-helper 2 (Th2) cells [56], respectively. Similarly, we showed reduced expression of T-bet transcription factors, a hallmark of Th1 cells [56]. Inversely, infected fish seem to activate Th17 cell differentiation, as suggested by simultaneous activation of the signal transducer and activator of transcription (STAT1-alpha/beta) and cytokine IL-17 [83]. Th17 cells are mainly dedicated to controlling bacterial and fungal entry [84]. In line with previous work [56,85], our results suggest a complex orchestration of T-cell differentiation via antigen communication and associated cytokine regulatory network in Platax during T. maritimum infection. However, we cannot rule out the possibility that changes in transcript abundance might also be indicative of cell migration. Complementary approaches, including cellular imaging [86] would clarify the presence of T-reg cells in fish and improve our knowledge of their regulation.

The genomic bases of resistance in P. orbicularis. Despite profound activation of the immune system response, most of the fish failed to contain the infection. At 96 hpi, less than 25% of the infected fish had survived the bacterial challenge. These surviving individuals did not display any skin lesions, suggesting that they resisted the T. maritimum penetration and/or limited initial bacterial adhesion. Considering the high bacterial concentration in the tanks during infection and the severity of the mortality event, it is very unlikely that resistant fish totally escaped contact with the pathogen. Indeed, T. maritimum were present in resistant_96h but not in control_96h individuals (see next section for details). Fish were thus able to maintain the integrity of their first barrier against pathogens, so we hypothesise that their differences in gene activity difference would reveal specific immune actors present and/or efficient enough (timing of regulation) to inhibit pathogen multiplication and entry in resistant fish [16]. We show that PRRs, specifically a C-type lectin receptor, was up-regulated in resistant_96h, together with a T-cell receptor and Low affinity immunoglobulin gamma Fc receptor II-like, while fibronectin-coding genes were down-regulated. These gene products are known for binding, agglutinating and neutralising bacteria [87], as well as triggering humoral immune response [88] or providing extracellular structure for pathogen adhesion through fibronectin-binding proteins [89,90]. Our experimental design did not permit us to tell whether there was a basal difference in expression in resistant_96h (genomic basis of resistance per se) or if the difference at 96 hpi was the result of a delayed adaptive immune response (timing of gene expression). Nonetheless, in catfish, lectin expression makes it possible to differentiate families resistant and susceptible to Flavobacterium colummnare, another gram-negative bacteria of the Flavobacteriaceae family [91]. Further longitudinal studies monitoring gene expression of resistant fish throughout the entire infection (including before the challenge) should prove useful in identifying resistance-specific responses to infection. Ideally, these studies should also simultaneously look at different immune-specific organ and tissue compartments and integrate genome-scale genetic variants (not limited to coding regions) to infer putative genetic bases of resistance.

Microbiome dynamics and host–pathogen communication. At 24 hpi, the microbiome community was dramatically affected by the over-dominance of T. maritimum. Abundance of T. maritimum was evident from metabarcoding data and contributed to significantly reducing species richness in infected_24h fish. We went further and compared expression of T. maritimum in vivo (during infection) compared to in vitro, with the hypothesis that key drivers of pathogenicity would be called upon to enable the bacteria to thrive and break host defence barriers. There are at least two major challenges that T. maritimum must to overcome to successfully infect the host: Pathogens need 1) to compete (at the intra and interspecific levels) for resources to metabolise from the local environment, and 2) to resist the host immune responses and stressful conditions. During infection, T. maritimum enhances its glucan catabolic activity. Although this might only reflect differences due to changes in the local environmental conditions (host mucus and skin) and/or resource availability, it also reveals some major mechanisms explaining the success of T. maritimum at growing on fish skin. Among the genes involved in glucan catabolism, we report several key components linking alternative food and mineral supplies and putative virulence-related functions, such as several genes involved in specifically degrading and uptaking sialoglycan and ions, mainly iron [19]. Sialidase activity explains why Capnocytophaga canimorsus bursts when in contact with host cells as this allows the pathogen to mobilise sugar directly from host phagocytes [92]. Similarly, the tonB-coding gene, regularly reported as a gene relevant for pathogenicity, confers virulence on Edwardsiella ctalurid by making it possible to maintain growth in an iron-depleted medium [93]. In parallel, several stress resistance-related genes were activated during experimental infection, all of which are also involved in the antibiotic catalytic functions. These genes include katA and katG, coding for two catalase-peroxidases involved in resistance to reactive oxygen species (ROS) by detoxifying exogenous H₂O₂ produced by host macrophages as a defence mechanism [94]. Obviously, the identification of virulence-related genes cannot be limited to those differentiating in vitro versus in vivo infectious status and other actors might be involved in making T. maritimum pathogenic. For instance, siderophore-coding genes are constitutively expressed in vitro or during experimental infection. These genes are a determining factor of host–pathogen and pathogen–pathogen interactions in the so-called ‘race for iron’ [95–97], which suggests that T. maritimum is highly efficient at mobilising iron independently of the local environment.

Finally, we mostly explored expression level variations in the light of an exclusive interplay between the host and T. maritimum, which might be effective considering the over-dominance of T. maritimum in fish mucus. However, most of the infection systems report several pathogen co-occurrences and the presence and/or activity of other opportunistic pathogens that might also play an important role in host fate [15]. In infected fish, we found that so-called ‘opportunistic’ bacteria were relatively largely represented. Opportunistic bacteria, which include V. harveyi, are known for their pathogenicity to fish. V. harveyi is a ubiquitous bacterium and one of the most common pathogens inducing major disease outbreaks in fish farming [60]. The enrichment of vibrio ASVs in resistant_96h mucosal communities and the absence of obvious associated physiological changes (cortisol, mortality, skin integrity) in this group, suggest that V. harveyi alone is not sufficient to induce mortality in Platax under our specific experimental conditions and associated bacterial burden.

Here we provide a comprehensive description of the interplay between host and T. maritimum under experimental infection conditions. Our results contribute to deciphering the complex orchestration of innate and immune responses of the host, but also suggest some promising avenues of research that could help to limit the impact of tenacibaculosis in fish farming. By taking an integrated ‘omic’ approach, we identified bacterial interactions as well as putative virulence-related genes in T. maritimum and candidate genes involved in fish resistance. Importantly, however, the detection of immune actors in fish and our comprehension of their regulation rely mainly on the quality of the annotation and the knowledge we have of their activity in other species, most of which are mammal model species. Consequently, further studies are now urgently needed to properly investigate and understand genetic and genomic bases of response to infection and possible resistance capacities in other non-model fish species.

Ethics approval. In vivo experiments complied with all the sections of deliberation n° 2001-16 APF from the Assembly of French Polynesia regarding domestic or wild animal welfare, issued in the Journal officiel de Polynésie française on 1^st February, 2001. In the absence of adhoc ethical committees, we followed European Commission DGXI [98] and ARRIVE [99] guidelines.

Consent for publication. Not applicable

Availability of data and material. Raw sequences have been deposited in the NCBI database under accession number (PRJNA656561). Codes are publicly available on Github repository https://github.com/jleluyer/metatranscriptomics_workflow

Competing interest. The authors declare no competing interest

Funding. The study was conducted as part of the Capamax project (Politique de site Ifremer) attributed to DS and JLL with financial support from the Aqua-Sana convention [Ifremer-Direction des Ressources Marines] for fish rearing and maintenance.

Authors’ contributions. DS, JLL and SP conceived the experiment. MC conducted fish husbandry and rearing. DS, MC, JLL and QS did the sampling. CB and CB made the DNA and RNA extractions. QS performed the cortisol level analyses. PA assembled the P. orbicularis transcriptome. JLL conducted the RNAseq and MiSeq analyses. JP did the Nanopore sequencing and JLL and QC the Nanopore analyses. All co-authors made substantial revisions to the manuscript.

Acknowledgements. We thank the Direction des Ressources Marines (DRM) for their help with fish rearing and especially Moana Maamaatuaiahutapu and Mike Teissier for the production of live preys. We thank Chedia Hamouna for her help with fish maintenance and monitoring during the experimental infection. We thank Céline Reisser for providing additional gonads tissues for the platax transcriptome assembly. We finally thank Eric Duchaud and Pierre Boudinot for providing valuable comments on previous version of the manuscript.

Daszak P. Emerging Infectious Diseases of Wildlife-- Threats to Biodiversity and Human Health. Science. 2000;287:443–9.
Celis JE, Kruhøffer M, Gromova I, Frederiksen C, Østergaard M, Thykjaer T, et al. Gene expression profiling: monitoring transcription and translation products using DNA microarrays and proteomics. FEBS Letters. 2000;480:2–16.
Kellam P. Post-genomic virology: the impact of bioinformatics, microarrays and proteomics on investigating host and pathogen interactions. Reviews in Medical Virology. 2001;11:313–29.
Kato-Maeda M, Gao Q, Small PM. Microarray analysis of pathogens and their interaction with hosts. Cellular Microbiology. 2001;7.
Casadevall A, Pirofski L. Host-Pathogen Interactions: Redefining the Basic Concepts of Virulence and Pathogenicity. Fischetti VA, editor. Infect Immun. 1999;67:3703–13.
Westermann AJ, Gorski SA, Vogel J. Dual RNA-seq of pathogen and host. Nature Reviews Microbiology. 2012;10:618–30.
Westermann AJ, Barquist L, Vogel J. Resolving host–pathogen interactions by dual RNA-seq. PLOS Pathogens. 2017;13:e1006033.
Rubio T, Oyanedel D, Labreuche Y, Toulza E, Luo X, Bruto M, et al. Species-specific mechanisms of cytotoxicity toward immune cells determine the successful outcome of Vibrio infections. Proc Natl Acad Sci USA. 2019;116:14238–47.
Huang L, Zhao L, Liu W, Xu X, Su Y, Qin Y, et al. Dual RNA-Seq Unveils Pseudomonas plecoglossicida htpG Gene Functions During Host-Pathogen Interactions With Epinephelus coioides. Front Immunol. 2019;10:984.
Zhang B, Zhuang Z, Wang X, Huang H, Fu Q, Yan Q. Dual RNA-Seq reveals the role of a transcriptional regulator gene in pathogen-host interactions between Pseudomonas plecoglossicida and Epinephelus coioides. Fish & Shellfish Immunology. 2019;87:778–87.
Valenzuela-Miranda D, Gallardo-Escárate C. Dual RNA-Seq Uncovers Metabolic Amino Acids Dependency of the Intracellular Bacterium Piscirickettsia salmonis Infecting Atlantic Salmon. Front Microbiol. 2018;9:2877.
Susi H, Barrès B, Vale PF, Laine A-L. Co-infection alters population dynamics of infectious disease. Nat Commun. 2015;6:5975.
Louhi K-R, Sundberg L-R, Jokela J, Karvonen A. Interactions among bacterial strains and fluke genotypes shape virulence of co-infection. Proc R Soc B. 2015;282:20152097.
Kinnula H, Mappes J, Sundberg L-R. Coinfection outcome in an opportunistic pathogen depends on the inter-strain interactions. BMC Evol Biol. 2017;17:77.
Kotob MH, Menanteau-Ledouble S, Kumar G, Abdelzaher M, El-Matbouli M. The impact of co-infections on fish: a review. Veterinary Research. 2016;47:98.
Avendaño-Herrera R, Toranzo A, Magariños B. Tenacibaculosis infection in marine fish caused by Tenacibaculum maritimum: a review. Diseases of Aquatic Organisms. 2006;71:255–66.
Rosani U, Varotto L, Domeneghetti S, Arcangeli G, Pallavicini A, Venier P. Dual analysis of host and pathogen transcriptomes in ostreid herpesvirus 1-positive Crassostrea gigas. Environ Microbiol. 2015;17:4200–12.
Reverter M, Saulnier D, David R, Bardon-Albaret A, Belliard C, Tapissier-Bontemps N, et al. Effects of local Polynesian plants and algae on growth and expression of two immune-related genes in orbicular batfish (Platax orbicularis). Fish & Shellfish Immunology. 2016;58:82–8.
Pérez-Pascual D, Lunazzi A, Magdelenat G, Rouy Z, Roulet A, Lopez-Roques C, et al. The Complete Genome Sequence of the Fish Pathogen Tenacibaculum maritimum Provides Insights into Virulence Mechanisms. Frontiers in Microbiology. 2017;8.
Rahman T, Suga K, Kanai K, Sugihara Y. Infection Kinetics of Tenacibaculum maritimum on the Abraded Skin of Japanese Flounder Paralichthys olivaceus. Fish Pathology. 2015;50:44–52.
Salinas I, Magadán S. Omics in fish mucosal immunity. Developmental & Comparative Immunology. 2017;75:99–108.
Parra D, Reyes-Lopez FE, Tort L. Mucosal Immunity and B Cells in Teleosts: Effect of Vaccination and Stress. Front Immunol. 2015;6.
Pérez T, Balcázar JL, Ruiz-Zarzuela I, Halaihel N, Vendrell D, de Blas I, et al. Host–microbiota interactions within the fish intestinal ecosystem. Mucosal Immunology. 2010;3:355.
Xu W, Yang L, Lee P, Huang WC, Nossa C, Ma Y, et al. Mini-review: perspective of the microbiome in the pathogenesis of urothelial carcinoma. Am J Clin Exp Urol. 2014;2:57–61.
Cadwell K. The virome in host health and disease. Immunity. 2015;42:805–13.
Cho I, Blaser MJ. The Human Microbiome: at the interface of health and disease. Nat Rev Genet. 2012;13:260–70.
Llewellyn MS, Leadbeater S, Garcia C, Sylvain F-E, Custodio M, Ang KP, et al. Parasitism perturbs the mucosal microbiome of Atlantic Salmon. Scientific Reports. 2017;7:srep43465.
Rawls JF, Samuel BS, Gordon JI. Gnotobiotic zebrafish reveal evolutionarily conserved responses to the gut microbiota. :6.
Kelly C, Salinas I. Under Pressure: Interactions between Commensal Microbiota and the Teleost Immune System. Front Immunol. 2017;8.
Bridel S, Bourgeon F, Marie A, Saulnier D, Pasek S, Nicolas P, et al. Genetic diversity and population structure of Tenacibaculum maritimum, a serious bacterial pathogen of marine fish: from genome comparisons to high throughput MALDI-TOF typing. 2020 [cited 2020 Aug 10]; Available from: https://pubag.nal.usda.gov/catalog/6938600
Therneau TM, Grambsch PM. Modeling Survival Data: Extending the Cox Model. Springer. New York; 2000.
Sadoul B, Geffroy B. Measuring cortisol, the major stress hormone in fishes. Journal of Fish Biology. 2019;
Carbajal A, Monclús L, Tallo‐Parra O, Sabes‐Alsina M, Vinyoles D, Lopez‐Bejar M. Cortisol detection in fish scales by enzyme immunoassay: Biochemical and methodological validation. Journal of Applied Ichthyology. 2018;34:967–70.
Walters W, Hyde ER, Berg-Lyons D, Ackermann G, Humphrey G, Parada A, et al. Improved Bacterial 16S rRNA Gene (V4 and V4-5) and Fungal Internal Transcribed Spacer Marker Gene Primers for Microbial Community Surveys. Bik H, editor. mSystems. 2016;1:e00009-15, /msys/1/1/e00009-15.atom.
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.
Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. DADA2: High resolution sample inference from Illumina amplicon data. Nat Methods. 2016;13:581–3.
DeSantis TZ, Hugenholtz P, Larsen N, Rojas M, Brodie EL, Keller K, et al. Greengenes, a Chimera-Checked 16S rRNA Gene Database and Workbench Compatible with ARB. Applied and Environmental Microbiology. 2006;72:5069–72.
McMurdie PJ, Holmes S. phyloseq: An R Package for Reproducible Interactive Analysis and Graphics of Microbiome Census Data. Watson M, editor. PLoS ONE. 2013;8:e61217.
Oksanen J, Blanchet FG, Kindt R, Legendre P, Minchin PR, Wagner RBO, et al. Vegan: community ecology package. R package. Version 2.0-3; 2012.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology. 2014;15:550.
Zhu A, Ibrahim JG, Love MI. Heavy-tailed prior distributions for sequence count data: removing the noise and preserving large differences. Stegle O, editor. Bioinformatics. 2019;35:2084–92.
Li H. Minimap and miniasm: fast mapping and de novo assembly for noisy long sequences. Bioinformatics. Oxford Academic; 2016;32:2103–10.
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;btu170.
Wu TD, Reeder J, Lawrence M, Becker G, Brauer MJ. GMAP and GSNAP for genomic sequence alignment: enhancements to speed, accuracy, and functionality. Statistical Genomics: Methods and Protocols. New York, NY; 2016;283–334.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics. 2009;25:2078–9.
Anders S, Pyl PT, Huber W. HTSeq—a Python framework to work with high-throughput sequencing data. Bioinformatics. 2015;31:166–9.
Varet H, Brillet-Guéguen L, Coppée J-Y, Dillies M-A. SARTools: A DESeq2- and EdgeR-Based R Pipeline for Comprehensive Differential Analysis of RNA-Seq Data. PLOS ONE. 2016;11:e0157022.
Paradis E, Claude J, Strimmer K. APE: Analyses of Phylogenetics and Evolution in R language. Bioinformatics. 2004;20:289–90.
Legendre P, Anderson MJ. Distance-Based Redundancy Analysis: Testing Multispecies Responses in Multifactorial Ecological Experiments. Ecological Monographs. 1999;69:1.
Legendre P, legendre L. Numerical Ecology, Volume 24 - 3rd Edition. Elsevier; 2012.
Klopfenstein DV, Zhang L, Pedersen BS, Ramírez F, Warwick Vesztrocy A, Naldi A, et al. GOATOOLS: A Python library for Gene Ontology analyses. Scientific Reports. 2018;8:10872.
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.
McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, et al. The Genome Analysis Toolkit: A MapReduce framework for analyzing next-generation DNA sequencing data. Genome Research. 2010;20:1297–303.
Auwera GAV der, Carneiro MO, Hartl C, Poplin R, Angel G del, Levy‐Moonshine A, et al. From FastQ data to high‐confidence variant calls: the genome analysis toolkit best practices pipeline. Current Protocols in Bioinformatics. 2013;43:11.10.1-11.10.33.
Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, et al. The variant call format and VCFtools. Bioinformatics. 2011;27:2156–8.
Mitra S, Alnabulsi A, Secombes CJ, Bird S. Identification and characterization of the transcription factors involved in T-cell development, t-bet, stat6 and foxp3, within the zebrafish, Danio rerio. FEBS Journal. 2010;277:128–47.
Brusgard JL, Passaniti A. RUNX2 Transcriptional Regulation in Development and Disease. In: Kumar R, editor. Nuclear Signaling Pathways and Targeting Transcription in Cancer [Internet]. New York, NY: Springer New York; 2014 [cited 2020 Sep 22]. p. 57–86. Available from: http://link.springer.com/10.1007/978-1-4614-8039-6_3
Le Luyer J, Deschamps M-H, Proulx E, Poirier Stewart N, Droit A, Sire J-Y, et al. RNA-Seq transcriptome analysis of pronounced biconcave vertebrae : a common abnormality in rainbow Trout (Oncorhynchus mykiss, Walbaum) fed a low-phosphorus diet. Journal of Next Generation Sequencing and Applications. 2015;2:1–13.
Kolde R, Kolde MR. Package ‘pheatmap.’ R Package. 2015;1.
Ina‐Salwany MY, Al‐saari N, Mohamad A, Mursidi F-A, Mohd‐Aris A, Amal MNA, et al. Vibriosis in Fish: A Review on Disease Development and Prevention. Journal of Aquatic Animal Health. 2019;31:3–22.
Offret C, Desriac F, Le Chevalier P, Mounier J, Jégou C, Fleury Y. Spotlight on Antimicrobial Metabolites from the Marine Bacteria Pseudoalteromonas: Chemodiversity and Ecological Significance. Marine Drugs. 2016;14:129.
Haas BJ, Chin M, Nusbaum C, Birren BW, Livny J. How deep is deep enough for RNA-Seq profiling of bacterial transcriptomes? BMC Genomics. 2012;13:734.
Montoya DJ, Andrade P, Silva BJA, Teles RMB, Ma F, Bryson B, et al. Dual RNA-Seq of Human Leprosy Lesions Identifies Bacterial Determinants Linked to Host Immune Response. Cell Reports. 2019;26:3574-3585.e3.
Avendaño‐Herrera R, Toranzo AE, Magariños B. A challenge model for Tenacibaculum maritimum infection in turbot, Scophthalmus maximus (L.). Journal of Fish Diseases. 2006;29:371–4.
Munck A, Guyre PM, Holbrook NJ. Physiological Functions of Glucocorticoids in Stress and Their Relation to Pharmacological Actions. Endocr Rev. Oxford Academic; 1984;5:25–44.
Johnson EO, Kamilaris TC, Chrousos GP, Gold PW. Mechanisms of stress: A dynamic overview of hormonal and behavioral homeostasis. Neuroscience & Biobehavioral Reviews. 1992;16:115–30.
Wingfield JC, Romero ML. Adrenocortical responses to stress and their modulation in free-Iiviing vertebrates. Comprehensive physiology. 2011;211–34.
Kulczykowska EZ. Stress response system in the fish skin - welfare measures revisited. Front Physiol. 2019;10:In press.
Li C, Yao C-L. Molecular and expression characterizations of interleukin-8 gene in large yellow croaker (Larimichthys crocea). Fish & Shellfish Immunology. 2013;34:799–809.
Basu M, Swain B, Maiti NK, Routray P, Samanta M. Inductive expression of toll-like receptor 5 (TLR5) and associated downstream signaling molecules following ligand exposure and bacterial infection in the Indian major carp, mrigal (Cirrhinus mrigala). Fish & Shellfish Immunology. 2012;32:121–31.
Palti Y. Toll-like receptors in bony fish: From genomics to function. Developmental & Comparative Immunology. 2011;35:1263–72.
Holland MCH, Lambris JD. The complement system in teleosts. Fish & Shellfish Immunology. 2002;12:399–420.
Sunyer JO, Tort L, Lambris JD. Diversity of the third form of complement, C3, in fish: functional characterization of five forms of C3 in the diploid fish Sparus aurata. Biochemical Journal. 1997;326:877–81.
Xia X, Wang X, Qin W, Jiang J, Cheng L. Emerging regulatory mechanisms and functions of autophagy in fish. Aquaculture. 2019;511:734212.
Ellis AE. Immunity to bacteria in fish. Fish & Shellfish Immunology. 1999;9:291–308.
Ganz T, Nemeth E. Iron homeostasis in host defence and inflammation. Nat Rev Immunol. 2015;15:500–10.
Dalmo RA, Bøgwald J. ß-glucans as conductors of immune symphonies. Fish & Shellfish Immunology. 2008;25:384–96.
Petit J, Bailey EC, Wheeler RT, de Oliveira CAF, Forlenza M, Wiegertjes GF. Studies Into β-Glucan Recognition in Fish Suggests a Key Role for the C-Type Lectin Pathway. Front Immunol. Frontiers; 2019;10.
Decostere A, Haesebrouck F, Van Driessche E, Charlier G, Ducatelle R. Characterization of the adhesion of Flavobacterium columnare (Flexibacter columnaris) to gill tissue. J Fish Diseases. 1999;22:465–74.
Ji L, Sun G, Li J, Wang Y, Du Y, Li X, et al. Effect of dietary β-glucan on growth, survival and regulation of immune processes in rainbow trout (Oncorhynchus mykiss) infected by Aeromonas salmonicida. Fish Shellfish Immunol. 2017;64:56–67.
Douxfils J, Fierro-Castro C, Mandiki SNM, Emile W, Tort L, Kestemont P. Dietary β-glucans differentially modulate immune and stress-related gene expression in lymphoid organs from healthy and Aeromonas hydrophila-infected rainbow trout (Oncorhynchus mykiss). Fish Shellfish Immunol. 2017;63:285–96.
Liu F, Weng D, Chen Y, Song L, Li C, Dong L, et al. Depletion of CD4+CD25+Foxp3+ regulatory T cells with anti-CD25 antibody may exacerbate the 1,3-β-glucan-induced lung inflammatory response in mice. Arch Toxicol. 2011;85:1383–94.
Secombes CJ, Wang T. The innate and adaptive immune system of fish. Infectious Disease in Aquaculture, Prevention and Control. Woodhead Publishing. Oxford, Cambridge, Philadelphia, New Delhi: Elsevier; 2012. p. 3–68.
Zhu J, Paul WE. Heterogeneity and plasticity of T helper cells. Cell Res. 2010;20:4–12.
Wang T, Secombes CJ. The cytokine networks of adaptive immunity in fish. Fish & Shellfish Immunology. 2013;35:1703–18.
Matheu MP, Othy S, Greenberg ML, Dong TX, Schuijs M, Deswarte K, et al. Imaging regulatory T cell dynamics and suppression of T cell priming mediated by CTLA4. Nat Commun. 2015;6:6219.
Weis WI, Taylor ME, Drickamer K. The C-type lectin superfamily in the immune system. Immunological Reviews. 1998;163:19–34.
Shishido SN, Varahan S, Yuan K, Li X, Fleming SD. Humoral innate immune response and disease. Clin Immunol. 2012;144:142–58.
Smani Y, McConnell MJ, Pachón J. Role of Fibronectin in the Adhesion of Acinetobacter baumannii to Host Cells. PLOS ONE. Public Library of Science; 2012;7:e33073.
Mosher D. Targeting the bacterial-host interaction. Virulence. 2012;3:349–50.
Beck BH, Farmer BD, Straus DL, Li C, Peatman E. Putative roles for a rhamnose binding lectin in Flavobacterium columnare pathogenesis in channel catfish Ictalurus punctatus. Fish & Shellfish Immunology. 2012;33:1008–15.
Mally M, Shin H, Paroz C, Landmann R, Cornelis GR. Capnocytophaga canimorsus: A Human Pathogen Feeding at the Surface of Epithelial Cells and Phagocytes. Cheung A, editor. PLoS Pathog. 2008;4:e1000164.
Abdelhamed H, Lawrence ML, Karsi A. The Role of TonB Gene in Edwardsiella ictaluri Virulence. Front Physiol. 2017;8:1066.
Hebrard M, Viala JPM, Meresse S, Barras F, Aussel L. Redundant Hydrogen Peroxide Scavengers Contribute to Salmonella Virulence and Oxidative Stress Resistance. Journal of Bacteriology. 2009;191:4605–14.
Cassat JE, Skaar EP. Iron in Infection and Immunity. Cell Host Microbe. 2013;13:509–19.
Choby JE, Skaar EP. Heme Synthesis and Acquisition in Bacterial Pathogens. J Mol Biol. 2016;428:3408–28.
Passalacqua KD, Charbonneau M-E, O’Riordan MXD. Bacterial metabolism shapes the host:pathogen interface. Microbiol Spectr. 2016;4.
Close B, Banister K, Baumans V, Bernoth E-M, Bromage N, Bunyan J, et al. Recommendations for euthanasia of experimental animals: Part 2. :32.
Kilkenny C, Browne W, Cuthill IC, Emerson M, Altman DG. Animal research: Reporting in vivo experiments: The ARRIVE guidelines: Animal research: reporting in vivo experiments the ARRIVE guidelines. British Journal of Pharmacology. 2010;160:1577–9.

LeLuyeretal.MicrobiomeSuppMat2clean.docx
Figure S1: Proportion of reads mapped to the microbial compartment. Read origin was dissociated in silico. Reads were considered to originate from the microbial compartment when no mapping was apparent in the fish transcriptome. Ctl-24h: control24h; Ctl-96h: control96h, Inf-24h: infected24h; Res-96h: resistant96h, groups. Different letters indicate significant differences, P < 0.05, Dunn’s test. Figure S2: Genetic variation and relatedness across established phenotypes. Principal component analysis for the total filtered dataset including 13,448 bi-allelic markers. Figure S3: Alpha-diversity estimates across groups. Alpha diversity was computed using A) Shannon (H’) and B) Fisher indexes. Ctl-24h: control24h; Ctl-96h: control96h, Inf-24h: infected24h; Res-96h: Resistant96h, groups. Different letters indicate significant differences, P < 0.05, Tukey’s HSD. Each dot represents a single individual. Figure S4: Taxa enrichment and beta-diversity dissimilarities across groups. A) ASV enrichment between infected24h (positive log2FC) and control24h (negative values). Colours represent different orders. The y-axis reports shrunken log2FC. Horizontal dash lines represent the log2FC threshold for significance (|log2FC| > 2). Nas represent missing taxonomic information for this ASV; B) ASV enrichment between resistant96h (positive log2FC) and control96h (negative values). Colours represent different orders. The x-axis represents associated genera; the y-axis reports shrunken log2FC. Nas represent missing taxonomic information for this ASV. C) PcoA based on Bray-Curtis distance values. Ellipses represent 95% CI intervals. Ctl-24h: control24h; Ctl-96h: control96h, Inf-24h: infected24h; Res-96h: Resistant96h, groups. Smaller points code for 24 hpi, larger points for 96 hpi groups. Figure S5: Relative abundance of bacterial genera in the infected24h group reported with Nanopore. 16S rRNA sequences and their relative abundance were generated from Nanopore reads (see methods). These sequences were blasted against the NCBI nt database (BLASTn; e-value<10-5). Dot size represents the number of unique sequences per genus. Figure S6: Total Illumina PE read counts for the most abundant species in the microbial compartment. The selected reference species were the most abundant species represented in the Nanopore 16S rRNA analysis (see methods). T.mar: T. maritimum; V.har: Vibrio harveyi; A.medit: Altermonas mediteranea; P.phe: Pseudoalteromonas phenolytica; V.algi: Vibrio alginolyticus; and S.yano: Sphingobium yanoikuyae. Each dot per species represents one individual. Figure S7: Plot of expression levels of in vitro and in vivo groups for all virulence-related genes previously identified. The gene list was obtained from whole-genome analysis in Tenacibaculum spp. [19]. Asterisks indicate genes with a significant difference between groups (Shrunken |log2FC| > 2; FDR < 0.01). Table S1: Individual mapping statistics against a combination of the P. orbicularis transcriptome and bacterial genomes. Transcriptome column refers to samples included in the metatranscriptome assembly. hpi = hours post-infection. PE = paired end. Mapping rate values are computed based on mapping against the combined reference transcriptome (host + microbiome). Table S2: Individual sequencing and mapping statistics for MiSeq Illumina reads. hpi = hours post-infection. PE = paired end. Table S3: P. orbicularis transcriptome statistics.
LeLuyeretalmicrobiome.TabS4.xlsx
Table S4: List of Fish and T. maritimum DEGs and GO term enrichments

Download PDF

Journal Publication

published 07 May, 2021

Read the published version in Animal Microbiome →

Version 1

posted 05 Oct, 2020

You are reading this latest preprint version

Dual RNAseq highlights the kinetics of skin microbiome and fish host responsiveness to bacterial infection

Status:

Journal Publication

Version 1

Abstract

Figures

Background

Methods

Results

Discussion

Conclusions

Declarations

References

Supplementary Files

Status:

Journal Publication

Version 1