Long‐lasting effects of chronic exposure to chemical pollution on the hologenome of the Manila clam

Abstract Chronic exposure to pollutants affects natural populations, creating specific molecular and biochemical signatures. In the present study, we tested the hypothesis that chronic exposure to pollutants might have substantial effects on the Manila clam hologenome long after removal from contaminated sites. To reach this goal, a highly integrative approach was implemented, combining transcriptome, genetic and microbiota analyses with the evaluation of biochemical and histological profiles of the edible Manila clam Ruditapes philippinarum, as it was transplanted for 6 months from the polluted area of Porto Marghera (PM) to the clean area of Chioggia (Venice lagoon, Italy). One month post‐transplantation, PM clams showed several modifications to its resident microbiota, including an overrepresentation of the opportunistic pathogen Arcobacter spp. This may be related to the upregulation of several immune genes in the PM clams, potentially representing a host response to the increased abundance of deleterious bacteria. Six months after transplantation, PM clams demonstrated a lower ability to respond to environmental/physiological stressors related to the summer season, and the hepatopancreas‐associated microbiota still showed different compositions among PM and CH clams. This study confirms that different stressors have predictable effects in clams at different biological levels and demonstrates that chronic exposure to pollutants leads to long‐lasting effects on the animal hologenome. In addition, no genetic differentiation between samples from the two areas was detected, confirming that PM and CH clams belong to a single population. Overall, the obtained responses were largely reversible and potentially related to phenotypic plasticity rather than genetic adaptation. The results here presented will be functional for the assessment of the environmental risk imposed by chemicals on an economically important bivalve species.


