Intestinal epigenotype of Atlantic salmon (Salmo salar) associates with tenacibaculosis and gut microbiota composition

It remains a challenge to obtain the desired phenotypic traits in aquacultural production of Atlantic salmon, and part of the challenge might come from the effect that host-associated microorganisms have on the fish phenotype


Introduction
It has become increasingly clear that host-associated microbes often play a central role in many aspects of the biology of their vertebrate hosts [1].Besides contributing new elements to many of the classic questions in ecology and evolution, manipulation of host-microbiota interactions is expected to improve crop and animal health in production systems, including aquaculture [2,3].However, we still lack a thorough understanding of which factors shape the microbiota, which is crucial for the successful exploitation of host-microbiota interactions in health and food sciences [2,4].For advancing this field, a promising start is to study species harbouring simple microbiota communities in controlled settings.
The Atlantic salmon (Salmo salar) is an important species in aquaculture, and as with many other piscivorous fish species, the diversity of their gut microbiota community is relatively low [5][6][7].Most bacteria in the salmon gut are considered transient bacteria, as their abundance in the gut microbiota follows their abundance in the surrounding environment [8,9].However, intracellular bacteria of the genus Mycoplasma are consistently found in the gut microbiota of wild and farmed salmon and appear to be resident [9].Recent studies have suggested a symbiotic relationship between salmonids and Mycoplasma spp.where a Mycoplasma-dominated microbiota is associated with increased host metabolism, micronutrient provisioning, and immune stimulation [10][11][12][13].
However, it remains unclear why the salmon gut microbiota is often dominated by one or a few mycoplasma species.
In salmon, the stability of the gut microbiota has been linked to intrinsic host factors such as the cortisol-related response to acute stress as well as diseases [14,15].Specifically, two independent studies have reported a similar displacement of Mycoplasma sp. by the opportunistic pathogen Aliivibrio sp. in the gut of Atlantic salmon that suffered from the ulcerous disease tenacibaculosis, compared to unaffected fish confined in the same tanks [16][17][18].The simple microbiota dynamics in Atlantic salmon and low environmental and genetic complexity in aquaculture production cohorts make such systems well-suited for studying host-microbiota interactions [1].A recent theoretical investigation found that the displacement of Mycoplasma sp. with Aliivibrio sp.fits well in a mathematical model where a disease alters the salmon's interaction with its microbiota [19].But it remains theoretical how intrinsic host factors associated with a disease can affect the microbiota at a molecular level.
Epigenetic analyses can provide mechanistic insights into complex responses modulated through the epigenome-transcriptome-proteome axis, and different epigenetic marks can be studied using various techniques [20][21][22].Methylation of the 5th carbon position in cytosine nucleotides (referred to as DNA methylation) is a common repressive epigenetic mechanism affecting gene expression levels in vertebrate species [23].Previous studies in salmon have reported changes in DNA methylation related to temperature stress [24], general stress [25], micronutrient supplementation [26], and the bacterial disease piscirickettsiosis [27,28].Among these studies, those that combined epigenomic and transcriptomic data found that the methylation level, especially around transcription start sites and transcription termination sites, was negatively correlated with gene expression levels [25,28].These findings suggest that DNA methylation links external stressors with intrinsic host responses, and that one may identify genes involved in complex responses by finding DNA methylation differences and placing them in a genetic context.However, the epigenotype is tissuespecific and since the connection is largely unexplored in the intestinal tissue, the potential role of host epigenetic factors on the microbiota remains theoretical.
The aim of this study was to examine DNA methylation in the distal gut tissue in a cohort of Atlantic salmon, that differed in being either visibly affected or unaffected by a spontaneous tenacibaculosis outbreak previously linked to a wider-scale gut microbiota displacement [18].These data offer the opportunity to study differences in the gut epigenotype as a potential mechanism to understand how host organisms can affect their intestinal microbiota.

