Comparative analysis of the caecal tonsil transcriptome in two chicken lines experimentally infected with Salmonella Enteritidis

Managing Salmonella enterica Enteritidis (SE) carriage in chicken is necessary to ensure human food safety and enhance the economic, social and environmental sustainability of chicken breeding. Salmonella can contaminate poultry products, causing human foodborne disease and economic losses for farmers. Both genetic selection for a decreased carriage and gut microbiota modulation strategies could reduce Salmonella propagation in farms. Two-hundred and twenty animals from the White Leghorn inbred lines N and 61 were raised together on floor, infected by SE at 7 days of age, transferred into isolators to prevent oro-fecal recontamination and euthanized at 12 days post-infection. Caecal content DNA was used to measure individual Salmonella counts (ISC) by droplet digital PCR. A RNA sequencing approach was used to measure gene expression levels in caecal tonsils after infection of 48 chicks with low or high ISC. The analysis between lines identified 7516 differentially expressed genes (DEGs) corresponding to 62 enriched Gene Ontology (GO) Biological Processes (BP) terms. A comparison between low and high carriers allowed us to identify 97 DEGs and 23 enriched GO BP terms within line 61, and 1034 DEGs and 288 enriched GO BP terms within line N. Among these genes, we identified several candidate genes based on their putative functions, including FUT2 or MUC4, which could be involved in the control of SE infection, maybe through interactions with commensal bacteria. Altogether, we were able to identify several genes and pathways associated with differences in SE carriage level. These results are discussed in relation to individual caecal microbiota compositions, obtained for the same animals in a previous study, which may interact with host gene expression levels for the control of the caecal SE load.


Introduction
Salmonella is a zoonotic pathogen that can cause human foodborne disease. In 2019, more than 80,000 human salmonellosis cases were confirmed in Europe with 140 reported deaths [1]. Among these cases, more than 9000 were associated with 926 food-borne outbreaks (FBOs), with a large majority (72.4%) caused by the serovar S. Enterica Enteritidis (SE). Eggs produced by infected layer hens seem to be the major food vehicle, representing more than 37% of the FBOs. In parallel, in spite of strict hygiene control in farms, the systematic detection of Salmonella serovars, and the use of vaccination, the prevalence of Salmonella in laying hen flocks increased from 2.07% in 2014 to 3.44% in 2019 [1]. In chicken, the carriage is asymptomatic. The bacteria can persist a long time in the gut and can quickly spread within a contaminated farm via oro-fecal recontaminations between birds [2]. Understanding the impact of factors such as host genetics or gut microbiota on Salmonella carriage, and even more their combined impact, could lead to innovative strategies to reduce Salmonella transmission and ensure human food safety.
The caecal tonsil is a major barrier controlling the entry of bacteria in the organism [3,4] and is therefore a tissue particularly relevant for identifying host factors potentially involved in the control of SE. Several studies have been conducted on the caecal tonsil transcriptome. In particular, they have helped to identify biological processes associated with resistance to S. Enteritidis [5][6][7], S. Typhimurium [8,9] and S. Pullorum [10]. Nevertheless, in these studies, the impact of host genetics on gene expression was not examined. The expression of specific immune genes has been compared between the two experimental inbred chicken lines 6 1 and 15I, but not whole transcriptome [11]. The impact of host genetic variations was considered in a recent study of the caecal tissue transcriptome after Campylobacter colonisation. Comparisons between the experimental White Leghorn inbred chicken lines 6 1 and N led to the identification of a large number of differentially expressed genes, which may underlie variation in heritable resistance to the pathogen [12].
The host genetic background is an important factor for the outcome of Salmonella infection in chicken. A number of quantitative trait loci (QTL) and candidate genes associated with Salmonella resistance have been identified [13]. In the inbred lines N (resistant) and 6 1 (susceptible) in particular, several QTLs with low to moderate effects were identified [14][15][16]. However, no causal gene could be pinpointed due to the large size of the QTL genomic regions. More generally, only a few genes have been identified for their direct implication in the control of SE load in chicken, and knowledge is lacking about the mechanisms leading to genetic resistance. For the studies conducted on Salmonella carriage in the N and 6 1 lines, birds were reared together on floor after infection, thus allowing Salmonella oro-fecal recontamination between birds. In the present study, we used another infection model, making use of isolators. Previously tested on the experimental White Leghorn line PA12, this model showed a strong reduction of oro-fecal recontaminations, leading to much increased Salmonella individual variation among birds [2]. It is therefore an interesting model to identify birds with highly contrasted carriage levels, in order to facilitate the identification of host genes involved in these differences.
In the current study, we performed an analysis using information about caecal tonsil gene expression in two distinct genetic lines: the inbred chicken lines N and 6 1 , respectively resistant and susceptible to SE infection. The objectives of this study were to: i. identify differentially expressed genes between genetic lines in the caecal tonsils after SE infection, in order to identify potential pathways involved in the genetic resistance to SE; ii. identify genes and pathways associated with SE resistance within line (low vs high carriers);

