An integrated transcriptomics and metabolomics study of the immune response of newly hatched chicks to the cytosine-phosphate-guanine oligonucleotide stimulation.

The immunological immaturity of the innate immune system during the first-week post-hatch enables pathogens to infect chickens, leading to the death of the animals. Current preventive solutions to improve the resistance of chicks to infections include vaccination, breeding, and sanitation. Other prophylactic solutions have been investigated, such as the stimulation of animal health with immunostimulants. Recent studies showed that administration of immune-modulators to one-day-old chicks, or in ovo, significantly reduces mortality in experimental bacterial or viral infection challenge models. Owing to a lack of molecular biomarkers required to evaluate chicken immune responses and assess the efficacy of vaccines or immune-modulators, challenge models are still used. One way to reduce challenge experiments is to define molecular signatures through omics approaches, resulting in new methodologies to rapidly screen candidate molecules or vaccines. This study aims at identifying a dual transcriptomics and metabolomics blood signature after administration of CpG-ODN (cytosine-phosphate-guanine oligodeoxynucleotides), a reference immune-stimulatory molecule. A clinical study was conducted with chicks and transcriptomics and metabolomics analyses were performed on whole-blood and plasma samples, respectively. Differentially expressed genes and metabolites with different abundance were identified in chicks treated with CpG-ODN. The results showed that CpG-ODN activated the innate immune system, within hours after administration, and its effect lasted over time, as metabolomics and transcriptomics profiles still varied 6 D after administration. In conclusion, through an integrated clinical omics approach, we deciphered in part the mode of action of CpG-ODN in post-hatch chicks.

Introduction

Early life is a period of elevated risk to infection as the immune system is still immature. Newly hatched chicks do not escape this observation. The chicken's innate immune system is already present at hatch, but not fully functional and displays lower cytokine production (Swaggerty et al., 2009; Lammers et al., 2010; Berghof et al., 2013). For example, heterophils become only fully functional after 2 to 3 wk of life (Swaggerty et al., 2009; Berghof et al., 2013) and chicks display a reduced humoral immune response until at least 2 wk of life (Lammers et al., 2010). Mounting a specific immune response after a primary infection or vaccination requires up to 2 wk, when broilers' life expectancy is about 6 wk. Thus, having an efficient innate immunity would help chickens to survive infections early in life. By contrast, cells of the innate immune system such as granulocytes, macrophages, or natural killer (NK) cells are fully functional post-hatch and capable to efficiently fight a broad range of pathogens.

Pathogens harbor specific pathogen-associated molecular patterns (PAMP) that are recognized by different classes of pattern recognition receptors (PRR), such as the toll-like receptors (TLR), that are expressed by immune cells, at their surface or intracellularly. In addition to providing early defense against infections, the innate immune response plays an essential role in triggering and driving the acquired immune system to respond effectively to infection (Medzhitov and Janeway, 2002; Akira, 2006). Several studies have shown that administration of immunomodulators, such as PAMPs, in one-day-old chicks or in ovo, significantly reduces mortality in experimental bacterial or viral infection challenge models. This suggests that they could be used in routine to boost the immune system and to prevent infections in chicken farms. Administration of synthetic oligodeoxynucleotides (ODN) containing CpG motifs (CpG-ODN), ligand of chTLR21, has been shown to efficiently protect different animal species, including poultry, against a variety of viral, bacterial, and protozoan pathogens (Klinman et al., 1999; He et al., 2007; Patel et al., 2008; Gunawardana et al., 2015). Studies in mammals and birds showed that administration of CpG-ODN alone increases antigen-presenting cells (APC) maturation, lymphocyte proliferation, heterophil degranulation, and production of antibodies, chemokines, Th1 cytokines, as well as inflammatory lipids and nitric oxide, resulting in decreased bacterial, viral, and intracellular parasitic infections (Zimmermann et al., 1998; Chen et al., 2001; Dalloul et al., 2004; Klinman, 2004; He et al., 2005, 2007; Linghua et al., 2007; Patel et al., 2008; Sanjaya et al., 2017). CpG-ODN has also direct stimulatory effects on chicken heterocytes, by enhancing production of both cytokines and chemokines (Kogut et al., 2006) and by inducing NF-κB-mediated leukotriene B4 and prostaglandin E2 production (Kogut et al., 2012). Treatment of NK cells with CpG-ODN increases their lytic activity and IFN-γ secretion (Ballas et al., 1996). Altogether, these findings suggest that treatment of chickens at day 1 post-hatch or in ovo with an immunostimulatory molecule such as CpG-ODN would improve their resistance to microbial infections.

A functional link between innate immunity and metabolism has already been established in animal models, as genes involved in immunity are also playing important biological functions in metabolism. For example, MyD88, involved in the signaling pathway of many TLRs, was shown to play a role in regulating metabolism, both in flies and mammals, as deletion of this gene in flies inhibits depletion of fat storage under starvation conditions (Ayyaz et al., 2013).

It has also been shown that inducing pro-inflammatory innate responses may alter metabolism and nutrient intake, putting chicken growth at risk (Matarese and La Cava, 2004; Afacan et al., 2012; Broom and Kogut, 2018). Thus, the complex interplay between metabolism and immunity emphasizes the need for studying both the metabolome and the transcriptome.

