Research Note: Snapshot of the transcriptome via RNA sequencing in the ileum of broiler chickens fed subtherapeutic concentrations of avilamycin.

Antibiotics have played a critical role in sustaining and improving livestock production in the past decades, but the emergence of antimicrobial resistance has led several countries to ban or limit their use. Since then, in-feed alternatives have gained a lot of attention but the development of efficacious alternatives implies a better understanding of the mode of action of antibiotic growth promoters (AGP) when administered at subtherapeutic concentrations. In the present study, 120 broiler chickens per group (8 pens/group) were fed for 35 d with either basal feed (control group) or feed supplemented with avilamycin (AGP group; 10 g/1,000 kg of feed). At the end of the trial, the ileum from the small intestine of 5 birds per group was sampled, and RNA were isolated for profiling their transcriptome via RNA sequencing (RNA-Seq). As expected, the growth of chickens in the AGP group was significantly higher than in the control group. Overall, 66 differentially expressed genes (false discovery rate ≤ 0.05 and fold change ≥ 2 or ≤ -2) were found in the ileum of chickens fed avilamycin in comparison with the control group. The functional analysis showed reduced activity of genes related to signaling by interleukins, with IL-22, SOCS3, and certain antimicrobial peptides found multiple times in these pathways in the AGP group at day 35. In addition, higher activity was predicted in a module of genes related to lipid metabolism and transport in the avilamycin group. The use of RNA-Seq allowed a snapshot of the whole transcriptome at day 35 and aimed at delivering additional data on the host-centric hypothesis regarding the mode of action of AGP (i.e. immunomodulation, reduction of the immunological stress).

Introduction

Soon after their discovery in the 1920s, antibiotics began to play critical role in contributing to the economic effectiveness of animal production. In addition to their therapeutic use against infections, supplementation at subtherapeutic doses is used to improve production and feed conversion efficiency (Castanon, 2007). However, the concern over antibiotic resistance and general consumer pressure have led to a decline in antibiotic growth promoter (AGP) use in some parts of the livestock sector including bans in some countries (e.g. EU since 2006). Achieving the practice of “No Antibiotics Ever” while maintaining conventional animal performance is challenging. Management changes include extra attention to the barn conditions (stocking density, downtime between flocks), breeders (sanitation and cleaning of the eggs), improved vaccination, and the use of in-feed alternatives to antibiotics (e.g. organic acids, enzymes, probiotics, phytogenic compounds, antimicrobial peptides). The in-feed alternatives have gained a lot of attention in the past years, and the identification and application of such alternatives is becoming an urgent issue for the global livestock industry. However, a lack of knowledge regarding the mechanisms of AGP action has hampered the development of efficacious alternatives. As reviewed by Brown et al., in 2017, current evidence on AGP mode of action has diverged into 2 main hypotheses: 1) bacteria centric and 2) host centric hypotheses. Most research has focused on the bacteria-centric hypothesis, suggesting that the changes in the microbiome create a more efficient system such as reducing competition for nutrients, preventing pathogen colonization, and selecting for bacteria better able to extract energy from the diet.

Given the extreme complexity of the interactions between gut bacteria, diet, and the host, it is likely that AGP also function by directly or indirectly affecting the intestinal physiology of the host. Evidence for nonmicrobial mechanisms of AGP action was reviewed by Niewold in 2007. These include the following: 1) the concentration of AGP used for growth promotion is normally lower than the minimum inhibitory concentration for enteric pathogens; 2) AGP consistently promote growth despite having variable effects on the microbiota; and 3) alternatives to AGP, although having known effects on the microbiota, have inconsistent effects compared with AGP.

Avilamycin is an antibiotic of the orthosomycin family known to inhibit the growth of gram-positive bacteria and has been widely used as AGP in poultry because it results in the best production responses compared with other antibiotics in poultry (Laxminarayan et al., 2015). Therefore, understanding the mechanisms of the most efficient AGP is expected to drive the development of rationale-based alternatives. With the emergence in the past years of next-generation sequencing technologies, scientists have now better tools to study the whole biological response to AGP.