Results
Two-hundred and forty animals from the two experimental White Leghorn inbred lines N and 6 1 were raised together on floor until 7 days of age. Then, chicks were challenged with Salmonella enterica Enteritidis (SE) LA5 by oral infection and separated into four isolators. Two independent replicates (n = 120) were conducted with a total number of 240 chicks. No clinical signs of disease were observed on the animals. Caecal contents and caecal tonsils were collected at 12 days post infection. The abundance of SE in caecal contents was measured by Droplet Digital PCR (ddPCR) and, as described previously, significant differences of Salmonella abundance were observed between lines for the two experiments [17]. The observed variability of carriage allowed us to identify extreme low and high carriers within each line, and 48 extreme animals were selected for the caecal tonsil RNA extraction, balancing the "experiment", "isolator" and "sex" factors. Different groups were defined as described in

Differentially expressed genes (DEGs) in caecal tonsils between lines
On average, more than 40M reads were sequenced for each of the 48 samples. After quality control, a total of 24,356 expressed genes were identified and used for the following analyses. Using all 48 samples, a principal component analysis (PCA) showed a distinct clustering between lines (Fig 2). Two ANOVA analyses on the PCA1 and PCA2 components showed that gene expression is significantly affected by line and sex (S1 File). A differential analysis with DESeq2 allowed the identification of 7,516 DEGs between the two lines (p-adj<0.05) among which 3,944 were up-and 3,572 were down-regulated (S1 Fig and S1 Table).

DEGs in caecal tonsils between low and high carriers within lines and experiments
Gene expression levels between low and high Salmonella carriers within each line and each experiment (L6L1/L6H1, L6L2/L6H2 and LNL1/LNH1, LNL2/LNH2) were compared. A PCA showed a clustering between the high and lower carriers within each of the N and 6 1 lines in experiment 1 (Fig 3; LNH1 vs LNL1 and L6H1 vs L6L1), but not in the equivalent groups of experiment 2. ANOVA analysis on the PCA1 and PCA2 components showed that gene expression is significantly (P<0.05) affected by low/high classes and by sex in experiment 1, but not in experiment 2 (S1 File).
A differential analysis with DESeq2 between low and high carriers within line 6 1 in experiment 1 (L6L1/L6H1) allowed the identification of 97 DEGs (p-adj < 0.05), among which 42 were up-and 55 were down-regulated (S2 Fig and S2 Table). A similar analysis performed between low and high carriers within line N in experiment 1 (LNL1/LNH1) allowed the identification of 1,034 DEGs (p-adj< 0.05) with 794 up-and 240 down-regulated genes (S3 Fig and S3 Table). Only 1 DEG was shared between these two comparisons (Fig 4). In experiment 2, no significant DEGs were found between low and high carriers regardless of the line (S4 and S5 Tables). The results are summarized in Fig 1. Merging both experiments in a single analysis, including a fixed effect for the experiment did not provide conclusive results.

DEGs common to intra-line and between-line analyses
When comparing DEGs identified between lines 6 1 and N, and within line 6 1 between low and high carriers, 38 genes were shared. Among these shared genes, 9 genes appeared to be When comparing DEGs identified between lines 6 1 and N and within line N between low and high carriers, 560 genes were shared. Among these shared genes, 58 appeared to be regu-