The goal of this work was to investigate the metabolomics and transcriptomics signatures associated with CpG-ODN immunostimulation in newly hatched chicks. The omics approaches appeared to be highly suitable to investigate deeply the global impact of immune-modulators on chicken health. To get a global view of the pathways targeted by CpG-ODN injection, NMR-MS metabolomics profiles and microarray transcriptomics signatures were investigated in blood from newly hatched chicks, 3 h, 3 D, and 6 D after injection, at day of hatch. The first objective of our study was to decipher pathways involved in CpG activity in vivo. The second objective aimed at developing an avian toolbox dedicated to the rapid screening of immune-stimulant candidates to be used as health boosters.

Material and methods

Reagents

Methanol, acetonitrile, LC grade water were obtained from Fisher. Phosphate buffer solution pH 7.4 of 0.1 mol was obtained from Sigma Aldrich. Deuterated methanol (MeOH-d4), deuterium oxide (D2O), and 4,4-dimethyl-4-silapentane-1-sulfonic acid (DSS) were obtained from Eurisotop (France). Protein precipitation 96-well plate Impact of 2 mL was obtained from Phenomenex (France). 10 kDa cutoff filter tube for ultrafiltration of 350 μL volume was obtained from VWR (France).

The Clinical Protocol

Animal handling was performed in strict accordance with good animal practice, as defined by the European relevant regulation, and was approved by the Ethical Committee of Merial. A total of 72 chicks from the JA957 conventional broiler chicken line were provided by a commercial supplier: Merial BI approved Duc Charolles Hatching Changy hatchery. They were received at the Boehringer Ingelheim animal facility at the day of hatch (D0 of hatch = D0 of the study). At the day of hatch, chicks were divided randomly into 2 groups. Birds were fed with Starvol Uni Initial (Dauphinoise de Nutrition Animale SA, France), with access to drinking water ad libitum. The birds were housed in a single box, battery cage, and at the same level of the battery. There was one cage per group, and groups were isolated from each other thanks to a Plexiglas wall. Chickens were then sampled randomly at the day of collection.

The first group was the control group. The second one received CpG ODN 2007 injected by subcutaneous route: 20 μg of CpG ODN 2007 (Sequence 5’-TCG TCG TTG TCG TTT TGT CGT T-3′), Avecia Biotechnology Inc. (USA). The dose of 20 μg/bird was selected based on literature and internal dose–response studies (data not shown). As described in Table 1, 3 h (D0+3h), 3 D (D3), or 6 D (D6) after injection, chicks (up to 8–12/group) were bled by intracardial route sampling and euthanized (by injection of Doléthal by intracardiac route [IC] in accordance with site Standard Operating Procedure). Blood was collected into Li-heparin tubes (BD Vacutainer Plus sterile N 368886). Chick sex was determined manually and confirmed by molecular biology.

Table 1

Final counts of chicks in each group and female (F)/male (M) ratios.

Group	D0+3h (F/M)	D3 (F/M)	D6 (F/M)
Control group	10 (4/6)	8 (5/3)	12 (6/6)
CpG ODN 2007 group	10 (3/7)	8 (5/3)	10 (3/7)

About 1 mL of blood was collected per chick; 50 to 150 μL was used for transcriptomics, about 800 μL for metabolomics, and 20 μL for chick-sex determination using PCR. For transcriptomics, blood samples were collected in Eppendorf containing 250 μL of Paxgene buffer (PreAnalytix). Tubes were turned around 10 times and then placed at room temperature for 2 h before being stored at −20°C. One day later, samples were stored at −80°C before processing. For metabolomics, after a rapid homogenization step, blood samples were stored in ice before centrifugation. The time between blood collection and centrifugation did not exceed 10 min. After a centrifugation step at 1,800 × g for 15 min at 4°C, plasma was collected, transferred to a cryotube, and stored at −80°C. Samples were collected and stored in MERIAL facilities and shipped to BIOASTER for analytical analysis. The final counts of chicks in each group and female (F)/male (M) ratios are given in Table 1. The design of the clinical protocol is given in Table 2.

Table 2

Study design.

Group	Chicks number	Treatment	Dose	Date
Control group	36	Control batch, no administration	-	-
CpG ODN 2007 group	36	Injection CPG subcutaneously	20 μg1	D0

Administered 20 μg/200 μL or 100 μg/mL per chicken.

RNA Extraction and Quantification