Experimental design and sample acquisition
Bozzi et al. 2021 described an experimental cohort of salmon reared in a flow-through system at the LetSea land facility in Bjørn, Norway (lat: 66.080, long: 12.586) [18].Two weeks before a planned growth trial salmon were transferred from a large holding tank with salinity around 18 ppt and randomly distributed into 12 smaller experimental tanks equilibrated to 34 ppt with water from the fjord next to the research facility.Transfer to seawater is a critical phase in the anadromous lifecycle of salmon associated with increased mortality in aquacultural production [29].After the seawater transfer a subset of the fish developed tenacibaculosis, a skin ulcerative disease diagnosed by a bacteriological analysis to be caused by Tenacibaculum dicentrarchi in the investigated outbreak [18].Critical to our study, and like in previous reports of tenacibaculosis, not all fish in each tank developed visible symptoms from the infection [16,30].The difference in resiliency resulted in two groups of fish that had been reared in identical conditions but differed in disease status at the time of sampling; hereafter referred to as sick and healthy (Fig. 1A, Fig. 1B).
Using bacterial 16S rRNA gene metabarcoding it was discovered that sick fish had not only a lower condition factor, but also a different gut microbiota composition as sick fish had a high relative abundance of Aliivibrio sp. in contrast to a high relative abundance of Mycoplasma sp.

Fig. 1.
Experimental setup; A) After adjusting salinity from 24 ppt to 34 ppt a subset of individuals developed tenacibaculosis in an industrial trial (red-coloured fish) while others appeared healthy and unaffected by the disease (blue-coloured).B) The daily count of dead fish across 12 replicate tanks per day during the outbreak.The arrow indicates the sampling nine days after the first signs of disease.C) Ten sick fish and ten healthy fish with distinct microbiotas were chosen for epigenetic profiling (microbiome profiles are adapted from Bozzi et al. 2021 [18]). in healthy fish (Fig. 1C) [18].We elected to benefit from these samples and prior results, to explore whether differences in the host epigenotype was in any way related to these differences.The distal part of the salmon intestine is particularly immunologically active and was the area investigated in Bozzi et al. 2021 [31].Therefore, distal gut tissue samples from the ten sick fish with the highest relative abundance of Aliivibrio sp. and the ten healthy fish with the highest relative abundance of Mycoplasma sp. were chosen for Whole Genome Bisulfite Sequencing (WGBS).The 20 samples selected from this criterion showed no systematic link between the disease status and replicate tank, as there were five occurrences of sick and healthy fish sampled from the same replicate tank (Fig. 1C, Table S1).

DNA methylation sequencing
Genomic DNA was extracted from the distal gut tissue of the 20 samples (Fig. 1) that had been stored in RNAlater at − 20 • C since sampling.Each thawed sample was washed twice in 1 ml of PBS buffer before 20-25 mg was cut into small pieces using sterile scalpels.DNA was isolated with the DNeasy Blood & Tissue kit (Qiagen) using the protocol for animal tissues.After quantification of DNA concentration using a Qubit 3.0 (Thermo Fisher Scientific), extracted DNA was provided to Novogene's commercial service for WGBS on the Illumina NovaSeq 6000 S4 platform for 150 bp paired-end sequencing aiming for 60 GB sequencing per fish corresponding to a coverage of 20×.
In addition to WGBS, we utilised Nanopore Cas9 Targeted Sequencing (nCATS) to validate methylation levels in ten regions each with a size of approximately 10,000 bp [32].DNA was extracted with the Blood & Cell Culture DNA Midi Kit (Qiagen), specifically targeting high molecular weight DNA with the below modifications.The remnant tissue samples weighing 30-40 mg from two sick and two healthy samples were prepared using the washing protocol described earlier in this section.Washed tissue were added to 19 μl RNAase A (100 mg/ml), 9.5 ml G2 buffer and 250 μl Proteinase K followed by overnight incubation at 50 • C. Samples were then extracted following the manufacturer's protocol with the only exception being that the DNA was eluted in 150 μl 1× TE buffer (pH 8.0, Invitrogen) after 2 h of incubation at 55 • C.
Similarly, DNA concentration and purity were quantified, and the fragment length was examined using Genomic DNA Screen tape in TapeStation 2200 (Agilent).
We used CHOPCHOP to locate potential target sites with the PAM motif NGG suitable for CRISPR-Cas9 cleavage around the ten regions [33].For each region, two crRNAs were designed for the sense strand and two were designed for the antisense strand respectively upstream and downstream the target region.The 40 crRNAs specified in Table S2 were designed based on the genomic location of the protospacer, the estimated efficiency of the crRNA, and the self-complementarity and were acquired from Integrated DNA Technology (IDT).Genomic DNA from the four samples were subsequently constructed to an nCATS library using the Ligation sequencing gDNA -Cas9 enrichment (SQK-CS9109) sequencing kit and protocol (Oxford Nanopore Technologies) and third-party consumables including the 40 crRNAs (IDT).Each library was sequenced on an individual minION flow cell (FLO-MIN106, Oxford Nanopore Technologies) and the sequencing reads were basecalled and quality-filtered using Guppy v5.0.11 with default parameters [34].