Functional enrichment analysis
Enrichment analyses were performed with topGO on DEGs between all animals from lines N and 6 1 and between low and high carriers within each line in experiment 1. The comparison between lines N and 6 1 led to the identification of 62 significantly (p-value < 0.05) enriched Biological Processes (BP) Gene Ontology (GO) terms (S6 Table). The comparison between low and high carriers of lines 6 1 (L6L1/L6H1) led to the identification of 23 BP GO terms (S7 Table). The comparison between low and high carriers of line N (LNH1/ LNL1) led to the identification of 288 BP GO terms (S8 Table).
groups infected with SE according to the line and experiment. Difference in mean between low and high carriers according to the line and experiment and ttest p-value. Difference in mean between lines and t-test p-value. Results of the DEG and BP GO term enrichment analyses. Results of the 16S analysis from the study [17].
https://doi.org/10.1371/journal.pone.0270012.g001 The results are summarized in Fig 1. Three BP GO terms were shared between the L6L1/ L6H1 and N/6 analyses (Table 5), fourteen between the LNL1/LNH1 and LN/L6 analyses ( Table 6) and four between the L6L1/L6H1 and LNL1/LNH1 analyses (Table 7). Interestingly, the 3 GO BP terms enriched and shared between L6L1/L6H1 and N/6 were related to the response to biotic stimulus and other organisms.

Discussion
We explored the caecal tonsil transcriptome after Salmonella Enteritidis infection by comparing samples from two genetic lines displaying contrasted levels of Salmonella carriage after infection: lines N and 6 1 . These inbred lines of chicken have been used in many infection studies, with results in apparent contradiction with ours. While N is more resistant to Salmonella Enteritidis in this study and in our previous experiments [14] in comparison to line 6 1 , it was more susceptible with Salmonella Typhimurium [18,19] or Campylobacter [12]. However, many factors may explain this difference: a potential genetic drift of the line (maintained for years at different experimental stations), environmental differences (different experimental farms), the infection model (the route of infection, the phenotype measured, the organ targeted, etc) [13]. However, to our opinion the most important factor is the pathogen itself: ST, Campylobacter or SE. Animals can display different resistance mechanisms to distinct pathogens, and even to distinct serotypes of the same pathogen. This was evidenced for instance in a  previous study comparing genetic locations of QTLs for resistance to SE and to ST: while two QTLs were probably common to both serotypes, several QTLs were specific of the serotype, thus revealing the probable existence of partially distinct mechanisms for those two serotypes of Salmonella [20].

A large difference of gene expression between lines
We showed a strong impact of the genetic background on gene expression in caecal tonsils, with 7,516 significant DEGs and 62 GO BP terms identified between the two lines. This strong difference in gene expression between lines is likely the result of genetic differences between the two lines. Some of these genes could explain the differences in susceptibility to Salmonella between lines. Thus, several genes display functions which may explain the higher resistance of line N, in which they were expressed to a greater extent. The latter genes code for the major histocompatibility complex I or II (MHCIBF2, MHCIYF5, MHCIIBLB1, MHCIIBLB2), antimicrobial peptides such as granzyme A and K, or the avian beta-defensins 10, 13 and 14. In the same way, the genes TLR 1,2,4,7,15 [21-24], NOS2, Gal 13, PSAP and IGL, which have already been associated with SE resistance in a genetic study in chicken [13], were more expressed in the resistant line N.
Interestingly, we showed previously that the caecal microbiota composition of these animals was highly different between individuals from lines N and 6 1 [17]. Do some of these DEGs in caecal tonsils indirectly impact Salmonella carriage through the modulation of the caecal microbiota composition? Conversely, do differences of microbiota composition indirectly impact Salmonella carriage through the modulation of gene expression in caecal tonsils? Further elements are needed to answer these questions.

A difference of gene expression between low and high carriers within line
The comparisons between low and high carriers within each line in experiment 1 (L6L1/L6H1 and LNL1/LNH1) revealed only one DEG and four BP GO terms in common, leading us to the conclusion that host pathways leading to a higher resistance could differ between lines. It