In the present study, RNA sequencing (RNA-Seq) was used for investigating the transcriptome of broiler chickens after 35 days of feeding with avilamycin and aimed at providing new data regarding the interaction of AGP with the host.

Materials and methods

Chicken Feeding Trial, Tissue Collection, and Digestibility Analysis

The animal experiment was approved by the office of the Lower Austrian Region Government, Group of Agriculture and Forestry, Department of Agricultural Law (approval code LF1-TVG-39/024-2015). The feeding trial was carried out at the Center for Applied Animal Nutrition (Waxenecker KEG, Austria), in compliance with the European Guidelines for the Care and Use of Animals for Research Purposes and the Austrian Law for Animal Experiments.

Day-old broilers (Ross 308) were obtained from a commercial hatchery (Schulz, Austria) and randomly assigned to 2 experimental treatment groups. Each treatment group comprised of 120 birds divided into 8 pens with 15 birds/pen. The AGP treatment group received avilamycin (Surmax 100, Elanco and the diagonal bar logo are trademarks of Eli Lilly and Company, Liverpool, UK) supplemented diet (10 g/1,000 kg of feed), whereas the control group was on basal diet. All birds were housed in pens with litter covering the floor and had ad libitum access to the experimental feed.

At the end of the trial at day 35, broilers were euthanized by CO2 asphyxiation, and the ileum was removed from 6 birds per pen. The ileum was then flushed with sterile saline water, and digesta samples were pooled by pen (8 pens per group). In addition, the middle part of the ileum tissue from 2 of those 6 birds per pen (16 birds per group) was placed on ice and processed. After cutting and rinsing the ileum with sterile saline water, small pieces were transferred into 1 mL of RNAlater (Ambion Inc., Austin, TX), stored overnight at 4°C, and transferred to −80°C until transcriptomics analysis. The same procedure was applied at day 14 on samples from the ileum of one bird per pen (8 pens per group). From day 25 to 35, the diets contained 0.5% titanium oxide to determine the apparent DM digestibility in ileal digesta samples, and titanium oxide was measured using a photometric assay following the modified protocol previously described in (Reisinger et al., 2011).

RNA Extraction and RNA-Seq at Day 35

Approximately 30 mg of ileal tissue was disrupted via a bead-beating step, and the total RNA was extracted and purified with the RNeasy Mini Kit (Qiagen, Hilden, Germany). The average RNA integrity number was >9.0 for all the RNA extracted. The RNA from 5 birds per group (each from a different pen at day 35) were selected, and 2 technical replicates per bird were used for the RNA-Seq analysis (n = 10 per group). Construction of RNA-Seq libraries was performed with the TruSeq Stranded Total RNA Illumina kit including Ribo-Zero Gold. The final multiplexed libraries were sequenced on an Illumina HiSeq 4000 sequencing instrument (150 bp paired-end).

Processing of RNA-Seq Reads, Differential Expression, and Functional Analysis

FASTQ files of demultiplexed raw data of the 20 RNA-Seq samples were checked using the FastQC tool. BBDuk, part of BBTools suite, v38.22, was used to remove the adapters from Illumina data, and the trimmed reads were aligned to the chicken genome GRCg6a with STAR aligner, v2.7.0e. Raw counts per transcripts were obtained with HTSeq-count, v0.11.2, and normalization was carried out by the relative log expression method implemented in the DESeq2, v1.18.1, package. DESeq2 was also used to filter out the low-expressed genes based on default settings. To determine differences in gene expression, the generalized linear model using a negative binomial distribution was used in DESeq2 between the 2 groups (n = 10 RNA-Seq samples per group, technical replicate pairs were combined, and comparison was performed 5 vs. 5). The P values obtained were adjusted for multiple testing using the Benjamini-Hochberg procedure to control the false discovery rate (FDR). Genes were considered differentially expressed (DEG) when the fold change was ≥2 or ≤ −2, and with FDR ≤ 0.05.