Data processing and filtering
Methylation sequencing reads were quality filtered before attaining methylation information.The quality of the WGBS data was examined using fastQC v0.11.8 before and after removing adapters and lowquality reads with AdapterRemoval v2.2.4 [35,36].Read pairs passing quality filtering were aligned to the Atlantic salmon reference genome ICSASG_v2 assembly version GCF_000233375.1 [37] using Bismark v0.22.1 [38].Uniquely aligned reads were deduplicated and used for methylation calling with bismark methylation extractor.All reports were gathered using MultiQC v1.9 [39].The output files from Bismark were imported into R v4.0.2 and analysed using MethylKit v1.18 [40,41].Only methylation calls from cytosines in CpG context with a coverage of minimum five reads in all samples were retained for downstream analysis.This is lower than the default setting of 10× in methylKit [39], however, ensuring a strand-specific coverage of 10× across 20 samples with WGBS require an immense sequencing effort.Additionally, a genome-wide coverage as low as 5× has been proven efficient to detect differentially methylated regions in a downsampling experiment [42].Furthermore, methylation calls from CpG sites sequenced at a depth above the 99.9% quantile were removed, as these likely represent low-complexity regions that attract low-complexity reads from multiple genomic regions.
The Nanopore reads were aligned to the same version of the reference genome using minimap2 v2.6 and the methylation levels in the targeted regions were called using Nanopolish v0.13.2 [43,44].To ensure robustness while retaining as much information as possible, only methylation levels from individual CpG sites with a strand-specific coverage of at least 5× were used for analysis and comparison.For specific program settings please refer to https://github.com/SrenBlikdal/Tenacibaculosis.

Analysis of general methylation patterns
The genome-wide methylation levels were used to examine general methylation patterns and differences between sick and healthy fish.The mean methylation level around the 79,000 transcription start sites in the reference annotation file was found in 50 bp windows.We tested if the sick and healthy groups had different centroids or dispersions based on the methylation level of all retained sites, with a PERMANOVA test implemented in the R-package vegan v2.6-4 [45,46].Similarly, a Principal Component Analysis (PCA) was performed based on the full dataset, and the quantile of sites with the highest variance in methylation level was used for hierarchical clustering based on Euclidean distance using prcomp and methylKit [41,47].The results of the analyses were plotted using ggplot2 [48].

Differentially methylated cytosines and regions
Significant Differentially Methylated Cytosines (DMCs) between sick and healthy fish were identified with methylKit's logistic regression, and adjusted for multiple testing [41,49].To oblige the concern of test statistic inflation, the inflation factor (lambda) was estimated using the Rpackage QCEWAS [50,51].Differentially Methylated Regions (DMRs) i. e. small regions with several DMCs were identified using the R-package edmr v0.6.4.1 [52].To identify DMRs, a predefined maximum distance of 150 bp between neighbouring CpG sites within a DMR was utilised instead of the 50 bp estimated by edmr, because 50 bp led to many short DMRs near one another.The candidate regions identified by edmr were filtered and only regions that contained at least 10 DMCs were defined as DMRs.

Genomic location of DNA methylation differences
To characterise the methylation differences relative to the genetic landscape, the reference genome was divided into genomic features, following the notation of Saito et al. 2020 with minor modifications [26].The number of DMCs within each type of genomic feature was counted and compared to the distribution of all CpG sites covered in the analysis using odds ratio (OR) and a G-test using Genomic Ranges v1.48.0 and Desctools v0.99.47 [53,54].We focused our functional genomic analysis on DMRs since they are indicated to have a stronger effect on gene expression than single DMCs and are less prone to false positive findings [55].All genes and putative promoters determined as the upstream 6000 bp region of genes, that overlapped with DMRs were S.B. Hansen et al. extracted using the reference annotation file and custom-made bash scripts available on the project GitHub site.The list of candidate genes was checked for updates in gene names and gene symbols on the NCBI gene database [56] (accessed December 5th, 2022).When functional information was used from gene orthologs in distantly related species, the amino acid gene sequences were downloaded from NCBI and compared for sequence homology using EMBOSS needle on EMBL-EBI's online server [57].