PLOS ONE
may be explained by their genetic differences. However, one of the four common BP GO terms may be directly related to the Salmonella infection process: the transmembrane receptor protein tyrosine kinase signalling pathway (GO:0007169), which is enriched in low carriers in both lines. Several studies have shown that the activation of protein tyrosine kinases and other specific transcription factors directly affects the innate immune response during SE infection [5,25]. Therefore, the transcriptional up-regulation of this pathway may be one of the mechanisms by which the animals from both lines control the infection. Interestingly, in spite of the smaller difference between low and high carriers observed in line N compared to line 6 1 , we identified many more DEGs and BP GO terms in line N, compared to line 6 1 (1034 DEGs and 288 BP GO terms in line N vs. 97 DEGs and 23 BP GO terms in line 6 1 ). Previously, we showed that the microbiota composition was significantly different between low and high carriers in line 6 1 but not between groups in line N [17]. One tentative hypothesis could be that the largest variability in Salmonella carriage in line 6 1 could be explained by differences in gut microbiota composition, as demonstrated with the PA12 chicken line [26]. The resistance against Salmonella in line 6 1 may thus be driven by mechanisms involving the intestinal microbiota, whereas in line N resistance mechanisms could be triggered by other factors, such as intra-line genetic variations.
Even in controlled environmental conditions with inbred lines, variability in the level of Salmonella carriage was identified within each line, as expected with the use of isolators preventing inter-individual recontamination [2]. This variability may derive from the residual genetic variability remaining within each line, or it may be caused by variations in the caecal microbiota composition, which is highly variable between lines and within line 6 1 [17,26]. However, this variability was much higher in experiment 1, allowing the identification of more contrasted low and high classes than in experiment 2. This higher variability, and thus greater contrast in host response to Salmonella infection between low and high carriers, is probably the reason why significant DEGs and BP GO terms could be identified in experiment 1 but not experiment 2. Experimental conditions were highly controlled and similar between experiments but some unchecked factors, such as the environmental microbial exposure at hatch, may differ and explain the observed differences between the two. Results in experiment 2 could be used as a "tendency validation" of the experiment 1. For example, in experiment 2, MUC4 was identified as more expressed in resistant animals, as in experiment 1 but with a higher pvalue, closer to the significance threshold (pvalue = 0.08 and padj = 0.9).

Identification of genes associated with the resistance to Salmonella infection
To highlight host genes that may be consistently associated with resistance to SE infection, we decided to focus on common DEGs between the intra-line and inter-line analyses. Indeed, genes more highly expressed in both the resistant line N and in low carriers in the susceptible line 6 1 , may be more likely to be associated with the resistance to SE infection. We further investigated if some of these genes have already been associated with the control of Salmonella infection in previous studies, supporting the reliability of our analysis. Four genes were identified as more highly expressed in both the resistant line N and in low carriers of the susceptible line 6 1 and are therefore associated with Salmonella resistance. Three of these genes have functions that are related to Salmonella: NOS2, DMXL1 and FUT2.  [38]. Finally, CNP (chitosan-nanoparticle) vaccination seems to increase NOS2 expression and protect against SE infection in chicken [39]. Second, DMXL1 is a gene involved in the phagosome acidification in macrophages. It seems to have an impact on innate immunity and macrophage bacteria killing after an activation by the TPL-2 kinase [40]. Finally, FUT2 is a gene coding for the α-1,2-fucosyltransferase enzyme involved in the glycosylation profile of the gastrointestinal tract. In mice and human, it has been shown that a "non-secretor" individual heterozygous for a loss-of-function mutation is more susceptible to chronic intestinal diseases such as Crohn's disease and to pathogen infection [41]. More specifically, it has been shown that "non secretor" mice (Fut -/-) show an increase of Salmonella Typhimurium in caecal tissue compared to wild-type mice [42,43]. Moreover, a FUT2 polymorphism was associated with the human faecal microbiota composition and diversity [44], which could explain host-microbe interactions and susceptibility to infection [45]. Finally, it has been shown that FUT2 was associated with the abundance of Christensenellaceae [46], a bacteria family we identified as associated with Salmonella Enteritidis resistance in the same animals [17]. Thus, FUT2 may be a gene indirectly associated with SE resistance through the modulation of the microbiota composition, modulating the abundance of competitive bacteria against Salmonella. Forty-two genes were identified as more expressed both in the resistant line N and in low carriers of the resistant line N and may therefore be associated with Salmonella resistance. Eight of these genes have functions of interest: MHCIBF2, CLDN7, SIRT5, ENTPD1, SYT7, SLC22A23, S100B and MUC4. The MHCIBF2 (Major histocompatibility complex class I antigen BF2) gene is the predominant ligand of cytotoxic T lymphocytes [47]. It has been shown that a particular MHC I haplotype may contribute to control the response to SE infection in chicken [48,49]. The CLDN7 (claudin 7) gene is involved in the formation of tight junctions between epithelial cells. It seems that a downregulation of CLDN7 by pathogen could facilitate translocation of invasive bacteria across the epithelium [50]. The SIRT5 (sirtuin 5) gene could have large impact on cellular homeostasis and is more expressed in colorectal cancer [51]. The ENTPD1 (ectonucleoside triphosphate diphosphohydrolase 1 or CD39) gene has an impact on inflammatory bowel disease (IBD) and regulation of pro-inflammatory responses and pathogen colonization [52][53][54], as does the SLC22A23 (Solute Carriers family) gene, which is associated with intestinal inflammation in human [55]. The SYT7 (synaptotagmin VII) gene is associated with the control of cytotoxic granule fusion in lymphocytes, and mice lacking syt7 have reduced ability to clear an infection [56]. The following two genes could be indirectly associated with SE resistance through the modulation of the microbiota composition. The S100B (S100 calcium binding protein B) gene codes for a signalling molecule which could be implicated in the communication mechanisms between microbiota and gut, and could explain differences between healthy and pathological microbiota in human [57]. Finally, it has been established that MUC4 (mucin 4, cell surface associated) and more generally MUC genes have an important role in preventing pathogens infections [58]. Many studies in chicken showed the association of the expression of mucin genes, especially MUC2, with SE infection [59][60][61]. Mucins favour the establishment and the maintenance of a commensal microbiota, and form a protection barrier against pathogens. In pigs, a decrease of MUC4 gene expression has been associated with ST infection [62], and genetic variants in this gene have been associated with the increase of gene expression relative to immune function and gut homoeostasis [63,64].