Owing to the limitations of functional analysis with a small list of DEG, the Gene Set Enrichment Analysis (GSEA) was applied for this data set with the Bioconductor R package CEMiTool. The GSEA eliminates the dependency on the gene selection criteria by taking all gene expressions into consideration. It assumes that phenotypic differences are manifested by small but consistent changes in a set of genes. CEMiTool allows the identification and analysis of coexpression modules and, with a sample annotation file, performs the GSEA to visualize which modules are induced or repressed in the different phenotypes or groups. In addition, CEMiTool performs an over-representation analysis (ORA) to determine the most significant module functions via the clusterProfiler R package. Canonical pathway gene sets derived from the Reactome pathway database were used as the gene pathway list (gmt file).

qPCR on Selected DEG

The quantitative reverse transcription polymerase chain reaction (RT-qPCR) analysis was conducted for 3 DEG (SOCS3, avian beta-defensin-1 [AvBD1], and IL22), and their expressions were assessed on RNA samples collected from the ileum of birds at day 14 (n = 8 birds per group), and at day 35 on the same RNA samples used for RNA-Seq as well as on an extended list of RNA samples (the ileum from other birds sampled at the same time but not used for RNA-Seq; n = 16 birds per group). The primers were designed with the software Primer3 (SOCS3, Fwd: GAGAGAAGATCCCGCTGGTG, Rev: GGTGCCCGTTGACAGTCTTA; AvBD1, Fwd: TGCCCTTCCCTCACTCTCAT, Rev: GCTTGGGATGTCTGGCTCTT; IL-22, Fwd: GACGTCCTTCTCACAGAGGC, Rev: TTGGGGCATCAGTACAAGCA; HPRT1, Fwd: GACTGTGGACTTCATTAGGTTGA, Rev: TGGCCACTTTCACCATCTTTG; RPL32, Fwd: AGAAGTTCATCCGCCACCAG, Rev: CCTTGAACCTCCTGCGAACT). After cDNA synthesis using the Maxima H Minus First Strand cDNA synthesis kit including dsDNase (Thermo Fisher Scientific), the qPCR reactions were conducted on the Mastercycler ep Realplex (Eppendorf, Hamburg, Germany) using SYBR green chemistry (Kapa SYBR Fast Universal, Sigma). Relative gene expression was calculated using the 2−ΔΔCt method, with the geometric mean of the cycle threshold from the housekeeping genes HPRT1 and RPL32 for normalization.

Statistics Other Than for RNA-Seq

Statistical analysis was performed using GraphPad Prism, version 8, for Windows (GraphPad Software, La Jolla, CA). The Shapiro-Wilk method was used to test the normality of the data distribution. The GraphPad “nested t test” (fitting a mixed model) was used for analyzing the data from BW gain for which each treatment was tested in experimental replicates and each experimental replicate was assessed a number of times (technical replicates). “Unpaired t test” was used for the FCR, apparent ileal digestibility, and qPCR analysis.

Results and discussion

As expected, the continuous feeding of the AGP avilamycin resulted in a significant increase in BW gain at day 35 (P = 0.0082; Table 1), together with a significant increase (P = 0.023) in the apparent ileal digestibility of the DM compared with the control birds. Although the growth-promoting effect of avilamycin started to be seen after 25 d of feeding (but not after 14 d, data not shown; in agreement with Crisol-Martínez et al., 2017), this effect was stronger in the final phase of production. Therefore, we used the approach of RNA-Seq to determine the gene expression profile of those birds at the end of the growing phase at day 35.

Table 1

Performance measurements and DM digestibility in broiler chickens.

Replicates (= pens)	BW gain (d1-35, in grams)		Feed conversion ratio (d1-35)		Apparent ileal digestibility (d35, in %)
	Control	Avilamycin	Control	Avilamycin	Control	Avilamycin