| INTRODUC TI ON
Highly polluted marine areas provide the opportunity to study the response of local populations to environmental variables. In fact, chronic exposure to pollutants may exert strong constant selective pressure on natural populations of aquatic animals. To cope with such pressure, marine organisms may respond by significantly altering the transcriptome profile of different organs, especially those involved in detoxification (e.g. liver). In the case of acclimatization, when phenotypic plasticity is responsible for the organism response to environmental stressors, transcriptional changes can be entirely reversible, temporary and repeatable (Hodgins-Davis & Townsend, 2009;Lopez-Maury et al., 2008;Whitehead & Crawford, 2006;Whitehead et al., 2017). In the case of adaptation, selection on genetic polymorphisms with permanent alterations of the transcriptome may be observed in exposed animals (Bélanger-Deschênes et al., 2013;Bozinovic & Oleksiak, 2010;Nacci et al., 2016;Reitzel et al., 2014;Whitehead et al., 2010Whitehead et al., , 2012. A third level of response is observed when epigenetic modifications (i.e. DNA methylation and histone modifications), without alteration of the DNA sequence, occur (Aluru et al., 2011;Rey et al., 2016;Suarez-Ulloa et al., 2015).
Alterations of the epigenome might be heritable, which makes epigenetic response as intermediate between acclimatization and genetic adaptation. In addition to these three levels of genomic response to environmental variables (phenotypic plasticity, genetic adaptation and epigenetic modifications), a fourth genomic component, the microbiome, is becoming increasingly recognized. In the last two decades, the relationships between resident microbial organisms and their eukaryotic hosts have been extensively explored, demonstrating that the host microbiota plays a vital role in host development and health, as well as in its responses to environmental perturbations and anthropogenic-induced changes in the marine realm (Gootenberg & Turnbaugh, 2010;McFall-Ngai et al., 2013;Bourne et al., 2016;Apprill, 2017;Milan et al., 2018;Milan, Maroso, et al., 2019). Concerning bivalve species, recent studies demonstrated that bivalve microbiota is shaped by the environment and strongly suggested the existence of long-term association between bivalve hosts and specific bacterial taxa (Dupont et al., 2020;Milan, Maroso, et al., 2019;Milan, Smits, et al., 2019;Offret et al., 2020).
The relative role of these four dimensions of response to environmental factors largely remains an open question. Here, we present the case of Manila clam (Ruditapes philippinarum) inhabiting two different sites within the Venice lagoon. The Manila clam is a filter-feeding organism living in seafloor sediment that represents the most important species for commercial clam landings in Europe and one of the major aquacultured species in the world. Its remarkable ability to adapt to a wide range of environmental conditions makes it an excellent model species for the study of adaptation to environmental changes, as well as for biomonitoring impacted marine areas (e.g. Ji et al., 2015;Matozzo et al., 2010;Moschino et al., 2010;Sacchi et al., 2013;Wang, Li, et al., 2010;Wang, Zhou, et al., 2010). Native to the Asian coasts of the Indian and Pacific oceans, the Manila clam was introduced to Europe a few decades ago (Chiesa et al., 2017;Cordero et al., 2017;Flassch & Leborgne, 1992;Utting & Spencer, 1992).
Localized areas characterized by high levels of pollution have been observed within the Venice lagoon since the early 1970s, when the rapid industrialization process of the Porto Marghera (PM) area led to severe chemical contamination of soil, groundwater and sediments. Despite the progressive reduction of industrial activities, the PM area still shows high levels of pollutants, specifically dioxins, halogenated compounds and trace metals, as testified by the ban on shellfish harvesting enforced by local authorities (DGR3366/2004;MODUS 4;Milan et al., 2015Milan et al., , 2018. In the past few years, we studied how such a contrasting situation in the same water body affects natural populations of the Manila clam. In a previous study, we showed specific transcriptome and biochemical profiles in clams collected in PM that reflected molecular signatures of chemical stress (Milan et al., 2011(Milan et al., , 2015. Subsequently, we explored the role of clam microbiota in response to environmental variables, demonstrating the recurrent presence of putatively detoxifying bacterial taxa in PM clams and the potential for host-microbial synergistic detoxifying actions (Milan et al., 2018). Finally, in a short-term common garden experiment, we observed that clams from PM showed a significantly less pronounced response to experimental copper exposure than clams originating from the clean area of Chioggia (CH), despite having measured similar concentrations of this metal in both PM and CH (Matozzo et al., 2013;Milan et al., 2016). The latter study suggested that chronic exposure to pollution (i.e. dioxin and PCBs) may influence how animals respond to other environmental stressors. Building on all the evidence discussed thus far, we set up a field experiment to test the hypothesis that chronic exposure to chemical pollutants may have substantial effects on the Manila clam hologenome (microbiome and host transcriptome) long after removal from the contaminated site. To further implement a highly integrative biological approach and increase the ecological relevance of the investigated responses, we performed transcriptome and microbiome analyses, supplemented by an evaluation of biochemical functional for the assessment of the environmental risk imposed by chemicals on an economically important bivalve species.

K E Y W O R D S
ecotoxicology, hologenome, host-microbiota interactions, phenotypic plasticity, Ruditapes philippinarum and histological profiles. Overall, this study provides new evidence of long-lasting effects of chronic exposure to pollution in a species of particular commercial interest. The results obtained will be functional to the scientific community and regulatory agencies to assess and predict the environmental risks imposed by chemicals on natural and farmed populations.

| Sampling activities
Approximately 200 adult Manila clams were collected in late January 2016 (T0) in a polluted site in inner Porto Marghera (PM_T0) and in the farming area of Chioggia (CH_T0), which is characterized by the absence of detectable chemical contamination (Figure 1). An additional 500 Manila clams collected in PM at T0 were transplanted to the CH site and placed into cages (100 × 100 × 25 cm) made of plastic netting. To provide duplicate assays, two cages were utilized per sampling time. These cages were partially buried into and anchored to the sediments. Investigated CH clams were sampled in one of the most important Chioggia clam farming sites and placed in cages in the same manner as the transplanted PM clams.
A total of 200 PM transplanted clams and 200 native CH clams were collected after one month (early March-T1: PM_T1 and CH_ T1) and after 6 months (early August-T2: PM_T2 and CH_T2) posttransplantation (see details in Table S1). To increase homogeneity across samples, Manila clams of similar size (3.5-3.8 cm shell length) and weight (9-12 g) were included in the study. After each sampling (T0, T1 and T2), clams were transferred to the laboratory, and for each sample set, the digestive glands of five individuals grown in different cages (2 and 3 clams/each) were harvested for transcriptome and microbiome analyses. In addition, for each sampling time/ site, haemolymph and digestive glands were rapidly removed from 30 specimens, pooled in 10 samples (each consisting of tissues from three specimens), frozen in liquid nitrogen and maintained at −80°C for biochemical analyses. An aliquot of haemolymph was also immediately processed for lysosomal membrane stability, and another aliquot was fixed with Carnoy's solution (3:1 methanol, acetic acid) for the microscopic evaluation of micronuclei frequency. In addition, 10 F I G U R E 1 Map of the Venice lagoon indicating the Manila clam sampling sites: Marghera (PM) and Chioggia (CH). Some individuals collected in PM were transplanted to CH and collected at 1 and 6 months post-transplantation replicates, each composed of whole tissues of at least 15 organisms, were also harvested and stored at −20°C for chemical analyses. The geographical coordinates, number of individuals and experimental analyses performed for each sampling site are summarized in Table   S1 and Table 1.

| Environmental parameters and chemical analyses
Water temperature and salinity were measured on-site at each sampling time with a field thermometer and a hand-held refractometer. Additional information on environmental parameters (salinity, temperature, dissolved oxygen, chlorophyll) was collected through monitoring stations located in close proximity to the two sampling areas from the Water Regional Authority (publicly available; Table   S1), allowing for definition of the annual trends of different environmental variables in the two areas of the investigated samples.
Aliphatic hydrocarbons (C10-C40), polycyclic aromatic hydrocarbons (PAHs), halogenated persistent organic pollutants and trace metals (As, Cd, Cr, Cu, Hg, Ni, Pb, V, Zn) were analysed in the whole body mass of clams by conventional procedures, based on gas chromatography with a flame ionization detector or mass spectrometry, high-performance liquid chromatography with diode array (DAD) and fluorimetric detection and atomic absorption spectrophotometry techniques, as described and validated elsewhere (Regoli et al., 2019). Details on the analytical methods and procedures for quality assurance/quality control are given in File S1. Additional information about bioaccumulation and concentrations of chemicals in the sediment has been provided by the Regional Monitoring Plan MODUS 4.

| Gene expression analyses
Total RNA was extracted from the digestive glands of five individuals for each investigated group (CH_T0, PM_T0, CH_T1, PM_T1, CH_ T2 and PM_T2) using the RNeasy Mini Kit (Qiagen) according to the manufacturer's instructions. Each extraction cycle included a sterility control. The same extracted RNA was examined for both gene expression (RNA-Seq) and microbiota analyses (16S). cDNA libraries were constructed from five individuals for each sampling site using a SureSelect Strand-Specific mRNA Library (Agilent Technologies) according to the manufacturer's protocol (details are reported in File S1). Library pools were sequenced using a HiSeq 4000 sequencer (Illumina) with a 50-bp single-end approach, yielding a total of 765,204,350 reads (details are reported in Table S1, sequences available in NCBI SRA; BioProject PRJNA612420).
A de novo version of the R. philippinarum digestive gland transcriptome was assembled starting from a library of two individuals (collected in PM and CH) sequenced a HiSeq 2500 sequencer (Illumina) with a 100-bp paired-end approach (details are reported in File S1). A total of 53,476 transcript contigs were obtained (reported in File S2) and used as the reference transcriptome for RNA-Seq read mapping. Transcriptome annotation was performed by a Blastx similarity search on SwissProt (UniProt), Homo sapiens protein Ensembl database, Danio rerio protein Ensembl database and Crassostrea gigas protein Ensembl database (Evalue < 0.0001). Of 53,476 unique sequences, 32,711 (61%) showed at least one significant match. Details and the annotation of each contig are reported in Table S1. Adapter trimming and read mapping were carried out on CLC Genomics Workbench v.10.1.1 (CLCbio). Detailed information is reported in File S1 and Table S1. Raw gene counts obtained from mapping (see File S1) were used as input in EdgeR (Robinson et al., 2010) to analyse differential gene expression. Samples were grouped according to sampling area/time and transcripts showing less than two counts per million (CPM) in 4 out of 5 samples per group (in all investigated groups) were discarded. Samples were then normalized by EdgeR using the Trimmed Mean of M-values (TMM) method. After filtering, a total of 29,303 transcripts were kept for subsequent analysis.
A likelihood ratio test was carried out in EdgeR to assess differentially expressed genes (DEGs) between PM and CH at T0, T1 and T2, with a significant fold change (FC) threshold set to >1.5 and a false discovery rate (FDR) set to <0.05. To detect differences in the transcriptome profiles of PM and CH across the three time points, we used Next maSigPro (Nueda et al., 2014), a software program developed to analyse time series experiments (see details in File S1).
Gene set enrichment analysis (GSEA) was employed to identify enriched functional categories. To focus on biological processes (BP) and KEGG pathways that were consistently differentially regulated in PM and/or consistently correlated to chemical concentrations in our previous studies (Milan et al., , 2015(Milan et al., , 2016, a panel of 27 gene sets was evaluated, which included pathways and BP related to cell cycle, DNA repair, detoxification processes, protein turnover, response to external stimulus, energy metabolism, immune response and cell processes (full list of gene sets is reported in Table 2).
An additional gene set named the 'Porto Marghera Gene Set' consisted of the six transcripts that were differentially expressed in PM clams compared to every other sampling site and across different seasons in our previous study (Glutathione sulfotransferase theta-1, Sulfotransferase 1C4, AChE, Cytochrome P450 3A9, Steryl sulfatase and Peroxisomal 2,4-dienoyl-CoA reductase; Milan et al., 2013Milan et al., , 2015 (see details in File S1).

| Histological analyses
Tissue samples collected from 5 Manila clams for each site/sampling time were placed in 10% neutral-buffered formalin for at least 48 h, routinely processed and embedded in paraffin. Specimens were cut into 5μm-thick sections and stained with haematoxylin and eosin for standard light microscopic evaluation. To investigate tissues morphology and the possible presence of organ lesions, a complete histological examination was performed on dorso-ventral section of the middle of the body of clams (Howard et al., 2004). Each section included the visceral mass, gonad, gills and mantle.

| Microbiota characterization
The same extracted RNA for gene expression (RNA-Seq) was used for microbiota characterization. A single sterility control pool was created by mixing 1 μl of each sterility control, and the resulting pool was included in the microbiome analysis. The RNA integrity number (RIN) index was calculated for each sample using Agilent 2100 Expert software. To reduce experimental biases in RNA-Seq analysis due to poor RNA quality, only RNA samples with an RIN >7 were further processed.
For microbiota analyses, 1 μg of RNA was reverse transcribed to cDNA using the Superscript IV Kit (Invitrogen, Life Technologies).
cDNA was diluted to 0.2 ng/μl and amplified in a 50 μl reaction that included 5 μl of diluted DNA and 1.5 μl of both the reverse and forward primers (10 μM) that specifically target the V3-V4 gene region of the bacterial 16S rRNA, as described by Milan et al. (2018).
Libraries were then pooled together based on their concentrations, and the final pool was quantified using a Bioanalyzer 2100 (Agilent Technologies) and sequenced by BMR Genomics with an Illumina MiSeq (2x300). Microbiome sequencing generated approximately 6 million reads (passing filters), averaging approximately 100,000 reads per pool (sequences available in the NCBI Sequence Read Archive (SRA) https://www.ncbi.nlm.nih.gov/sra; SUB7159366).
Raw reads were first analysed using CLC Genomic WorkBench version 10.1.1 to select high-quality sequences. Reads with a Phred score lower than 20 were discarded (https://www.qiage nbioi nform atics.com/). Sequences were uploaded in QIIME 2 (Quantitative insights into microbial ecology; Bolyen et al., 2019) using DADA2 (Callahan et al., 2016) to filter adaptors, low-quality sequences and to merge forward and reverse reads of each fragment and obtain high-quality representative sequences. After the quality-filter step, read merging and removal of chimaeric fragments, a total of 1,528,388 reads were retained, represented by 1,802 features.
Sequence alignment was performed using MAFFT software (Katoh & Standley, 2013) and then classified using the Python library Scikit-Learn on QIIME2. Taxa assignment was carried out using the SILVA- and a two-way ANOVA was carried out to identify different taxa between experimental groups at OTUs level (p-value < 0.05). The resulting VCF (Variant Call Format) file was further filtered to retain variants present in at least 80% of samples using Bcftools v1.11, according to the following criteria for biallelic SNPs and insertions/ deletions (indels): min-alleles = 2, max-alleles = 2, type = snp/indel, min-af = 0.01, exclude-min-quality<20, exclude-max-missing>0.2.

| Genetic data analysis
Next, genotypes for SNPs and indels (in 0/1 format) were extracted from the filtered VCF file using the Genome Analysis ToolKit (GATK) v4.1.9.0, and genotype counts by population (n = 2; 15 individuals for each population) were used as input for the BayPass package v2.2, a population genomics software primarily aimed at identifying genetic markers subjected to selection and/or associated TA B L E 2 GSEA performed at each sampling time. A panel of 28 gene set were analysed. Columns T0, T1 and T2 report the upregulated site and corresponding FDR-q value for the significant gene sets identified for at least one sampling time. N.S.: not significant. Red and green colour indicates up-and down-regulated molecular pathways, respectively to population-specific covariates (Gautier, 2015). Candidate loci under selection that were identified by BayPass as being significantly contrasted (p-value < 0.001 and Bayes Factor > 20) between CH and PM were then functionally annotated using Annovar (Wang, Li, et al., 2010). For all SNP loci, pairwise (CH vs PM) FST statistics were estimated using the adegenet (v 2.1.3) and hierfstat (v 0.5-7) R packages.

| Environmental parameters and chemical analysis
Negligible differences in temperature, salinity and chlorophyll concentrations were recorded between PM and CH (  and T2, while 2 DEGs coding for putative Cathepsin L (CATL1) and

| Transcriptome analysis
Oncoprotein-induced transcript 3 protein (OIT3) were commonly found upregulated in PM clams at the first and the last sampling time.
To consider the contributions of all genes, a GSEA was carried out. The results obtained are summarized in Table 2. Gene sets involved in the cell cycle, RNA processing, ribosome biogenesis and energy metabolism were significantly enriched in CH clams compared to PM clams at T0, while the pattern was reversed after translocation (significant upregulation in PM samples at T1 and/or T2).
Gene sets involved in xenobiotic metabolism followed an opposite trend, resulting in upregulation in PM at T0 and one month after translocation. In addition to precompiled gene sets, we also tested a custom set composed of six transcripts that were consistently upregulated in PM according to our previous research, irrespective of either the sampling season or sampling year (Milan et al., , 2015. This analysis showed significant upregulation in PM only at

| Biochemical analyses
Variations in biomarker responses according to site and period are reported in Table 3

| Histology
A total of 5 Manila clam samples for each site were histologically examined. The cross section taken included the whole body of each individual. Normal microscopic structures were found, and no anatomical abnormalities or histological lesions were underlined.  indels were retained after filtering and further analysed. Genetic differentiation expressed as FST (estimated according to Weir and Cockerham (1984)) was −0.00037 (bootstrap confidence interval with p-value 0.99). Contrast analysis using BayPass showed that no indels were significantly contrasted between the two sites, while 17

| Microbiota characterization
SNPs were highly significantly contrasted (File S8). Annotation of the variants using ANNOVAR indicated that 10 SNPs were located in annotated regions, whereas the remaining 7 lacked annotation information (File S5). None of the 17 contrasting SNPs were located in the DEGs, which had been identified in the transcriptomic analysis.

| Chemical contamination and its effects on the clam hologenome
Levels of chemical pollutants were confirmed to be higher in clams living in PM than in CH, especially for trace metals. Detected concentrations were within the range of previously reported values for this area and typical of areas exposed to anthropogenic disturbance (Boscolo et al., 2007;Moschino et al., 2012). Concentrations of PCDD/Fs, PCB and HCB were not measured in this study, as such data were available from the Regional Monitoring Program (MODUS 4) for the same year and the same locations (PM and CH), further highlighting higher levels of these organic pollutants in PM clams (Table S1). Likewise, PCDD/Fs, PCB and HCB were not analysed at T1 and T2 because the content of these compounds in clam tissue rapidly decreases to 1/10 of the original level after just one week in uncontaminated water (Milan et al., 2016;Raccanelli et al., 2008).  (Mezzelani et al., 2021). Moreover, the inhibition of ACOX activity in native PM clams, which catalyses the first reaction of β-mitochondrial oxidation of fatty acids, may be related to changes in the metabolism of energy resources driven by several classes of environmental pollutants (Mezzelani et al., 2021). At the transcriptome level, most DEGs at T0 overlapped with those observed in previous studies (Milan et al., , 2015(Milan et al., , 2016. In analogy with biochemical and cellular biomarker results, functional analysis of DEGs, genes in clusters 1-4 and GSEA all showed a clear indication of upregulated genes in PM involved in xenobiotic metabolism, protein turnover, lysosome, immune response and tissue damage in PM clams. In addition to SULT, GST and CYP450 playing key role in detoxification and widely discussed in previous studies (Milan et al., , 2015(Milan et al., , 2016, the upregulation of AQP-8 (FC > 4) should be also highlighted as it was already found permanently overexpressed in PM under controlled conditions and differentially regulated compared to CH clams in response to heavy metal exposure study (Milan et al., 2016). This protein belongs to a family of integral membrane proteins functioning as water-selective channels and playing an important role in cellular osmotic balanced volume regulation (Agre et al., 2002;King et al., 2004). Recent studies demonstrated an in- oxidative phosphorylation and citrate cycle suggesting higher energy budget in CH clams compared to PM. Downregulation of these key biological process in PM clams could be related to reduced energy input from food because of reduced feeding activity, with the net effect of energy trade-offs that may affect population dynamics by altering growth performance, development, reproductive output and immune function (Akberali & Trueman, 1985;Sendra et al., 2021;Sussarellu et al., 2016;Widdows, 1985). In the future, experiments measuring ecophysiological parameters (e.g. growth, respiration) under controlled conditions might provide confirmation to this hypothesis.
Overall, clam-associated bacterial communities were similar to those found in a previous study (Milan et al., 2018), indicating a rather species-specific composition of microbial communities in the clam hepatopancreas. The large variability in bacterial diversity observed across different PM individuals may indicate less stable, potentially dysbiotic bacterial communities associated with clams challenged by high chemical contamination, although such a hypothesis would require a larger number of individual microbiome data to be tested. Overall, the evidence supports the hypothesis that Manila clams chronically exposed to chemical pollution in PM respond consistently across different years at all investigated biological levels (transcriptome, biochemical and microbiome).

| Same season, no chemical pollution: The clam hologenome one month post-transplantation
In accordance with the rapid decrease of PCDD/Fs, PCBs and HCB to 1/10 of the original level after just one week in uncontaminated water (Raccanelli et al., 2008)  represented taxa), a potentially opportunistic pathogen that is frequently associated with unhealthy animals (Fan et al., 2013;Lokmer & Wegner, 2015;Tanaka et al., 2004) and has already been described in our previous study (Milan et al., 2018). Several of the other overrepresented taxa in PM clams have been frequently associated with polluted environments, including Porto Marghera. These include Pseudomonas (8 OTUs overrepresented in PM at T1), Acinetobacter (3 OTUs overrepresented in PM) and Zoogloea (Gao et al., 2016;Milan et al., 2018;Mollaeia et al., 2010;Ouyang et al., 2016;Sohn et al., 2004;Wasi et al., 2013). Key functions of microbiota in the physiology and health of vertebrate species have been widely investigated (Ruby et al., 2004), as well as how genetics, diet and other external factors influence the gut microbial community. Knowledge of bivalve microbiota is still underdeveloped, with the exception of a few studies related to either public health or mortalities that involve bivalve species of commercial interest (Romalde et al., 2013;Li et al., 2018;King et al., 2019;Milan, Maroso, et al., 2019).
It has been recently reported that the microbiota of filter-feeding animals may be significantly modified through the acquisition of microbial species directly from the environment, as well as by seasonal environmental fluctuations in chemical-physical parameters that characterize marine and lagoon environments (Beleneva & Zhukova, 2009;Dubilier et al., 2008;Meisterhans et al., 2016;Milan et al., 2018;Tanaka et al., 2004). However, stable host-microbiota associations have also been reported (e.g. Baldi et al., 2013;Desriac et al., 2014;Freese & Schink, 2011;Harris, 1993;Moriarty, 1997;Wegner et al., 2013;Zurel et al., 2011). In this respect, it has been shown that the microbiota composition of the clam hepatopancreas is siteand season-specific (Milan et al., 2018), but consistent within sites and seasons even across different years (Milan, Maroso, et al., 2019;Milan, Smits, et al., 2019). One month in the same environment was not sufficient to equalize the PM clam microbiota to that of CH animals, although it likely altered the microbial community, leading to the emergence of opportunistic pathogens such as Arcobacter spp.
Such modifications to the resident microbiota may be related to the upregulation of several immune genes in the PM clams (T1), which could represent the host response to the increased abundance of deleterious bacteria. In particular, overexpression of BDEF, RUNX1 and FKBP1A (cluster 5) suggests specific responses to pathogens.
BDEF is a well-known AMP with remarkable microbicidal activity against bacteria in bivalve species (Zhao et al., 2010); RUNX represents a transcription factor involved in the immune response (Song et al., 2019;Zhang et al., 2014), while FKBP1A belongs to a family of proteins implicated in the response to infectious pathogen (Chen et al., 2011). Other DEGs that suggest a targeted response against pathogenic microorganisms are those encoding mannose receptors (MRC1 and MRC2), lactose-binding lectin l-2 (CLEC10A) and perlucinlike protein (PLCL).

| Different seasons and different stressors: The clam hologenome 6 months after transplantation
After 6 months of translocation, only Pb tissue levels were still significantly higher in PM clams. These clams also exhibited a significant upregulation of the putative heavy metal-binding protein (HIP) (FC > 90), which has a potential role in metal binding and detoxification, as a carrier of divalent cations in the plasma (Hattan et al., 2001;Yin et al., 2005). In addition to chemicals, seasonality exerts a profound effect on the Manila clam hepatopancreas transcriptome, its associated microbiome, along with a marked modulation of several biological parameters normally used as biomarkers (Milan et al., , 2018. Several results from the present study confirm this evidence. The summer season is expected to induce higher oxidative challenge in clams due to higher water temperatures, higher metabolic demand and lower oxygen availability. An additional challenge is represented by energy investment in reproduction, as the Venice lagoon clams reproduce between June and September. Evidence from oxidative stress markers appears to support this scenario. The effect of season is quite evident on the microbial community composition irrespective of the clam origin, as the first axis in the PCoA clearly separated summer samples from winter samples (Figure 3a). However, significant differences still existed between PM and CH after 6 months, which was potentially related to the strong resilience of the clam microbiota, but also to a different response to seasonality in the two groups of animals. As the host is able to modulate resident microbiota, persistent differences in its physiology may translate into differences in microbial communities. The transcriptome appeared to be markedly different during the summer in the CH and PM clams, and the time-course analysis clearly showed three clusters (7-9) where a large set of genes were not differentially expressed at T0 and T1 but  . All this evidence may suggest that PM clams have a lower ability to respond to environmental/physiological stressors (higher temperatures and reproduction) even after a relatively long acclimation to a pollutant-free area. Similar evidence has been observed for the response to copper exposure, which was less marked in PM clams after removal from a contaminated site (Milan et al., 2016), although the time frame (three weeks) was shorter than that reported here.
In coming years, it will be fundamental to perform further analyses aimed at confirming (or rejecting) the hypothesis that the observed temporal variation in CH and PM clam hologenome is mostly attributable to season effects.

| Genetic analysis
Analysis of anonymous genetic variants from 2bRAD sequencing data obtained in a previous study from 62 CH individuals and 56 PM individuals collected in 2009 and 2011 (Milan, Maroso, et al., 2019;Milan, Smits, et al., 2019)  It would also be interesting to further investigate the interactions between the host and associated microbiota, focusing on changes in the microbial community and immune response in the host or the effects of chemicals on host microbiota. Such interactions may have short-term consequences on the host but also long-term effects, as observed in this study.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
16S sequence data are deposited in the SRA database (SUB7159366).