Identification of genes associated with the susceptibility to Salmonella infection
Five genes were identified as less expressed both in the resistant line N and in low carriers of the susceptible line 6 1 and may therefore be associated with Salmonella susceptibility. Among these genes, three belong to the cytochrome P450 family: CYP4B7, CYP2D6 and CYP2AC1. CYP enzymes play a key role in metabolic processes in the intestine as the metabolism of xenobiotic substances. Metabolism of CYP enzymes is closely connected with infection, inflammation and intestinal microbiota in human [65,66].
Fifteen genes were identified as less expressed both in the resistant line N and in low carriers of the resistant line N and may therefore be associated with Salmonella susceptibility. Four of these genes have immune functions: LECT2, DEFB4A, LOC121107850 and EXP. The LECT2 (leukocyte cell derived chemotaxin 2) gene has an antibacterial function, is expressed in chicken heterophils, and increases in abundance in macrophage after SE infection [53]. It is also more expressed in vaccinated chicken again SE. The DEFB4A (defensin beta 4A) gene has an important role in innate immunity in mucosal tissues through its antimicrobial activity against various microorganisms. It has been shown that beta-defensin plays a role in immunoprotection against Salmonella Enteritidis in in vitro embryonic chicken cell model [54] and in the development of innate immunity in gastrointestinal tract of newly hatched chicks [55]. The LOC121107850 (T-cell-interacting, activating receptor on myeloid cells protein 1-like) gene codes for an activating receptor on myeloid cells protein 1-like, and the EPX (eosinophil peroxidase) gene codes for an enzyme in myeloid, which has a function in bacterial destruction [67].
The increase in the expression of these genes is probably related to the high level of Salmonella infection. However, with this type of experimental design it is not possible to decipher whether host gene expression is responsible for high levels of Salmonella infection, or whether it is a consequence of the high infection.
Comparing the results obtained in the analyses between genetic lines on the one hand, and between high and low carriers within each line on the other hand, seems interesting to highlight genes having a meaningful impact in the response to SE infection. All of the genes discussed in these two last sections appear to be associated with SE infection and a modulation of the expression of these genes could prevent the colonization of the pathogen.

Identification of BP GO terms associated with the response to Salmonella infection
Some of the BP GO terms identified are implicated in immune response: notably immune system process (GO:0002376), defense response (GO:0006952) or cell surface receptor signalling pathway (GO:0007166). These BP GO terms have been already identified in another study working on the jejunum of a commercial genetic line and comparing high and low Salmonella infection using a kinome peptide array [68]. Thus, despite differences in experimental designs and methods used, similar results were found.
The three BP GO terms identified both between lines N and 6 1 and between low and high carriers in line 6 1 (LN/L6 and L6L1/L6H1) have interesting links to immunity: response to biotic stimulus, response to external biotic stimulus and response to other organisms ( Table 5). Seven of the fourteen BP GO terms identified both between lines N and 6 1 and between low and high carriers in line N (LN/L6 and LNL1/LNH1) are related to the regulation of the polymerisation or depolymerisation of proteins as actin (Table 6). It has been shown that a cytoskeletal actin rearrangement is induced by SE invasion in host cells [69]. Indeed, the actin cytoskeleton is targeted by Salmonella to promote its invasion, survival and growth in cells [70]. Genes implicated in these pathways may contribute to response to Salmonella infection and could be under genetic control.

