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-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 61 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 61 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 61 (susceptible) in particular, several QTLs with low to moderate effects were identified [14-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 61 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 61, respectively resistant and susceptible to SE infection. The objectives of this study were to:

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;

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 61 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 Fig 1 according to the line, class and experiment factors. Means, standard deviations and p-values according to these groups are also given in Fig 1.

Fig 1

Summary of the results of comparisons between lines and between low and high carrier classes within lines according to the experiment.

Mean (red points) and standard deviation (SD; red bars) of the Salmonella Enteritidis abundance at 12 dpi in caecal contents (log10/g of caceal contents) of chicken 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 t-test 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].

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).

Fig 2

Principal component analysis of the gene expression in caecal tonsils in all samples.

The first principal component contains 29% of the variance and may be attributed to the line effect and the second principal component contains 25% of the variance and may be attributed to the sex effect.

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 61 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).

Fig 3

Principal component analysis of the gene expression in caecal tonsils from experiment 1 in each line.

(A) Principal component analysis (PCA) of the gene expression in caecal tonsils within line 61 in experiment 1. (B) PCA of the gene expression in caecal tonsils within line N in experiment 1.

A differential analysis with DESeq2 between low and high carriers within line 61 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.

Fig 4

Venn diagram.

Differentially expressed genes identified by comparing different groups of animals: (A) L6L1/L6H1 and LNL1/LNH1 (low vs high carriers of lines 61 and N, respectively, in experiment 1), (B) L6L1/L6H1 and LN/L6 (low vs high carriers of line 61 in experiment 1, and line N vs line 61), and (C) LNL1/LNH1 and LN/L6 (low vs high carriers of line N in experiment 1, and line N vs line 61).

DEGs common to intra-line and between-line analyses

When comparing DEGs identified between lines 61 and N, and within line 61 between low and high carriers, 38 genes were shared. Among these shared genes, 9 genes appeared to be regulated in the same direction: 4 DEGs (CEMIP, DMXL1, FUT2, NOS2) were upregulated in both the resistant line N (low carriage) and in low carriers in line 6, and 5 DEGs (PLEKHS1, RPS6KB2, CYP4B7, CYP2D6, CYP2AC1) were downregulated in both line N and in low carriers in line 61 (Fig 4 and Tables 1 and 2).

Table 1

Common differentially expressed genes downregulated between LN/L6 analysis and L6L1/L6H1 analysis (potentially associated with resistance to Salmonella).

			Between LN/L6		Between L6L1/L6H1
symbol	Description	id	LFC	padj	LFC	padj
FUT2	fucosyltransferase 2	ENSGALG00000001806	1.28	3.30e-03	1.11	0.05
DMXL1	Dmx like 1	ENSGALG00000002227	0.23	0.02	0.47	0.05
CEMIP	cell migration inducing hyaluronan binding protein	ENSGALG00000006413	0.52	6.03e-03	0.85	5.81e-04
NOS2	nitric oxide synthase 2	ENSGALG00000038096	0.46	7.17e-03	0.76	9.53e-03

Table 2

Common differentially expressed genes upregulated between LN/L6 analysis and L6L1/L6H1 analysis (potentially associated with susceptibility to Salmonella).

			Between LN/L6		Between L6L1/L6H1
symbol	Description	id	LFC	padj	LFC	Padj
CYP4B7	cytochrome P450, family 4, subfamily B, polypeptide 7	ENSGALG00000010469	-1.12	0.03	-2.51	5.73e-04
CYP2D6	cytochrome P450 family 2 subfamily D member 6	ENSGALG00000011894	-0.66	1.37e-05	-0.83	0.05
CYP2AC1	cytochrome P450, family 2, subfamily AC, polypeptide 1	ENSGALG00000016690	-1.49	0.01	-2.91	9.71e-04
PLEKHS1	pleckstrin homology domain containing S1	ENSGALG00000020679	-0.64	0.01	-0.91	0.02
RPS6KB2	ribosomal protein S6 kinase B2	ENSGALG00000031629	-0.17	0.04	-0.38	0.04

When comparing DEGs identified between lines 61 and N and within line N between low and high carriers, 560 genes were shared. Among these shared genes, 58 appeared to be regulated in the same direction: 43 DEGs were upregulated in both the resistant line N (low carriage) and in low carriers in line N, and 15 DEGs were downregulated in both line N and in low carriers in line N (Fig 4 and Tables 3 and 4).