Rep. 1	1,560 ± 348	1,891 ± 352	1.98	1.72	78.8	82.8
Rep. 2	1,727 ± 476	1,920 ± 331	1.73	1.85	81.3	82.2
Rep. 3	1,877 ± 313	2,004 ± 253	1.80	1.63	81.2	83.6
Rep. 4	1,740 ± 338	1,910 ± 343	1.64	1.53	84.4	80.8
Rep. 5	1,899 ± 299	1,887 ± 265	1.85	1.61	81.7	85.2
Rep. 6	1,809 ± 243	1,831 ± 269	1.87	2.06	81.8	84.8
Rep. 7	1,807 ± 211	1,932 ± 212	1.69	1.64	80.5	85.6
Rep. 8	1,826 ± 277	1,808 ± 495	1.57	1.86	83.2	84.7
Average1	1,779a ± 329	1,897b ± 320	1.77 ± 0.14	1.74 ± 0.17	81.6a ± 1.7	83.7b ± 1.7
AverageRNA-Seq	1,847 ± 176	1,951 ± 141	BWG average (±SD) of the 5 birds per treatment used for RNA-Seq analysis (each coming from separate pens)

Abbreviations: BWG, BW gain; FCR, feed conversion ratio; AID, apparent ileal digestibility.

Mean average of all the birds for BWG (average ± SD) and the mean average of the pens for the FCR and AID considering there was no individual bird measurement for feed intake and DM digestibility. Significant differences between treatments are indicated with dissimilar superscripts a,b (P < 0.05).

The libraries were successfully sequenced, generating 28.9 million to 52.6 million raw reads per sample, and with approximately 92% of the reads aligning to the chicken assembly GRCg6a. The snapshot of the whole transcriptome at day 35 revealed 66 DEG (with FDR ≤ 0.05 and fold change ≥ 2 or ≤ −2) in the ileum of birds fed avilamycin (41 DEG were downregulated, whereas 25 were upregulated). The full list of DEG from the GRCg6a analysis can be found at the indicated data repository on Gene Expression Omnibus. Because ORA takes a list of DEG as the input and that many tools using ORA recommend to import a minimum of 200 DEG, the GSEA, implemented in CEMiTool, was also performed. All the genes can be used in the GSEA, therefore, making it possible to detect situations in which all genes in a predefined set change in a small but coordinated way. CEMiTool generated eight modules containing 992 genes. Module 1 (M1) showed a significant and negative normalized enrichment score after the GSEA for the avilamycin group, suggesting reduced activity of the genes belonging to M1 (Figure 1B). The genes in M1 (254 genes) were associated with significant enrichment of IL 20 and IL-4 and IL-13 signaling (Figure 1C). Interestingly, this was in agreement with the ORA using only DEG showing similar signaling by IL (ORA performed in Reactome and Ingenuity Pathway Analysis; data not shown). The entities found multiple times in these pathways were the genes coding for IL-22 (9.46-fold decrease), suppressor of cytokine signaling 3 (; 4.67-fold decrease), matrix metallopeptidase 9 (MMP9; 2.45-fold decrease), IL-18 (2.10-fold decrease), and bactericidal/permeability-increasing protein (BPI; 6.25-fold decrease). This latter gene, bactericidal/permeability-increasing protein, was also found in the pathway associated with antimicrobial peptides, together with AvBD1 (6.14-fold decrease), defensin beta 4A (DEFB4A; 6.39-fold decrease), cathelicidin 3 (CATHL3; 7.18-fold decrease), and calgranulin C (S100A12; 3.44-fold decrease). Interestingly, the strongest effect in the RNA-Seq data set was seen on IL-22 expression. The IL-22 is a potent inducer and regulator of antimicrobial peptides from epithelial cells, including not only beta-defensins, lipocalin-2, S100, and regenerating protein family but also SOCS3 (reviewed in Ouyang and Valdez, 2008) and has been well characterized in avian species (Kim et al., 2012). This effect on IL-22 expression was very consistent when performing the RT-qPCR analysis on an expanded list of birds, namely all the sixteen birds initially sampled (2 birds/pen, 8 pens/group) and fed avilamycin (6.21-fold decrease, P < 0.001). The same effect in RT-qPCR was observed for 2 other DEG, SOCS3 and AvBD1 (3.24-fold decrease, P < 0.001, and 4.67-fold decrease, P = 0.005, respectively).

Figure 1