Conclusion
The two experimental genetic lines N and 6 1 , displaying contrasted levels of Salmonella carriage after infection, showed a large difference of caecal tonsil gene expression associated with the outcome to SE infection. The comparison of resistant chicks (from line N or from low carriers within both lines) and susceptible chicks (from line 6 1 or from high carriers within both lines) allowed us to identify several genes and pathways associated with Salmonella resistance. Different mechanisms seem to be involved in the response to SE between these two experimental lines. A lower number of DEGs is associated with a larger inter-individual variability in line 6 1 compared to line N.

Experimental design
All animal procedures were authorised by the Ethic committee: APAFIS#5833-2016062416362298v3. Animals from the two experimental White Leghorn inbred lines N and 6 1 were provided by the experimental unit PEAT (Pole d'Expérimentation Avicole de Tours, Nouzilly, France). As described previously [17], animals were raised together on the floor until infection at the PFIE unit (Plateforme d'Infectiologie Expérimentale, INRA, Nouzilly, France) with free access to food and water. At 7 days of age, chicks were orally infected with Salmonella enterica Enteritidis (Strain 775 [LA5 Nal20Sm500], 5.10 4 cfu/0.2 mL/chick) and immediately separated into isolators. Four isolators were used for each experiment: two isolators for chicks from line N and two others for chicks from line 6 1 , with 30 birds per isolator. Caecal contents and caecal tonsils were collected at 12 days post infection (19 days of age) after the animal sacrifice and were immediately frozen in liquid nitrogen and stored at -80˚C until use. Two experiments were conducted with a total number of 240 chicks.

DNA extraction, Salmonella count by Droplet Digital PCR and choice of low and high carriers
As described previously [17], individual caecal DNA was extracted from an average of 200 mg of frozen caecal contents and DNA samples were stored at -20˚C. Individual abundances of Salmonella Enteritidis in caecal contents were obtained by Droplet Digital PCR (ddPCR) using the QX200 Droplet Digital PCR system (Bio-Rad) at the @bridge platform (INRAE, Jouy-en-Josas, France). The ddPCR method has been proven to reliably quantify the amount of salmonella spp. in a sample [71]. We targeted and amplified a region of the InvA gene specific to SE [72]. Data were analysed with a log transformation of the copies of Salmonella. Analyses of variance (ANOVA) were performed to test the significance of differences of the copies of Salmonella according to different factors (line, sex, experiment or isolator) using the anova function from base R (Type I sum of squares). We selected 48 chicks, either low or high Salmonella carriers based on data from ddPCR, balancing the experiment, sex, isolator and line factors.

RNA extraction and sequencing
Caecal tonsils from the 48 low/high carrier chicks were first grinded using ULTRA-TURRAX T25 (IKA). RNAs were then extracted using the NucleoSpin RNA Kit (MACHEREY-NAGEL) according to the manufacturer's protocol. The quantity of RNA was measured using a Nanodrop spectrophotometer (Thermo Scientific) and its quality was assessed using a 2100 Bioanalyzer Expert system using a total RNA nano Kit (Agilent), all samples displayed an RNA integrity number (RIN) > 7. RNAs were sent to the genomic platform GeT-Plage (Toulouse, France) for the cDNA library preparation (TruSeq Stranded mRNA kit, Illumina) and the sequencing (NovaSeq 6000 Sequencing System, Illumina). The sequencing data analysed during the current study are available in the NCBI Sequence Read Archive (SRA) database under the Bioproject accession number PRJNA649900.
The metadata and table of gene counts have been included as Additional files 6 and 7, respectively. PCA and ANOVA tests on PCA components were performed with the stats package (3.6.1). DESeq2 (version 1.26.0) was used with R version 3.6.1 to perform differential gene expression analyses [75]. In our model, the sex and the isolator effects were fixed. Genes with an adjusted p-value < 0.05 were considered to be differentially expressed genes (DEGs). The ggplot2 packages (3.2.1) were used for visualisation. The functional enrichment analyses were performed with topGO package version 2.38.1 [76].