WGBS data quality and summary statistics
The accuracy and completeness of the epigenotypic data impact the quality of the downstream analysis.After the removal of adapters and reads with low-quality base calls, a mean (±SD) of 230 (±35) million WGBS read pairs with an average length of 149 bp remained per sample.A total of 65.85% (±2.67%) of these reads mapped uniquely to the salmon genome, and a mean (±SD) of 134 (±15) million reads remained per sample after deduplication.This corresponded to a mean (±SD) genome-wide coverage of 13.2× (±1.5).One healthy sample (D7, Table S1) was removed from subsequent analyses because of a low methylation level exclusively on reverse reads, which indicated a technical issue in the preparation of that specific sample [36].In total, 16.6 million CpG sites met the coverage requirement in all remaining samples, which represented 17% of all CpG sites in the salmon genome available for downstream analysis.

Strong genome-wide DNA methylation associated with tenacibaculosis
DNA methylation levels depended on the genomic context and disease status.DNA was mainly methylated in the CpG nucleotide context, where the mean methylation level was 77.5%, as opposed to <1% in other sequence contexts.The methylation level was generally low near transcription start sites and increased gradually towards the genomewide mean levels in both directions (Fig. 2A).The hierarchical clustering indicates six out of nine healthy samples clustering together and similarly, in the PCA the healthy and sick samples are separated along PC1 (Fig. 2B and Fig. 2C).Interestingly, both analyses indicate that the same three samples from healthy fish exhibit methylation patterns resembling that of sick fish (Fig. S3).Despite the three outlier samples, the sick and healthy groups showed differences in genome-wide methylation patterns, which was further supported by a significant PERMANOVA test (p-value<0.001).

Tenacibaculosis is associated with genome-wide hypermethylation
With observable differences across millions of sites, we targeted the sites with the largest methylation difference between the groups.Of the 16.6 million sites tested, we identified 19,012 DMCs (FDR < 0.01), where 67% (12,692) of the DMCs were significantly more methylated in the sick fish (hypermethylated) and conversely 33% (6320) were significantly less methylated in sick fish (hypomethylated).A total of 1337 DMCs were located within 89 DMRs which also showed a tendency towards hypermethylation as 54 were hypermethylated and 35 were hypomethylated in sick fish.The mean length of the DMRs was 784 bp, and the mean CpG methylation difference between the sick and healthy fish was 23% (Fig. 3, Table S5).Together, this indicated the DNA methylation differences in the gut tissue were unevenly distributed as 7% of the DMCs were located in the 89 DMRs which in total covered <0.0025% of the genome.Additionally, both DMCs and DMRs were mainly hypermethylated in sick fish.

Validation of DMR methylation level using Nanopore sequencing
As experimental validation of the obtained methylation differences, ten regions containing eleven DMRs were sequenced to a mean coverage (±SD) of 4.7× (±4.0) in four samples using nCATS.The high standard deviation reflected a mean coverage of 1.3× in the two healthy samples  and 8.2× in the two sick samples.Furthermore, the coverage varied among the regions, and consequently, we were only able to estimate the methylation level in respectively eight and seven DMRs in the two sick samples, D1 and D9 (Table 1).Despite the limited number of sites, the mean methylation levels were very consistent when comparing the mean of the two methods across DMRs.Further, a comparison of the methylation estimates between all the 391 sites with at least 5× coverage in both methods revealed a Pearson correlation coefficient of 0.91 (Fig. S4).Together demonstrating that the methylation estimates were robust, even when samples were extracted and sequenced with different methods.

Methylation differences associated with specific genomic features and genes
To explore the potential genetic consequence of the DNA methylation differences, the salmon genome was divided into genomic features relative to annotated genes (Fig. 4A).The number of DMCs was enriched in specific genomic features, as significantly more (p-value <0.001) were located in distal promoters, medial promoters, and gene starts compared to the empirical distribution of all tested sites (Table 2).The DMRs were similarly connected to the genetic landscape, and 56 out of the 89 DMRs were found in putative promoter regions or within gene boundaries.Conversely, gene features of different genes can overlap, which means that one DMR can overlap several genes in our model.Therefore 68 different genes had putative promoter regions or gene boundaries overlapping DMRs (Fig. 4A, Fig. S5).