Modular coexpression analysis of the RNA-Seq data set with CEMiTool. (A) Profile plots of the modules M1 to M6. The black line represents the mean expression of all genes for that module across the samples. Samples are shown in the x-axis and colored by treatments (control and avilamycin). (B) Gene set enrichment analyses showing the module activity for each treatment. A positive normalized enrichment score (NES) indicates higher activity or expression of the genes from this module, and a negative NES indicates lower activity or expression of the genes from this module. (C) Over-representation analysis of modules M1 and M4. Bar graphs show the -log10 (adjusted P-value) of the enrichment between genes in modules and gene sets from Reactome pathway database. (D) Gene networks of modules M1 and M4. The most connected genes (hubs) are labeled and colored based on their “origin”: if originally present in the CEMiTool module, they are colored blue; if inserted from the interactions file, they are colored red. The size of the node is proportional to its degree.

Collectively, our findings showed that the activity of pathways related to the innate immune system (antimicrobial peptides, IL in host defense) was decreased by avilamycin at the end of the growing phase. This is in line with the hypothesis on host-centric mode of action and anti-inflammatory properties of AGP (Brown et al., 2017). By decreasing immunological stress in the intestinal mucosa, AGP would reduce the catabolic cost to the host, thereby increasing the energy available for muscle development and improving growth performance. However, as formulated for the first time by Niewold in 2007, this growth-promoting effect is reported for amount of AGP normally lower than the minimum inhibitory concentrations for enteric pathogens (therapeutic usage). The growth promoters are usually administered in relatively low concentrations (i.e. subtherapeutic), ranging from 2.5 to 125 mg/kg (10 mg/kg in the present study), depending on the compound and animal species (Laxminarayan et al., 2015). Subtherapeutic administration of AGP used in livestock species has been shown in mice to modulate the response to the bacterium Citrobacter rodentium (used as colitis model; Costa et al., 2011).

The functional enrichment analysis with CEMiTool has also revealed other small changes not seen with the typical ORA using the DEG. Indeed, the genes from module 4 (144 genes) showed higher activity in the avilamycin group and were associated with lipid metabolism and transport (fatty acid, peroxisomes, cholesterol, chylomicron; Figures 1B and 1C). The majority of lipid absorption occurs in the small intestine, and bile acids conjugated to fats aid in their transport. Increased ileal absorption of fatty acids has been already reported in broiler chickens using the same dosage of avilamycin (10 mg/kg of feed; Knarreborg et al., 2004). This was linked to reduced bile salt hydrolase (BSH) activity of Clostridium perfringens. Bile salt hydrolase is known to deconjugate bile salts and decrease lipid absorption, and Lactobacillus species are known to be major BSH-producers in the intestine. Bile salts are actively transported in the ileum, and therefore, the pathways from module 4 found enriched in the ileum of birds fed avilamyin suggest a connection to bile salts (they are synthesized from cholesterol, act as steroid hormones, and peroxisomes also play a role in the production of bile acids). Although we did not investigate the microbiota composition in the ileum, recent studies using 16S rRNA sequencing have reported contrasting results of subtherapeutic administration of avilamycin (15–25 mg/kg of feed) on both the ileum and ceca content (Crisol-Martínez et al., 2017; Choi et al., 2018). Interestingly, Choi et al. (2018) showed that unlike the ceca, the bacterial communities were influenced by avilamycin in the ileum. However, 2 known BSH-encoding Lactobacilli, Lactobacillus crispatus and Lactobacillus reuteri, showed opposite effects such as decreased and increased, respectively, in the avilamycin group. Lower abundance of the genus Lactobacillus was reported in the ceca of birds fed zinc bacitracin but not fed avilamycin (Crisol-Martínez et al., 2017). Discrepancies in the effects of several AGP on the microbiota in livestock species have been summarized by Brown et al. (2017), showing that the effects were not consistent.

The effects of AGP on host response have not been extensively studied, but it does not seem that microbiota changes alone (as initially pointed out) are fully explaining the growth-promoting effect of AGP. Investigations at earlier time points (RT-qPCR at day 14 in the present study did not show significant differences in the expression of SOCS3, IL-22, and AvBD1 [-1.57, -1.75, and +2.23, respectively]) would be needed to fully understand the effects of AGP over time.