Table 3

Common differentially expressed genes downregulated between LN/L6 analysis and LNL1/LNH1 analysis (potentially associated with resistance to Salmonella).

			Between LN/L6		Between L6L1/L6H1
symbol	Description	id	LFC	padj	LFC	padj
SPDEF	SAM pointed domain containing ETS transcription factor	ENSGALG00000002792	0.54	0.01	0.99	0.03
FNDC5	fibronectin type III domain containing 5	ENSGALG00000003567	0.48	0.01	0.71	0.04
CPD	carboxypeptidase D	ENSGALG00000004295	0.21	0.01	0.44	0.04
RBM41	RNA binding motif protein 41	ENSGALG00000004832	0.23	4.51e-03	0.31	0.04
S100B	S100 calcium binding protein B	ENSGALG00000006217	3.80	4.59e-65	1.01	8.66e-03
TRPM5	transient receptor potential cation channel, subfamily M, member 5	ENSGALG00000006521	0.75	1.72e-03	0.87	0.02
XKRX	XK related, X-linked	ENSGALG00000006637	0.41	1.87e-03	0.70	0.01
BCAS1	breast carcinoma amplified sequence 1	ENSGALG00000007796	0.43	3.63e-03	0.60	0.02
-	Uncharacterized	ENSGALG00000008309	1.16	5.10e-14	0.59	0.05
TCERG1L	transcription elongation regulator 1 like	ENSGALG00000010470	0.97	1.06e-14	0.41	0.04
TTC39A	tetratricopeptide repeat domain 39A	ENSGALG00000010540	0.41	0.02	0.67	0.02
RAB3B	RAB3B, member RAS oncogene family	ENSGALG00000010567	0.51	2.27e-03	0.80	0.04
SIRT5	sirtuin 5	ENSGALG00000012692	0.59	9.80e-20	0.41	0.02
SLC22A23	solute carrier family 22 member 23	ENSGALG00000012816	0.28	0.01	0.43	0.02
UROS	uroporphyrinogen III synthase	ENSGALG00000013688	0.14	8.61e-03	0.31	5.64e-03
ENDOUL	endonuclease, polyU-specific-like	ENSGALG00000014126	0.17	0.04	0.35	0.03
MOGAT2	monoacylglycerol O-acyltransferase 2	ENSGALG00000014170	0.48	1.42e-03	0.76	0.02
SIDT1	SID1 transmembrane family member 1	ENSGALG00000014812	0.41	0.02	0.75	0.02
GLDC	glycine decarboxylase	ENSGALG00000015053	1.44	1.e-12	0.66	0.01
SLC37A1	solute carrier family 37 member 1	ENSGALG00000016185	0.32	0.01	0.41	0.02
SH3YL1	SH3 and SYLF domain containing 1	ENSGALG00000016362	0.26	0.04	0.60	0.02
LOC427778	ectonucleoside triphosphate diphosphohydrolase 1	ENSGALG00000021274	0.34	0.02	0.62	0.04
MIR458	gga-mir-458a	ENSGALG00000021908	1.02	1.30e-10	0.61	0.03
HOXA4	homeobox A4	ENSGALG00000022622	0.26	0.01	0.48	0.03
ST3GAL4	chromosome 20 open reading frame 173	ENSGALG00000029747	4.86	0.01	0.66	1.68e-03
PAK6	p21 activated kinase 6	ENSGALG00000034605	0.34	0.04	0.78	0.02
AOC1	amine oxidase, copper containing 1	ENSGALG00000036190	0.73	0.03	1.07	0.05
CFAP43	cilia and flagella associated protein 43	ENSGALG00000036896	0.77	3.68e-08	0.46	0.04
ENPP3	ectonucleotide pyrophosphatase/phosphodiesterase 3	ENSGALG00000037587	1.51	3.79e-14	1.01	0.02
MYRF	chromosome 5 open reading frame, human C11orf9	ENSGALG00000041253	0.38	1.30e-03	0.56	0.04
BF2	Major histocompatibility complex class I antigen BF2	ENSGALG00000041380	0.89	2.16e-14	0.32	0.03
FAM8A1	family with sequence similarity 8 member A1	ENSGALG00000042858	0.13	0.04	0.30	0.02
GJC3	gap junction gamma-1 protein-like	ENSGALG00000042929	0.49	0.03	0.93	0.01
SYT7	synaptotagmin VII	ENSGALG00000043139	0.48	4.25e-04	0.63	0.04
-	Uncharacterized	ENSGALG00000043598	0.47	3.49e-05	0.53	0.01
RAPGEF3	Rap guanine nucleotide exchange factor 3	ENSGALG00000043775	0.65	4.05e-04	0.70	0.04
-	Uncharacterized	ENSGALG00000046455	2.49	1.62e-44	0.74	0.04
KLF1	Kruppel-like factor 1 (erythroid)	ENSGALG00000047208	0.61	6.69e-04	0.85	0.04
MUC4	mucin 4, cell surface associated	ENSGALG00000048048	2.40	5.84e-14	1.45	1.52e-04
LRRC66	leucine rich repeat containing 66	ENSGALG00000048601	0.38	5.60e-03	0.67	0.02
CLDN7	claudin 7	ENSGALG00000049326	1.27	4.86e-08	1.02	0.02
-	Uncharacterized	ENSGALG00000053586	2.20	3.38e-27	0.92	0.02
GZF1	GDNF inducible zinc finger protein 1	ENSGALG00000053659	0.47	8.90e-07	0.44	0.05
-	Uncharacterized	ENSGALG00000054057	0.17	0.03	0.36	0.03