Genes differentially methylated in sick fish with microbiota displacement
Relating the differential methylation in the 68 genes to phenotypic consequences is not a simple task, which is further complicated by the fact that only a few genes in the salmon genome have completely known functions.The functional annotation of the salmon genome is constantly improving and as a direct consequence hereof, 16 of the 68 genes had HGNC-style gene symbols or aliases on NCBI database (Table 2), whereas the remaining had the less informative automatically generated gene names with the prefix LOC (Fig. S5) [58].The functions of these genes suggest that blood-related functions may be differential epigenetically regulated in sick fish, as the distal promoter of the heme transporter solute carrier family 48 member 1a (slc48a1a) and the gene end of erythropoietin receptor (epor) were both hypermethylated (Fig. 4B).Additionally, hypomethylation of the distal and medial promoter of the appetite-regulating gene thyrotropin-releasing hormone (trh) could indicate an altered feeding behaviour in sick fish, which may affect the microbiota (Fig. 4B).The list also included genes with functions related to the environment in the gut potentially affecting the microbiota.The start of LOC106590732 with alias ramp1 was hypermethylated and is suggested to encode an ortholog of the receptor activity-modifying protein 1 (ramp1), which in mice is involved in mucus secretion (Fig. 4B, Fig. S6) [59].Similarly, the distal promoter of protein kinase C delta a (prkcda) was hypermethylated, and in mice, prkcd has been shown to be crucial for macrophage killing of bacteria (Fig. 4B, Fig. S6) [60].

Table 1
Validation sequencing with nCATS; The number of sites in the DMR sequenced to at least 5× with both nCATS and WGBS (N).The mean of the methylation estimates with respectively WGBS and nCATS from sites within the DMR.The standard deviation of the difference in methylation calls between the methods for individual sites (SD).Only data from samples and DMRs with comparable sites are shown.

Discussion
This study finds profound DNA methylation differences in both individual sites and larger regions in the distal gut tissue of Atlantic salmon associated with tenacibaculosis and microbiota displacement.While stress and disease-related epigenetic changes have previously been reported in Atlantic salmon, this study is the first to characterise differences in the intestinal tissue serving as a barrier between the host and its gut microbiota.

Tenacibaculosis is linked to profound methylation differences in the gut
The separation of sick and healthy fish by their genome-wide methylation patterns indicates that tenacibaculosis is strongly associated with changes in the gut epigenotype.Overall >19,000 sites were significantly differentially methylated, and although an inflation factor of 1.26 indicates this number might be inflated, there was a consistent tendency towards hypermethylation in sick fish.Correspondingly, hypermethylation was also found to be the primary response to piscirickettsiosis in the liver and head kidney in Atlantic salmon, which could indicate hypermethylation and decreased gene expression is an energypreserving systemic disease response [28].
Our analysis in section 3.2 indicated three samples classified as healthy and with a Mycoplasma-dominated microbiota but with epigenotypes resembling sick fish.This observation is consistent with a model where the epigenotype responds to internal systemic disease symptoms before the development of external ulcers and gut microbiota displacement.Future studies with temporal sampling may elucidate if our hypothesised model of epigenetic changes affecting the microorganisms is more common than the alternative model of microbial modulation of the host epigenome [27,61].However, with a single time point sample, we can only speculate on the sequence of events.Future studies with temporal sampling may also investigate the role of standing epigenetic variation on disease susceptibility.Regardless the limitations of our experimental design, the association between disease and the epigenotype is compelling and offers novel perspectives on using epigenetic biomarkers for diagnostic purposes.

Methylation differences aggregate in specific regions
The DMCs had a convincing non-random genomic distribution specifically clustering in DMRs and around transcription start sites.Methylation levels in neighbouring CpG-sites often correlate, and although the reason is not fully understood, it adds more statistical and biological significance to DMRs than single DMCs [55].The majority of the DMRs were found near or within genes, which is confounded by high CpG levels around transcription start sites [62].However, the number of Fig. 4. DMRs associated with tenacibaculosis and microbiota displacement; A) The gene model used for defining genomic features, where the upstream sequence is divided into proximal promoter, medial promotor, and distal promoter being respectively the first 250 bp, the next 750 bp, and the following 5000 bp from the transcription start site.The genes were divided into gene start, gene end, and gene body being respectively the first 2500, the last 2500 bp, and the remaining part of a gene.For the non-overlapping model, conflicts were resolved using the hierarchy: Proximal promoter > Medial promoter > Distal promoter > Gene start > Gene end > Gene body > Intergenic regions.B) The mean methylation level of the N sites in five DMRs and their position relative to the specified genes.Individual samples are coloured based on the relative abundance of Mycoplasma sp.Created with BioRender.com.