The extraction and stabilization steps of the RNAs were carried out using the PAXgene tube (Ozyme # BD762165) and the HT simply RNA Kit (Promega # AX2420), automated on the TECAN robotic platform. PAXgene was originally developed for humans, where a PAXgene reagent/blood ratio of 2.8 is recommended by QIAGEN (2.5 mL of blood for 7 mL of PAXgene solution). Not knowing whether these results could be transposed to chickens in particular with low volumes of blood, different amounts of blood (50–500 μL) and different PAXgene reagent/blood ratios (1; 2, 8; 5; 10; 20) were tested in a first step (Supplementary Data). To improve the quality of the RNAs, a DNase step was added after the extraction (RNA clean-up Nucleospin, Macherey-Nagel # 740948), to the reach the quality threshold of RIN >7. As aggregates were also observed for high values of PAXgene reagent/blood ratio, it was decided to use an intermediate ratio of 2.8. Finally, the appropriate protocol for the study was 50 μL of blood instead of 2.5 mL and a ratio equal to 2.8, with an additional DNase step (Supplementary Data). Using this protocol, all the 58 samples, except one (quantity = 4.2 ng/μL), passed the quality (7 ≥ RIN) and quantity (≥16.6 ng/μL RNA) criteria recommended in Affymetrix guidelines. The sample that did not reach Affymetrix requirements was nevertheless kept in the microarray experiment. Sample removal was based on Affymetrix QCs at the bioinformatics stage. The microarray selected in the project is the Affymetrix ChiGene-1_0-st that consists of 472,009 probes, 165,995 probesets (exon), and 20,804 transcripts, which corresponds to about 23 probes per transcript. Samples were then randomly distributed and analyzed in 11 batches.

Metabolomics Sample Preparation

The metabolites from plasma samples were recovered after protein precipitation with deuterated methanol in a 1:3 v/v ratio using 96-well filter plate Strata (Phenomenex) (Vuckovic, 2012). On contrary to protonated methanol, deuterated methanol allows preparing samples adapted for both NMR and LC-MS, as it does not affect metabolites identification by the LC-MS technique (Beltran et al., 2012). For protein precipitation, plasma to solvent ratios between 1:2 and 1:3 are commonly applied (Dunn et al., 2011; Sarafian et al., 2014). As there is no definitive consensus regarding plasma to solvent ratio in the literature, the 1:3 ratio was applied, in agreement with the recommendations of the precipitation plate supplier. The sample preparation was performed with an automatic liquid handler Freedom EVO150 (TECAN), equipped with 8 pipetting channels (4 with disposable tips and 4 steel needles), griper, T-vacs aspiration unit, and thermoshaker unit. A 96-well filter plate was filled with cold MeOH-d4 (900 μL) and then plasma (300 μL) was added. Next, the filtrate was divided into 2 aliquots, for NMR analysis (750 μL) and for MS analysis (200 μL). Filtrates were evaporated to dryness and the resulting residues were dissolved in appropriate solvents. The residues for MS analysis were dissolved in aqueous methanol solution (150 μL, 75%). The NMR residues were reconstituted in phosphate buffer solution (250 μL, 0.1 mol, pH 7.4) and further cleaned by ultrafiltration using a 10 kDa cutoff filter tube (VWR Centrifuge Unity). The filtrates (165 μL) were then mixed with phosphate buffer (55 μL, 0.1 mol, pH 7.4) in D2O containing 1.7 mmol/L DSS. The resulting mixture was transferred to 3 mm NMR SampleJet tubes (180 μL) for analysis.

LC-HRMS Bioanalyses

LC-HRMS analyses were performed on an Orbitrap Q-Exactive Mass Spectrometer (SN: 021192, Thermo Scientific) combined with a Dionex Ultimate 3000 RLS-HPLC system (SN: 8069641 & 8087153). Chromatographic separation was performed on an Xselect CSH C18 column (1 × 150 mm, 3.5 μm i.d) using water/acetonitrile (60/40, v/v) as solvent A, isopropanol/acetonitrile/water (88/10/2) as solvent B and water (100%) as solvent C. All 3 solvents contained 10 mmol/L ammonium acetate and 0.1% acetic acid as modifiers. The flow was set to 100 μL/min, the column temperature to 55°C, and the gradient was applied as follows: 50% A & 50% C for 1 min, 100% A at 8 min, 60% A and 40% B at 10 min, 30% A and 70% B at 25 min, 1% A and 99% B at 30 min, which was kept during 5 min and returned to initial conditions. This gradient was chosen to separate pro- and anti-inflammatory lipids as described by Jónasdóttir et al., 2015 on one hand and more non-polar lipids in a untargeted approach as described by Sarafian et al., 2014 on the other hand. Aliquots of 20 μL of each sample were mixed to prepare quality control (QC) samples that were injected regularly during the sequence to follow and correct intra- and inter-batch drifts. The injection volume was 5 μL and 2 separated injections were performed to acquire spectra in positive and negative ionization modes. In each mode, 2 batches were created to limit intra-batch drift during the acquisition. Parameters of the mass spectrometer were set as follows: Spray voltage: 3.0 kV in positive and −3.0 kV in negative mode, capillary temperature: 320°C, sheath gas: 45, aux gas: 5, aux gas heater temperature: 150°C. Data were acquired over an m/z from 100 to 1,125 at a resolution of 70,000. Furthermore, MS/MS spectra were acquired in a data-dependent acquisition mode, by fragmenting the 10 most intense signals at a resolution of 35,000.

NMR Bioanalyses

The NMR analyses were performed on 600 MHz Avance III HD spectrometer equipped with 5 mm QCI probehead and SampleJet autosampler from Bruker (Germany). To preserve the sample's metabolites integrity before analysis, the autosampler was thermostated at 5°C to 6°C. All samples were analyzed using noesygppr1d acquisition sequence for one dimensional proton spectra and sequentially were processed with the help of the Topspin 3.2 software from Bruker. For each analysis, 128 scans were cumulated within 14 ppm spectral width and 32 k time domain points. The mixing time and the relaxation delay were 80 μs and 4 s, respectively. The FID data were processed into 64 k final data points using 0.3 Hz exponential apodization and zero filling. Then, all NMR spectra were phase- and baseline-corrected before metabolites quantification.

Statistical Analysis

Transcriptomics Data Preprocessing

Four standard quality criteria were examined: (i) the amplification and hybridization Affymetrix controls, (ii) the median absolute deviation vs. the intensity median (MAD-Med) representation, (iii) the Normalized Unscaled Standard Error (NUSE), and (iv) the Relative Log Expression (RLE) (McCall et al., 2011). The 56 arrays were then simultaneously normalized using RMA (Irizarry et al., 2003). Probesets in the background noise, defined as those with intensities under 25 in at least one array, were excluded from the analysis of differentially expressed loci. This threshold corresponds to the intensity for which CVs across biological replicates were under 20%. Technical effects are commonly corrected using COMBAT (Johnson et al., 2007), a method treating batch effects as the combination of additive and multiplicative effects. However, because variance was homogeneous across batches, a simple linear model was used. After filtration, 19,666 transcripts (94.5% of the transcripts) were kept for the final analysis. The distribution of raw and normalized data stratified by gender and batches is given in the Supplementary Text.

Metabolomics Data Preprocessing

For LC-HRMS, the extraction and annotation of MS features were performed using Biotracs-OpenMS, an open-source bioinformatics tool based on OpenMS software (Sturm et al., 2008; Röst et al., 2016) that we recently published on GitHub (https://github.com/bioaster/biotracs-m-openms). Metabolite extraction was based on peak picking algorithms to annotate the most significant peaks using their mass-to-charge ratio (m/z) and retention time (RT). Peaks below a predetermined background-intensity threshold of 1e4 were removed from data. This threshold was chosen in agreement with the spectrometer specifications. The resulting list of features was then used for QC analysis to correct other analytical factors, in particular intra-batch and inter-batch drifts in BioTracs-Mimosa, a Matlab application designed for the processing of metabolomics data (https://github.com/bioaster/biotracs-m-mimosa). For intra-batch drift correction, the analytical drift is modeled using QC samples (Wehrens et al., 2016). Linear correction models were built without the first QC samples that are injected to stabilize the chromatographic separation (Gika et al., 2007). Inter-batch drift correction was performed by centering and scaling the 4 batches (Johnson et al., 2007) to remove additive and multiplicative inter-batch effects. Next, the coefficient of variation (CV = standard deviation/mean) of each feature was computed to assess their analytical stability in QC samples. Features having CVs >0.2 were considered as unstable in our analytical conditions and were removed from the study resulting in 509 m/z-RT features. Drift corrections was implemented using in-house bioinformatics solutions developed with Matlab. We assumed that correlated co-eluting features (RT shift < 10 s, correlation > 0.80) across all samples potentially correspond to the same metabolite and were denoted “isofeatures” here. Subsequently, the average signal intensity of each group of isofeatures was calculated for further statistical treatment. This allowed reducing the initial feature set to 334 isofeature groups. Within each isofeature group, putative adduct forms were assigned according to exact masses. The identification of the lipids was performed using the LipidMatch software that we wrapped in the BioTracs-LipidMatch application. This application is available online on GitHub (https://github.com/bioaster/biotracs-m-lipidmatch). Among all the putative metabolites extracted, 116 lipids (i.e. 22.8%) were unambiguously identified, mainly in the positive mode.

For NMR, the spectra were exported into the Chenomx NMR Suite software (version 8.2) for metabolites quantification. Based on Chenomx library spectra, the absolute concentrations of polar metabolites were quantified according to DSS internal standard and exported for all samples into an Excel table. All quantified metabolites concentrations were normalized and corrected according to the dilution factors applied during the sample preparation. For the ambiguous metabolites confirmation, additional analyses were performed by spiking authentic standards.

Transcriptomics Data Analysis

After correction of confounding effects, DEGs were identified using LIMMA, which relies on a moderated t-statistic (Smyth, 2004) and is well suited for small studies. CpG-ODN samples were compared to control samples at each time point. In this study, because we are looking for omics variations in whole blood, reflecting subcutaneous CpG-ODN administration, effect sizes are not expected to be large. For this reason, no threshold was set on log fold-changes, whereas a usual 0.05 threshold was set on the q-values (Q-value ≤ 0.05). These criteria were chosen so that a large enough number (>100) of probesets are found differentially expressed for functional analysis.

Metabolomics Data Analysis

Non-polar and polar data were analyzed separately. Multivariate partial-least-square discriminant analyses (PLS-DA) (Wold et al., 2001; Barker and Rayens, 2003) and univariate t-test analyses were assessed. The first approach allows extracting metabolic features and co-regulated metabolic features that contribute in discriminating the control from the CpG group while the last one does not account for co-regulations. A P-value threshold P ≤ 0.05 was chosen for the univariate t-test. The multivariate PLS-DA method and the univariate differential analysis (t-test) gave similar results (data not shown). The univariate analysis was therefore selected as it provides a more direct criterion, the P-value criterion, which can be fixed to select features of interest, and a t-score that represents the magnitude of the regulation.

Integrated Multi-Omics Analysis

The integrative methods are commonly classified into 3 categories: matrix factorization, Bayesian modeling, and network-based methods, as previously reviewed (Wei, 2015; Huang et al., 2017). In this study, we chose a matrix factorization approach, DIABLO (Singh et al., 2016), because of its ability to (i) perform variable selection, (ii) group correlated variables into a small number of latent variables and (iii) build these latent variables such that they maximize covariance between response and explanatory variables. DIABLO aims at building linear combinations of variables from each omics that are maximally correlated as specified in a design matrix. Because the present study is explanatory-oriented (the purpose being to select highly correlated variables regardless of their discriminatory power), we selected a full design matrix where all explanatory and responses data are connected. The level of sparsity (lasso penalty), allowing the selection of the most important features, was determined by leave-one-out cross-validation.

Pathway Enrichment and Functional Analysis

Owing to the lack of custom chicken pathway database, the functional analysis was done by orthology from well-characterized species. In this way, pathway enrichment was performed using Ingenuity Pathway Analysis (IPA, QIAGEN) on differentially expressed genes (DEG) while differentially expressed metabolites (DEM) were mapped onto the KEGG metabolic pathways. The choice of IPA was motivated by its ability to directly carry out enrichment from Affymetrix probesets. In our study, 10,975 out of 20,804 transcripts (about the half) were mapped in IPA and subsequently used for enrichment. The pathway enrichment step was performed using the Fisher exact test with a standard P-value threshold P ≤ 0.05.

Results

Clinical Study and Sample Analysis

CpG-ODN 2007, a prototypal immunomodulator, described as one of the most efficient agonists of chTLR21, which is the counterpart of mammalian TLR9 (Keestra et al., 2010; Sanjaya et al., 2017), was injected at day of hatch and systemic changes were monitored in blood at different time points post-injection. The design of the clinical protocol is presented in Table 1 and Table 2. An early time point of 3 h (D0 +3h) after injection was chosen to secure early signature after CpG-ODN injection (St. Paul et al., 2011). As the immunomodulatory properties of CpG-ODN have been detected up to 4 D after injection in the spleen, 2 additional time points, at day 3 (D3) and day 6 (D6), were chosen to assess its functional effects over time (He et al., 2007; Patel et al., 2008). The impact of CpG-ODN administration was monitored by transcriptomics and metabolomics analyses. It was compared to a control group that did not receive the immunostimulant.

From D0+3h to D6, chick's health was monitored daily. Ten chicks were used per group and no chick mortality was observed at D0+3h. Similar mortality rates were observed at D3 in the CpG-ODN-treated and control groups (3 per groups) and no additional death was noticed up to D6 in any group, with a general good health status of the birds.

At each time point, transcriptomics analyses were performed on whole blood samples and metabolomics studies on plasma samples. As shown in Table 1, males and females were not equally represented in each group. To avoid any bias due to animal gender in this study (Garcia-Morales et al., 2015), its impact was corrected during the analysis, as described in the Material and Methods section. A detailed description of collected samples is given in Supplementary Table 1. The sample collection and preparation were optimized to both extract RNAs and metabolites for untargeted analyses of transcriptomics and metabolomics data. Because birds' red blood cells are nucleated, and low volumes of blood were available, customized sample preparation was designed for this study. Our protocol, following data preprocessing, allowed the extraction of 19,666 transcripts, 47 polar metabolites by NMR, and 509 metabolic features by LC-HRMS.

Assessment of Confounding Effects

We first evaluated the relative weight of 3 experimental parameters (time, gender, and injection) before determining the omics signatures. Principal component analyses (PCA) were performed to assess putative confounding effects. Time, gender, and injection were projected onto the score plot to analyze how these variables shape molecular data. These analyses showed that time dominated the first component for both transcriptomics and metabolomics (Figure 1), highlighting the importance of the dynamical regulation of the metabolome and the transcriptome during the first days after hatch. Moreover, gender had an important effect on transcriptomics data, namely on the second principal component, while this confounding effect was not observed on metabolomics data. Gender effects were therefore corrected on transcriptomics data before computing differentially expressed genes. The correction was done using COMBAT (Johnson et al., 2007), a method treating batch effects as the combination of additive and multiplicative effects.

Figure 1

PCA score plot with transcripts (top), polar metabolic data (middle) and lipids (bottom). Each point represents a sample data collected from one chick. From left to right: projection of the sampling day (D0+3h, D3, D6 in yellow, red and blue, respectively), chick sex (male or female in blue and red, respectively) and group (CpG-ODN or control, in blue and red, respectively). Abbreviation: PCA, principal component analyses.

Transcriptomics Signatures

We found that 315, 2,783, and 3,246 transcripts were differentially expressed at D0+3h, D3, and D6 respectively. The DEG at time D0+3h, D3, D6 were subsequently mapped onto IPA (Ingenuity Pathway Analysis), leading to the identification of 205, 1,496, and 2,101 transcripts, respectively, that were successfully mapped in the Ingenuity Knowledge Base, with approximately 600 transcripts shared between the 3 time points (Supplementary Text, Supplementary Figure 5). There were 13, 109, and 207 pathways significantly modulated (P < 0.05) at D0+3h, D3, and D6, respectively. A heat map of the comparison analysis for the top 40 canonical pathways is presented in Figure 2. Canonical pathways are ranked from the highest total score (the negative log of the P-value derived from the Fisher's exact test) across the set of observations. Most of the identified pathways were related to immunological functions with a growing enrichment score (−log10 (q-values)) over time, and preferentially associated with a predicted pathway activation. As expected, many signaling pathways identified were associated to innate immunity (DC, macrophages) and their functions (phagocytosis, pathogen-related signaling [TLR]). Pathways of adaptive immunity were also highlighted, such as cellular immune response (TCR and BCR signaling), cytokine signaling, humoral immune response, apoptosis, and intracellular and second messenger signaling. Except apoptosis, all the observed pathways were upregulated.

Figure 2

Comparison analysis heat map for the top 40 canonical pathways. Canonical pathways are ranked to the top with the highest total score (the negative log of the P-value derived from the Fisher's exact test) across the set of observations (the circle size is proportional to the negative log of the P-value). The Z-scores (predicted pathway activation or inhibition scores) are indicated by the shade of color within the circles (activation: red; inhibition: blue). The criterion P ≤ 0.05 corresponds to −log10(P) ≥ 1.3.

Metabolomics Signatures

Following univariate differential analysis, we identified metabolites significantly regulated when chicks are injected with CpG-ODN. Globally, the lipid signatures can be stratified into 2 different groups according to the dynamical regulation of the lipids. A first group is composed of Fatty Acid Hydroxy Fatty Acid (FAHFA) that is upregulated early at D0+3h, silent at D3, and strongly downregulated at D6 (red box in Figure 3). These lipids were composed of FAHFA (2:0/16:2); FAHFA (2:0/18:2); FAHFA (2:0/20:2); FAHFA (5:0/17:2). The second lipid group displayed a different dynamical shape. These lipids were globally downregulated or silent at D0+3h, upregulated at D3, and downregulated or silent at D6. These were composed of 2 subgroups of lipids. The first subgroup was composed of lysophosphatidylcholine (LPC), plasmanyl-PC, plasmanyl-LPC mainly (blue box in Figure 3). They were strongly upregulated at D3 compared to the other lipids in this group but globally silent ad D6. The second subgroup of lipids was composed in majority of phosphatidylcholine (PC), phosphatidylethanolamine (PE), phosphatidylinositol (PI), plasmanyl-PE. They were characterized by a strong downregulation at D6. The list of lipids with their P-values and Z-scores are given in the Supplementary Materials.

Figure 3

LC-MS “non-polar” features regulated in the CpG-ODN group at time D0+3h, D3, D6, compared to the control group at the corresponding time points. The threshold of significance is P ≤ 0.05 (t-test). Negative (in blue) and positive (in red) Z-scores represent downregulation and up-regulation of the metabolites in the CpG-ODN group compared to the control group. Abbreviations: FAHA, Fatty Acid Hydroxy Fatty Acid; LPC, lysophosphatidylcholine; PC, phosphatidylcholine; PE, phosphatidylethanolamine; OxLPC, oxidized lysophosphatidylcholine. Pls-PC refers to plasmanyl-PC and plasmenyl-PC families, Pls-LPC refers to plasmanyl-LPC and plasmenyl-LPC families. The criterion P ≤ 0.05 corresponds to −log10(P) ≥ 1.3.

Few polar metabolites were upregulated at D0+3h after administration of CpG-ODN: hypoxanthine, lactate, pyruvate, and threonine (Figure 4). At D3, after administration of CpG-ODN, 3 fatty acids (3-hydroxybutyrate, acetate, isobutyrate) and one metabolite from the Krebs tricarboxylic acid cycle (fumarate) were downregulated. The highest differences on polar metabolites were found at D6 post-treatment: several amino acids were observed to be downregulated, which was already shown in literature (Kidd, 2004). One polar metabolite (creatine) was found to be upregulated at D6.

Figure 4

NMR polar metabolites in the CpG-ODN group at time D0+3h, D3, and D6, compared to the control group at the corresponding time points. The threshold of significance is P ≤ 0.05 (t-test). Negative (in blue) and positive (in red) Z-scores represent downregulation and upregulation of the metabolites in the CpG-ODN group compared to the control group. The criterion P ≤ 0.05 corresponds to −log10(P) ≥ 1.3.

Discussion

The functional interpretation of the molecular signatures allowed to cover the expected pathways and functions involved in innate immunity after CpG-ODN administration. In contrast to targeted analysis approaches, untargeted analyses cover a broader range of molecular mechanisms that could be involved in the regulation of immunity and host-pathogens interaction upon stimulation of the immune system. However, such large-scale untargeted omics data require specific preprocessing and adapted methods to extract meaningful biological information. Our analysis pipeline was therefore designed to reduce omics signatures, to the most relevant ones, without losing meaningful information. This strategy allowed extracting and clarifying transcriptomics and metabolic pathways functions in newly hatched chick upon CpG-ODN injection. It was shown that chick gender substantially affects the molecular signatures, but despite these confounding effects, the developed analysis methods successfully extracted the most significant transcriptomics and metabolomics signatures that correlated with the effect of CpG-ODN on the immune system. These signatures were therefore used to infer putative candidate biomarkers of the innate immunity, adaptive immunity, host/virus interaction, or anti-viral response.

Isolating RNA from avian whole blood is challenging because of the presence of nucleated red blood cells (RBC) that increase DNA and protein contamination. This was indeed observed in other species having nucleated red blood cells (Chiari and Galtier, 2011). The protocol we set up in this study allowed us to successfully extract RNA from whole chick blood. A DNase treatment was required to extract RNA, but this is in agreement with previously published works (Chiari and Galtier, 2011). For metabolomics, precipitation with methanol allowed covering a broad range of polar (e.g., amino-acids, organic acids, sugars) and nonpolar metabolites (lipids), from a single sample. With this method, all lipid classes were detectable, but recovery of lipids with low polarities, such as di- and triacylglycerides and ceramides, were compromised (Sarafian et al., 2014). This monophasic extraction was characterized by a higher reproducibility compared to biphasic extractions, which are classically used to separate lipids and polar metabolites from a single sample (Sarafian et al., 2014). Furthermore, methanol precipitation was described as an optimal extraction method for pro-resolving lipid mediators (Jónasdóttir et al., 2015), highly relevant in the context of our study. Hence, the method chosen was the best compromise between reproducibility and coverage of the relevant metabolome. This approach succeeded in addressing part of the mechanisms involved in innate immunity upon stimulation with CpG-ODN.

Interestingly, there is a close evolution of the differentially expressed genes associated with immune response over time. We found a signature of innate immunity together with a signature correlating with the activation of adaptive immunity at both D3 and D6. Several transcription factors, such as IRFs, AP-1, and NF-kB families, were differentially expressed, with NF-kB overexpressed at only D3 and IkB only at D6, confirming the activation of these transcription factors (Oeckinghaus and Ghosh, 2009). In agreement with this transcriptional activation, several genes encoding proteins involved in cell signaling were also differentially expressed. This included early transduction adaptors and/or kinases associated with T/B cell activation (such as Fyn, Lyn, Vav, SLP76, ZAP70, TEC, SYK) or with PRRs activation (such as Myd88, TRAF2, TRAF3, TRAF6, FADD, TOLLIP, RIP2). Downstream signal transduction kinases expressed in both adaptive and innate immune cells are upregulated, resulting ultimately in NF-kB activation through IKK complex. Consequently, we thought to find NF-kB targeted genes upregulated such as cytokines and chemokines (Akira, 2006). However, we failed to detect pro-inflammatory cytokines as a result of this innate immune activation, except BAFF and IL-1 which are downregulated.

We also evidenced the upregulation of several PRRs such as TLRs and RLRs involved in RNA detection as well as PKR (protein kinase R). Interestingly, TLRs recognizing others PAMPs are also upregulated, such as TLR2, TLR3, TLR4, TLR5, and TLR6. Of note, TLR21, the chicken TLR9, was detected but surprisingly not differentially expressed (Chrząstek and Wieliczko, 2014).

This immune activation was also accompanied by the upregulation of genes involved in apoptosis and rearrangement of the actin cytoskeleton (Wickramarachchi et al., 2010), which may reflect the recruitment of immune cells to the injection site or from the injection site to the draining lymph node.

Regarding metabolomics, the main lipids often associated with immunity are metabolites derived from arachidonic acid, eicosapentaenoic acid, and docosahexaenoic acid (e.g., prostaglandins, leukotrienes, resolvins). We therefore looked specifically for their signature in our data, but they were not detected in the analyzed plasma samples. It is possible that the concentrations in plasma were too low to be detected with our untargeted method. Targeted approaches able to detect those metabolites could therefore be needed to study the pro/anti-inflammatory response and it could be interesting to further investigate it in the context of CpG-ODN injection. In the end, the lipids significantly regulated after CPG-ODN injection were fatty acids and phospholipids, mainly the choline and ethanolamine forms.

At D0+3h, modulations were low in intensity, but we noticed the upregulation of a very specific lipid family, the Fatty Acid Esters of Hydroxy Fatty Acid. These lipids, discovered in 2014 by Yore et al., have been shown to stimulate insulin production and reduce inflammation in mammals (Yore et al., 2014). Interestingly, they have been found in egg yolk (Ma et al., 2015), which is internalized into the abdominal cavity of chickens during the last days of incubation and is an early source of lipids in the first days of growth (Noy and Sklan, 2001). Thus, by inducing inflammation, CPG-ODN may have induced an increased use of FAHFA, improving the metabolic state of young chicks.

At D3, plasmanyl-PC, plasmanyl-LPC, and LPC, which are part of the phospholipids, were significantly increased. These lipids are the main source of choline, also referred as vitamin B4, which has recently been claimed as a performance promoter in poultry (Igwe et al., 2015). Plasmanyl phospholipids, or plasmalogens, are an essential component of cells membrane. In chickens, they are important for growth of muscles, the choline form being overrepresented in the breast muscle (Miake et al., 2014). Moreover, in mammals, they are known to play a role as precursors of inflammatory mediators and have antioxidative properties (Wallner and Schmitz, 2011). Lysophosphopatidylcholine (LPC) is also interesting in the context of our study, as it has been both associated to growth performance and immunity stimulation in broilers (Allahyari-Bake and Jahanian, 2017). In the context of immunity, C16:0-LPC and C18:0-LPC, which were increased in our study, have been demonstrated to trigger significant and specific humoral-mediated immunity and therefore exhibit adjuvant activity (Bach et al., 2010). In addition, LPC induces, via the autotaxin (ATX) enzyme, the production of extracellular lysophosphatidic acid (LPA), which is an important phospholipid mediator in inflammation and immunity. It has also been shown in vitro that ATX is induced by TLR agonists, such as the TLR9 ligand CpG (Song et al., 2015). One adjuvant mechanism of CPG-ODN could therefore be directly related to the metabolism of LPC.

The results obtained at D6 showed that differential lipids were decreased, indicating that the metabolic effect induced by CPG-ODN was limited in time in the first days after treatment. Overall, the results on lipids showed that CPG-ODN induces the overexpression of lipids strongly involved in growth and immunity which is relevant with its adjuvant function and impact on growth.

For polar metabolites, at D0+3h, only some amino acids were increased. Threonine supplementation has been clearly associated to increased productivity in broilers (Al-Hayani, 2017), when pyruvate, lactate, or hypoxanthine are direct markers of an increased of metabolism. Thus, CPG-ODN has some direct early effect on global metabolism of chickens.

At day 3, organic acids (isobutyrate, hydroxybutyrate, acetate) were negatively correlated with CpG-ODN. These acids are known to play important roles in immunity. Isobutyrate was described to be closely linked to butyrate, which was shown to moderate the immune response by controlling NO production and some cytokines, mainly in the presence of agents that stimulate the immune cells. It was also shown to reduce tissue damage and improve the growth performance in chickens under stress (Zhang et al., 2011). In addition, isobutyrate originates from food (or microbiota) and can only be converted in butyrate by microbiota. Moreover, at day 6, the signatures were mainly associated with key amino acids that were described to be related to broiler immunity (Kidd, 2004). These amino acids are methionine, lysine, threonine, branched-chain amino-acids (valine, isoleucine, and leucine). However, the molecular mechanisms involved in the regulation of these metabolites in the immunity of chickens are not characterized. Our observations on lysine are particularly of interest as, in chickens, it is known to be a modulator of immunity (Kidd, 2004), to decrease lipogenesis in the liver when added to diet (Grisoni et al., 1991) and also to be the main amino acid represented in muscle (Munks et al., 1945). Thus, in chickens, amino acids modulation will affect the immunity response, the lipids metabolism and the muscle protein synthesis, which is important for growth. Overall, we can hypothesize that the decrease at D3 and D6 of all these immunity-related metabolites is a consequence of overstimulation happening at early time points (D0+3h and D3).

Lipids are known to be direct products of lipogenesis. But in contrary to mammals, where fatty acids are produced mainly by adipocytes, in chickens more than 90% of lipogenesis occurs in the liver (Leveille et al., 1975). This is in agreement with other findings (Dupont et al., 1999) showing that liver is the main genetic and metabolic site of chicken lipogenesis and ultimately fattening. As such, amino acids are an important substrate or regulator for lipid synthesis, as it has been confirmed in vitro by Caffin, 1994. In this context, CPG-ODN, by modulating the amino acids metabolism, could have an impact on lipid metabolism and therefore explains the shift in the response between amino acids and lipids.

Finally, creatine was the only metabolite upregulated by CpG-ODN at D6. It is strongly related to energy metabolism and is a key non-essential amino acid stored in muscles that have been previously associated to growth performance in broilers (Carvalho et al., 2013). It seems to be a direct translation of the early increase of the metabolic activity.

This study proved that systems vaccinology approaches enable to efficiently decipher the molecular mechanisms involved in the regulation of chick's immunity by CpG-ODN in the early days after hatch. The use of untargeted omics approaches allowed covering a broad range of molecules and then to improve our understanding of both transcriptomics and metabolomics processes involved in innate and adaptive immunity upon stimulation by CpG-ODN. However, integration of both omics would require to take into account the dynamical interplay between the transcriptome and the metabolome, to better extract co-regulations between these omics (Supplementary Text). This is particularly important on clinical data, where omics signatures are generally collected over days, on contrary to in vitro experiments, which are generally performed over hours (Jha et al., 2015). As transcripts at a given time point can be relevantly correlated to metabolites at another time point, longitudinal data could also allow to address this difficulty, as each individual is monitored all along the study, allowing to model the dynamical co-regulation of the transcripts and metabolite events. However, this approach may be difficult to implement for chicks, because of the low blood volumes available.

To our knowledge, this study is unique as it shows for the first time that it is possible to infer key processes of immunity after immunostimulation, by applying untargeted transcriptomics and metabolomics analyses on clinical blood samples of newly hatched chicks. It paves the way to the design of toolboxes that will allow accelerating the screening of health boosters. This work also highlights the benefit of using CpG as a health booster. By increasing the expression of genes involved in innate immunity and pathogen sensing, or by inducing the early production of lipids stimulating immunity, it lowers the threshold of the immune system to efficiently detect and eliminate infecting microorganisms.