Table 4

Common differentially expressed genes upregulated between LN/L6 analysis and LNL1/LNH1 analysis (potentially associated with susceptibility to Salmonella).

			Between LN/L6		Between L6L1/L6H1
symbol	Description	id	LFC	padj	LFC	padj
TGM3	transglutaminase 3	ENSGALG00000004804	-1.43	8.55e-04	-0.94	0.04
LECT2	leukocyte cell derived chemotaxin 2	ENSGALG00000006323	-2.48	5.47e-08	-0.91	0.04
JAKMIP1	janus kinase and microtubule interacting protein 1	ENSGALG00000015528	-0.30	0.04	-0.36	0.03
DEFB4A	defensin beta 4A	ENSGALG00000016669	-2.72	2.31e-07	-1.29	0.02
MTMR2	myotubularin related protein 2	ENSGALG00000017200	-0.20	4.51e-06	-0.21	0.02
MMP27	matrix metallopeptidase 27	ENSGALG00000019060	-0.54	1.23e-03	-0.48	0.03
CHST9	carbohydrate sulfotransferase 9	ENSGALG00000019489	-1.09	1.49e-05	-1.16	0.04
SERTAD4	SERTA domain containing 4	ENSGALG00000026359	-0.27	0.05	-0.57	3.22e-03
-	Uncharacterized	ENSGALG00000030245	-0.38	0.04	-1.41	5.34e-03
-	T-cell-interacting, activating receptor on myeloid cells protein 1-like	ENSGALG00000033375	-0.91	7.37e-04	-1.33	0.04
FKBP5	FK506 binding protein 5	ENSGALG00000042148	-0.44	2.10e-04	-0.36	0.04
EPX	eosinophil peroxidase	ENSGALG00000043254	-0.47	0.02	-0.64	0.03
SHLD1	shieldin complex subunit 1	ENSGALG00000049137	-0.32	2.25e-03	-0.56	0.04
-	Uncharacterized	ENSGALG00000050793	-0.66	1.11e-03	-0.75	0.04
-	uncharacterized LOC107049467	ENSGALG00000052009	-1.61	1.67e-04	-2.99	0.04

Functional enrichment analysis

Enrichment analyses were performed with topGO on DEGs between all animals from lines N and 61 and between low and high carriers within each line in experiment 1. The comparison between lines N and 61 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 61 (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).

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.

Table 5

BP GO terms shared between LN/L6 and L6L1/L6H1 analyses.

GO.ID	Term
GO:0009607	response to biotic stimulus
GO:0043207	response to external biotic stimulus
GO:0051707	response to other organism

Table 6

BP GO term shared between LN/L6 and LNL1/LNH1 analyses.

GO.ID	Term
GO:0003413	chondrocyte differentiation involved in endochondral bone morphogenesis
GO:0003416	endochondral bone growth
GO:0003417	growth plate cartilage development
GO:0003418	growth plate cartilage chondrocyte differentiation
GO:0006644	phospholipid metabolic process
GO:0030835	negative regulation of actin filament depolymerization
GO:0030837	negative regulation of actin filament polymerization
GO:0032272	negative regulation of protein polymerization
GO:0043242	negative regulation of protein complex disassembly
GO:0043244	regulation of protein complex disassembly
GO:0051494	negative regulation of cytoskeleton organization
GO:0051693	actin filament capping
GO:0098868	bone growth
GO:1901880	negative regulation of protein depolymerization