Table 2
Genomic distribution of DMCs and DMRs; The number of DMCs located in each genomic feature in a non-overlapping gene model (Fig. 4).In parenthesis the odds ratio (OR) between the observed number of DMC and the expected number from the empirical distribution of all tested sites.The result of a g-test, testing if the OR is different from one, is indicated by (*) where (***) equals a p-value <0.001.The number of hypermethylated (↑) and hypomethylated (↓) DMRs located in each genomic feature in an overlapping model.Differentially methylated genes were the genes, including promoters overlapping DMRs with a gene symbol or alias in the NCBI database on December 5th, 2022.Different promoter types and Gene starts are concatenated in the last two rows, because many DMRs overlapped several of these features of the same genes.DMCs which is unaffected by this bias was also enriched around transcription start sites.Together it suggests that the methylation differences were predominantly found in the regions regulating gene expression [25,28].The non-random distribution of methylation differences also emphasises the potential of using targeted approaches such as nCATS to optimise the level of information obtained with a given sequencing effort.Here, the method was successfully used to validate the methylation levels in specific DMRs in two samples.The two samples with insufficient coverage put forward inconsistency as a challenge in our validation.However, when sufficient coverage was obtained, the methylation level from the nCATS approach was highly similar to the WGBS methylation estimates, even though the nCATS data were based on independent DNA extractions and a limited number of sites.Therefore, we argue that nCATS could offer a useful approach for future studies aiming to validate significant differences in an independent set of samples without the need for comprehensive WGBS.Further, nCATS also offers information on hitherto understudied epigenetic modifications such as 5-hydroxymethylcytosine and N-6 methyladenosine [63].

Differential regulation of epor and slc48a1a associated with tenacibaculosis ulcers
We identified interesting links between the epigenotype, specific genes, and their functions related to the phenotype.We note that these are exploratory rather than conclusive and vulnerable to confirmation bias.However, the skewed nature of DNA methylation complicates quantitative analyses of these data, limiting us to a more qualitative interpretation of the results [64].Hypermethylation of the gene end of epor encoding an erythropoietin receptor and the distal promoter of slc48a1a encoding a heme transporter could be connected to the ulcers of the sick fish.Both genes are associated with anaemia in fish, as epor has been found upregulated during acute anaemia in Atlantic salmon and the slc48a1a ortholog in zebrafish was found upregulated following hemolysis [65,66].However, the fact that epor and slc48a1a were both hypermethylated and therefore putatively downregulated in sick fish contradicts our expectations, since the sick fish had ulcers and were exposed to anaemia.On the contrary, their functions may not be as relevant in the gut compared to other tissues closer to the ulcers.

Indirect epigenome-microbiota interactions through altered feeding behaviour
Feeding behaviour may also be altered in the sick fish as the putative appetite-controlling gene trh was hypomethylated.Sick fish in Bozzi et al. 2021 were significantly smaller in size than the healthy ones, which seems generic to most diseases [18].Mechanistically it may be related to the regulation of anorexic and orexigenic genes as part of an inflammatory response [67].In mammals, trh signalling has long been recognised to suppress feeding activity and to be a key regulator of brain-gut interactions, also expressed in peripheral tissues [68][69][70].The role of trh in teleosts is however less understood, but it can be speculated that altered feeding behaviour has consequences for the gut microbiota [71].