Table 7

BP GO term shared between L6L1/L6H1 and LNL1/LNH1 analyses.

GO.ID	Term
GO:0003008	system process
GO:0007169	transmembrane receptor protein tyrosine kinase signaling pathway
GO:0042592	homeostatic process
GO:0065008	regulation of biological quality

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 61. 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 61, 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 TLR1,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 61 [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 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 61, we identified many more DEGs and BP GO terms in line N, compared to line 61 (1034 DEGs and 288 BP GO terms in line N vs. 97 DEGs and 23 BP GO terms in line 61). Previously, we showed that the microbiota composition was significantly different between low and high carriers in line 61 but not between groups in line N [17]. One tentative hypothesis could be that the largest variability in Salmonella carriage in line 61 could be explained by differences in gut microbiota composition, as demonstrated with the PA12 chicken line [26]. The resistance against Salmonella in line 61 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 61 [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 61, 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 61 and are therefore associated with Salmonella resistance. Three of these genes have functions that are related to Salmonella: NOS2, DMXL1 and FUT2. First, NOS2 is a gene coding for the inducible nitric oxide synthase (iNOS) protein involved in macrophage inflammatory response [27]. Its function in the innate immunity against bacteria, viruses, fungi and parasites is well established, especially against Salmonella Typhimurium in mice models [28-33] and in in vitro models [34, 35]. In chicken, a transcriptome analysis of caeca showed an increase of NOS2 expression in the SE infected group [6], while genetic studies showed an association between NOS2 gene alleles and spleen SE bacterial load [36, 37]. An in vitro test showed the implication of iNOS between pathogen and macrophage cells during SE infection [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-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-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 61 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 61 and between low and high carriers in line 61 (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 61 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 61, 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 61 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 61 compared to line N.

Materials and methods

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 61 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.104 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 61, 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.

Transcriptome analysis

The quality control of raw reads was performed using the FastQC program. Trimming was performed using Sickle (version 1.33). Reads were then mapped on the Gallus reference genome (Gallus_gallus.GRCg6a.98.gtf) with STAR (version 2.5.3a) [73]. Genes were counted with htseq-count (version 0.12.4) [74].

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].

ANOVA analysis results on PCA1 and PCA2 components.

(TXT)

Click here for additional data file.

Metadata associated with all samples.

(TXT)

Click here for additional data file.

Gene count table.

(TXT)

Click here for additional data file.

DESeq2 results for the comparison of the N and 61 lines.

Distribution of raw p-values (A), dispersion plot (B), and MA plot highlighting 1492 and 1285 up- and down-regulated genes with a p-adj < 0.05 (C).

(TIFF)

Click here for additional data file.

DESeq2 results for the comparison of the low and high carriers within line 61 in experiment 1.

Distribution of raw p-values (A), dispersion plot (B), and MA plot highlighting 42 and 55 up- and down-regulated genes with a p-adj < 0.05 (C).

(TIFF)

Click here for additional data file.

DESeq2 results for the comparison of the low and high carriers within line N.

Distribution of raw p-values (A), dispersion plot (B), and MA plot highlighting 794 and 240 up- and down-regulated genes with a p-adj < 0.05 (C).

(TIFF)

Click here for additional data file.

Differentially expressed genes between LN and L6 groups.

(XLSX)

Click here for additional data file.

Differentially expressed genes between L6H1 and L6L1 groups in experiment 1.

(XLSX)

Click here for additional data file.

Differentially expressed genes between LNL1 and LNH1 groups in experiment 1.

(XLSX)

Click here for additional data file.

Differentially expressed genes between L6H2 and L6L2 groups in experiment 2.

(TXT)

Click here for additional data file.

Differentially expressed genes between LNL2 and LNH2 groups in experiment 2.

(TXT)

Click here for additional data file.

Enrichment analyse results performed on DEGs between line N and 61.

(XLSX)

Click here for additional data file.

Enrichment analyse results performed on DEGs between low and high carriers in line 61—experiment 1.

(XLSX)

Click here for additional data file.

Enrichment analyse results performed on DEGs between low and high carriers in line N—experiment 1.

(XLSX)

Click here for additional data file.