Genes involved in host-microbiota interactions
We examined gut tissue because it is the adjacent host tissue to the gut microbiota and differential epigenetically regulated genes therefore may affect the microbiota.However, only few genes are known to directly alter the microbiota composition by well-described mechanisms.One intensively studied gene in this context is prkcda whose promoter in our study was hypermethylated in the sick fish and associated with a reduction of Mycoplasma sp. and an increase in Aliivibrio sp.In a mice model, prkcd has been shown critical for immunity against various pathogens including Listeria monocytogenes and Mycobacterium tuberculosis [60,72].The protective mechanism in mice occurs through macrophage killing and repression of prkcd leads to fewer macrophages and a higher load of pathogenic bacteria in a tuberculosis mice model [60].In zebrafish, the function of macrophages is dependent on the irf8 gene, and the lack of a functioning irf8 gene has a severe effect on the gut microbiota community [73,74].Hence, the hypermethylation of the distal promoter of prkcda could suggest less macrophage activity in the gut of sick fish.This may explain the higher relative abundance of the potential pathogen Aliivibrio sp. as a consequence of less efficient macrophage-based control of the intestinal microbiota.
We also identified hypermethylation at the very start of LOC106590732, a suggested ortholog of the mammalian ramp1 gene.In mice, ramp1 is involved in mucus secretion in intestinal goblet cells affecting the mucosal thickness, and disruption of ramp1 is linked to microbiota dysbiosis [59].While the effect of hypermethylation of LOC106590732 could be a thinner mucus layer in sick fish and potential dysbiosis, the function of ramp1 in fish is less established.In addition to the mentioned genes, the 68 genes associated with DMRs also included transmembrane proteins, transcription factors, etc. that may be important for host-microbiota interactions but with functions yet to be described (Table 2, Table S5).

Final remarks
It is worth stressing that transcriptional regulation is controlled by a broad number of mechanisms not included in our analyses, and DNA methylation alone is insufficient to elucidate the full epigenometranscriptome-proteome landscape.Technological and analytical development continues to actualise comprehensive epigenotyping of multiple markers in varying tissues at a large scale.In parallel, the resolution of the techniques moves from bulk-sample input towards single-cell sequencing, which allows for separating the bulk signal into each cell type's contribution [75].Utilising these developments in an efficient manner remains a challenge but offers unprecedented opportunities for analysing the full epigenome-transcriptome-proteome in the forthcoming years.
Our results suggest differences in DNA methylation associated with tenacibaculosis in the distal gut tissue, which may give functional insights to understand the gut microbiota displacement seemingly connected to the disease.We highlighted genes in the gut that showed significant differential methylation and had functions linked to hostmicrobiota interactions.We acknowledge the importance of confirming the biological significant association between specific genes and the microbiota.This can be done by using suitable data sets with DNA methylation and microbiota information for independent replication, or it could be done by also incorporating transcriptomic or proteomic data for a more mechanistic understanding.
In the concrete system of salmon aquaculture, it is already well established that tenacibaculosis should be avoided, and neither the associated microbiota displacement nor epigenetic differences change this.However, our study does put forward the potential of using DNA methylation data for disease monitoring and in the future potentially also for microbiota manipulation [76].This adds value to understanding interactions between the host epigenome and the microbiota and becomes increasingly interesting as the epigenotype can be manipulated by environmental factors and may be passed down from generation to generation.

Funding
Fig. 1.Experimental setup; A) After adjusting salinity from 24 ppt to 34 ppt a subset of individuals developed tenacibaculosis in an industrial trial (red-coloured fish) while others appeared healthy and unaffected by the disease (blue-coloured).B) The daily count of dead fish across 12 replicate tanks per day during the outbreak.The arrow indicates the sampling nine days after the first signs of disease.C) Ten sick fish and ten healthy fish with distinct microbiotas were chosen for epigenetic profiling (microbiome profiles are adapted from Bozzi et al. 2021[18]).Figure created with BioRender.com.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2 .
Fig. 2. Genome-wide CpG methylation associated with tenacibaculosis; A) Mean methylation level of CpG sites in 50 bp windows located from 6000 bp upstream to 4000 bp downstream of the 79,000 annotated transcription start sites.B) Unsupervised hierarchical clustering of samples based on the Euclidean distance of the DNA methylation level from the 25% of CpG sites with the highest variance.C) PCA of DNA methylation levels where individuals are plotted on the two first principal components and coloured by disease status.

Fig. 3 .
Fig. 3. Significantly differentially methylated sites and regions associated with tenacibaculosis; Manhattan plot showing the -log10(p-value) of DMCs plotted against their position in the genome.Hypermethylated and hypomethylated DMCs are plotted respectively above and below the dotted red lines indicating a q-value of 0.01.DMCs located in DMRs are plotted in solid circles whereas all other DMCs are plotted as transparent circles.The inflation factor of raw p-values was 1.26.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) This project was funded by The Danish National Research Foundation award to M.T.P.G. (CEH -DNRF143) and the European Unionʼs Horizon 2020 Innovation Action grant to M.T.P.G. and M.T.L. (Hol-oFood -817729).