Microbiota and Transcriptomic Effects of an Essential Oil Blend and Its Delivery Route Compared to an Antibiotic Growth Promoter in Broiler Chickens

This study evaluated the effect of the delivery of a commercial essential oil blend containing the phytonutrients star anise, cinnamon, rosemary, and thyme oil (via different routes) on broiler chickens’ ileal and ceca microbiota and liver transcriptome compared to an antibiotic growth promoter. Eggs were incubated and allocated into three groups: non-injected, in ovo saline, and in ovo essential oil. On day 18 of incubation, 0.2 mL of essential oil in saline (dilution ratio of 2:1) or saline alone was injected into the amnion. At hatch, chicks were assigned to post-hatch treatment combinations: (A) a negative control (corn-wheat-soybean diet), (B) in-feed antibiotics, (C) in-water essential oil (250 mL/1000 L of drinking water), (D) in ovo saline, (E) in ovo essential oil, and (F) in ovo essential oil plus in-water essential oil in eight replicate cages (six birds/cage) and raised for 28 days. On days 21 and 28, one and two birds per cage were slaughtered, respectively, to collect gut content and liver tissues for further analysis. Alpha and beta diversity differed significantly between ileal and ceca samples but not between treatment groups. In-feed antibiotic treatment significantly increased the proportion of specific bacteria in the family Lachnospiraceae while reducing the proportion of bacteria in the genus Christensenellaceae in the ceca, compared to other treatments. Sex-controlled differential expression of genes related to cell signaling and tight junctions were recorded. This study provides data that could guide the use of these feed additives and a foundation for further research.


Introduction
Over the years, the sub-therapeutic supplementation of antibiotic growth promoters (AGPs) has been used to preserve gut health, intestinal microbiota balance, and growth performance in the poultry industry [1]. This trend has now triggered both consumer and public health concerns bordering on the emergence of antibiotic resistance and antibiotic residues in the food chain [2,3]. Accordingly, a few country-specific restrictions on the use of AGPs are already in place, including in the EU [4], the US [5], and Canada [6]. The potential elimination of AGPs could exacerbate the risks of intestinal dysbiosis and bacterial diseases in poultry [7]. In the post-AGP era, understanding the complex interplay between the host and its intestinal microbiome signifies a critical step to achieving optimum gut health and intestinal microbiota balance in poultry.
The gastrointestinal tract (GIT) of poultry, populated by microorganisms in constant interaction with the host and digesta, is known to play a critical role in the host's growth and health. It is now evident that the proliferation of a balanced and beneficial gut microbiota population is vital to ensuring host protection against pathogenic bacteria and enhancing gut integrity and immunity [8][9][10]. Several factors, including environmental stressors [11][12][13], bird age [14], and nutrition [15][16][17] can modify the gut microbiota profile. Of all these factors, nutrition (including the type of diet and time of feeding) has been regarded as the main factor influencing poultry gut microbiota dynamics [9]. Given this modulatory role, several feed

Birds, Housing, and Diets
Hatchlings were weighed and randomly assigned to 6 new treatment groups ( Figure 1). Chicks (straight run) from the initial non-injection group were randomly allocated into 3 new treatment groups consisting of (A) chicks fed a basal corn-soybean meal-wheat-based diet (negative control treatment, NC), (B) chicks fed NC + 0.05% bacitracin methylene disalicylate (infeed antibiotics), and (C) chicks supplied the same commercial blend of EOs as earlier described via the water route (in-water essential oil) at the recommended dosage of 250 mL/1000 L of drinking water. The initial in ovo saline and in ovo essential oil groups were placed on the control diet to form treatments (D) (in ovo saline treatment) and (E) (in ovo essential oil treatment), respectively. The last treatment group, (F), consisted of chicks from the in ovo essential oil treatment group also supplied EO via the water route (in ovo + in-water essential oil treatment). All treatment groups had 48 birds each. Birds were placed in battery cages (0.93 m × 2.14 m), there were 6 birds per cage, and 8 replicate cages per treatment. Birds were reared for 28 d under uniform controlled environmental conditions in line with Cobb Broiler Management Guide recommendations. The room temperature was set at 31 • C on day 0 and gradually reduced to 23 • C on day 28, and relative humidity ranged between 45 and 55%. The ingredient and nutritional compositions of the basal diet used in the study are available in Oladokun et al. [47] and Supplementary Table S1. Birds were provided with feed and water ad libitum and diets were fed as mash throughout the rearing period which included the starter (0-14 d) and grower (15-28 d) phases. Diets met or exceeded the NRC [48] nutritional requirements for broiler chickens.
compositions of the basal diet used in the study are available in Oladokun et al. [47] and Supplementary Table S1. Birds were provided with feed and water ad libitum and diets were fed as mash throughout the rearing period which included the starter (0-14 d) and grower (15-28 d) phases. Diets met or exceeded the NRC [48] nutritional requirements for broiler chickens.

Sample Collection
On day 21, 1 bird per cage (8 replicate birds per treatment group) was randomly selected, weighed, and euthanized by electrical stunning and exsanguination. After slaughter, the small intestinal segment-the ileum (1.5-cm length mid-way between Meckel's diverticulum and the ileocecal junction)-was longitudinally opened, and digesta content was collected into microcentrifuge tubes. Aside from being the most studied small intestinal segments, the ileum microbiota was evaluated because reported trends suggest increasing microbial density in the distal region of the small intestine compared to the proximal regions as a result of longer digesta transit times [17,37].
Similarly, on day 28, 2 birds per cage (16 replicate birds per treatment group) were randomly selected and euthanized by electrical stunning and exsanguination. After euthanasia of the birds, digesta content from the pair of ceca was mixed and divided into two subsamples. One part was stored in plastic RNase-and DNase-free tubes placed in liquid nitrogen to analyze gut microbiota. The other part was placed in bio-freeze kits (Alimetric Diagnostics, Espoo, Finland) for the determination of short-chain fatty acids following published protocols [46]. Liver tissues (50-100 mg) using 1 mL TRIzol™ (Qiagen, Hilden, Germany) from 8 replicate birds/treatment were also rapidly collected on day 28 and promptly frozen in liquid nitrogen. All samples were stored at −80 °C until further analysis.

DNA Extraction, Qualification, Library Preparation, and Sequencing
The Qiagen DNeasy ® PowerSoil Pro Kit (50) (Cat. No./ID: 47014) was used to extract DNA from both ileal and ceca digesta contents. Digesta contents were allowed to thaw briefly at room temperature before subsequent DNA extraction, following the manufacturer's protocol. Briefly, 250 mg of digesta content was added to PowerBead Pro Tubes and then subjected to cell lysing steps involving vortexing and centrifugation. The retrieved lysate was then captured onto an MB Spin Column, followed by a series of purification and centrifugation steps. The MB Spin Column was then carefully placed into the provided 1.5 mL elution tubes from which extracted DNA was recovered. The concentration and purity of extracted DNA were subsequently determined by spectrophotometry (Nanodrop ND1000, Thermo Scientific). Extracted DNA samples (volume-50 µL,

Sample Collection
On day 21, 1 bird per cage (8 replicate birds per treatment group) was randomly selected, weighed, and euthanized by electrical stunning and exsanguination. After slaughter, the small intestinal segment-the ileum (1.5-cm length mid-way between Meckel's diverticulum and the ileocecal junction)-was longitudinally opened, and digesta content was collected into microcentrifuge tubes. Aside from being the most studied small intestinal segments, the ileum microbiota was evaluated because reported trends suggest increasing microbial density in the distal region of the small intestine compared to the proximal regions as a result of longer digesta transit times [17,37].
Similarly, on day 28, 2 birds per cage (16 replicate birds per treatment group) were randomly selected and euthanized by electrical stunning and exsanguination. After euthanasia of the birds, digesta content from the pair of ceca was mixed and divided into two subsamples. One part was stored in plastic RNase-and DNase-free tubes placed in liquid nitrogen to analyze gut microbiota. The other part was placed in bio-freeze kits (Alimetric Diagnostics, Espoo, Finland) for the determination of short-chain fatty acids following published protocols [46]. Liver tissues (50-100 mg) using 1 mL TRIzol™ (Qiagen, Hilden, Germany) from 8 replicate birds/treatment were also rapidly collected on day 28 and promptly frozen in liquid nitrogen. All samples were stored at −80 • C until further analysis.

DNA Extraction, Qualification, Library Preparation, and Sequencing
The Qiagen DNeasy ® PowerSoil Pro Kit (50) (Cat. No./ID: 47014) was used to extract DNA from both ileal and ceca digesta contents. Digesta contents were allowed to thaw briefly at room temperature before subsequent DNA extraction, following the manufacturer's protocol. Briefly, 250 mg of digesta content was added to PowerBead Pro Tubes and then subjected to cell lysing steps involving vortexing and centrifugation. The retrieved lysate was then captured onto an MB Spin Column, followed by a series of purification and centrifugation steps. The MB Spin Column was then carefully placed into the provided 1.5 mL elution tubes from which extracted DNA was recovered. The concentration and purity of extracted DNA were subsequently determined by spectrophotometry (Nanodrop ND1000, Thermo Scientific, Waltham, MA, USA). Extracted DNA samples (volume-50 µL, concentration-10-200 ng/µL) were then sent to the Integrated Microbiome Resource (IMR), located at Dalhousie University in Halifax, Nova Scotia, for library preparation and sequencing. Libraries of the V4-V5 hypervariable region of the 16S rRNA gene were prepared using universal primers 515 F (Illumina adapters + 5 GTGYCAGCMGCCGCGGTAA3 ) and 926 R (Illumina adapters + 5 CCGYCAATTYMTTTRAGTTT3 ) following protocols described by Comeau et al. [49].
Each sample was amplified with a different combination of index barcodes to allow for sample identification after multiplex sequencing. Library preparation and sequencing for all samples were performed with the Illumina MiSeq at the Integrated Microbiome Resource (http://imr.bio/, accessed on 15 July 2021) of Dalhousie University.

Short-Chain Fatty Acid Concentration and Total Bacteria Density
Ceca samples were collected using BioFreeze™ sampling kits (Alimetrics Diagnostics Ltd., Espoo, Finland) following the manufacturer's protocol. Samples were then subsequently submitted to Alimetrics Diagnostics 20007-1 (Espoo, Finland) for both SCFA concentration and total bacterial density quantification. The SCFA profiles were analyzed by gas chromatography (Agilent Technologies, Santa Clara, CA, USA) using pivalic acid as an internal standard. The acids quantified included acetic, propionic, butyric, valeric, and lactic acids. To quantify the total bacteria density, submitted samples were initially washed to remove solid particles and complex polysaccharides that may disturb subsequent DNA purification processes and downstream qPCR applications. The liquid phase was subjected to differential centrifugation for collecting the bacterial cells. The cell walls of the microbial cells were disrupted, and the chromosomal DNA was quantitatively extracted and purified using optimized protocol (Alimetrics Diagnostics 20007-1, Espoo, Finland). All measurements were performed with 16 replicates per treatment group.

RNA Extraction, Qualification, Library Preparation, and Sequencing
Total RNA in liver tissues was extracted on a QIAcube Connect using RNeasy Plus Universal Mini Kit (Qiagen, Cat. No. ID: 73404) following the manufacturer's instructions after disruption and homogenization were performed with a TissueLyser system. RNA elution volume was 30 µL. Total RNA was quantified, and its integrity was assessed on a LabChip GXII (PerkinElmer). Libraries were generated from 250 ng of total RNA and mRNA enrichment was performed using the NEBNext Poly(A) Magnetic Isolation Module (New England BioLabs). cDNA synthesis was achieved with the NEBNext RNA First-Strand Synthesis and NEBNext Ultra Directional RNA Second Strand Synthesis Modules (New England BioLabs). The remaining library preparation steps were performed using the NEBNext Ultra II DNA Library Prep Kit for Illumina (New England BioLabs). Adapters and PCR primers were purchased from New England BioLabs. Libraries were quantified using the Kapa Illumina GA with Revised Primers-SYBR Fast Universal kit (Kapa Biosystems). Average fragment size was determined using a LabChip GXII (PerkinElmer) instrument. The libraries were normalized, pooled, and then denatured in 0.05 N NaOH and neutralized using HT1 buffer. The pool was loaded at 225 pM on an Illumina NovaSeq S4 lane using Xp protocol as per the manufacturer's recommendations. The run was performed for 2 × 100 cycles (paired-end mode). A phiX library was used as a control and mixed with libraries at 1% level. Base calling was performed with RTA v3.4.4. The program bcl2fastq2 v2.20 was then used to demultiplex samples and generate fastq reads.

Bioinformatics and Statistical Analysis
The analysis of microbiota data was carried out using the Microbiome Helper pipeline (https://github.com/LangilleLab/microbiome_helper/wiki, accessed on 29 July 2021), based on QIIME2. This uses amplicon sequence variants (ASVs) created with Deblur. Primer sequences were removed from sequencing reads using cutadapt (v 1.14) [50], and primer-trimmed files were imported into QIIME2 (v. 2019.4.0) [51]. Reads (forward and reverse paired ends) were joined using VSEARCH (v 2.9.0) [52] and inputted into Deblur [53] to correct reads and obtain amplicon sequence variants (ASVs). Taxonomic assignment was performed with the SILVA database (v.1.3.2) using a naive Bayes approach implemented in the scikit learn Python library [49]. Rarefaction curves were used to examine the individual alpha diversity for all samples (with the default observed OTUs as the metric). Alpha diversity comparisons for the treatments were explored using boxplots and the Kruskal-Wallis statistical test set at p < 0.05. Beta diversity was visualized using weighted UniFrac PCoA plots. The relative abundance at different taxonomic levels was visualized using stacked bar charts, while significant microbiota proportions were determined in the Statistical Analysis of Metagenomic Profiles (STAMP) software [54] with an ANOVA test using the Benjamin-Hochberg false discovery rate as multiple test correction and then sorting by Corrected p-value (p < 0.05). Data on SCFA concentrations and total bacteria density were subjected to one-way ANOVA analysis in the Minitab statistical package (v.18.1). Data were analyzed in a completely randomized design and the analyzed data are presented as means ± SEM and probability values. Values were considered statistically different at p ≤ 0.05.
For the RNA-Seq analysis, adaptor sequences and low-quality scores containing bases (Phred score < 30) were trimmed from reads using Trimmomatic [55]. The resulting reads were aligned to the GRCg6 genome using STAR [56]. Read counts were obtained using HTSeq [57]. The R package DESeq2 [58] was used to identify differentially expressed genes between the groups. Nominal p-values were corrected for multiple testing using the Benjamini-Hochberg method. Gene ontology (GO) enrichment analysis was performed using the R package GOSeq [59]. Kyoto Encyclopedia of Genes and Genome (KEGG) Pathway Enrichment analyses of differentially expressed genes were performed on the PANTHER platform (http://pantherdb.org, accessed on 26 September 2021) [60].

Results
The 16S rRNA V4-V5 sequencing resulted in 8,774,523 quality read counts at an average of 60,934 counts per sample after quality filtering and demultiplexing. A total of 554 operational taxonomic units (OTUs) at the 97% sequence similarity level were obtained from all samples.

Microbiota Diversity
Internal sample α-diversity was estimated using the number of observed features (richness) and Shannon's index (diversity). Rarefaction curves of observed features and Shannon's index values reached a plateau in all samples, demonstrating that sequencing depth was adequate to cover the bacterial diversity in both ceca and ileal samples (Supplementary Figures S1 and S2).
Alpha diversity inspection revealed significant (p < 0.001) diversity between the ileal and ceca samples but not between treatment groups (Figure 2a-c). Ceca samples recorded a higher Shannon diversity index compared to ileal samples. While Shannon's diversity index showed a similar profile between the treatment groups in both ceca and ileal tissues, the ovo EO treatment recorded numerically higher alpha diversity in the ileum, the same as the NC treatment in the ceca.
To determine beta diversity, a principal coordinate analysis (PCoA) based on unweighted UniFrac distances was conducted. The PCoA plot showed unique cluster separation between the ileal and ceca microbiota; contrastingly, no difference in microbial community structure between treatments in both the ileum and ceca was observed (Supplementary Figure S3).

Microbiota Composition
The relative abundance of the predominant bacteria phyla and genus in both the ileum and ceca are presented in Figures 3 and 4, respectively. At the phyla level, ileum microbiota was dominated by Firmicutes (range of 99.5-99.8%), Proteobacteria (range of 0.03-1.79%), and Actinobacteria (range of 0.03-0.12%) for all treatments. Ceca microbiota phyla taxa showed a similar trend as the ileal microbiota, as the relative abundance of Firmicutes (range of 98.3-99.6%) was found higher than Proteobacteria (range of 0.38-0.81%), which was also higher than Actinobacteria (range of 0.01-0.24%) across the treatment groups. At the genus taxa, the ileal microbiota was 96% dominated by Lactobacillus, Clostridium sensu_stricto_1, Enterococcus, Romboutsia, and Lachnospiraceae_unclassified species, with Lactobacillus species being the prevalent species (occurring > 64% in all treatments, except for the in-water EO treatment, which recorded a 46.2% Lactobacillus relative abundance). Faecalibacterium was the most abundant genus in the ceca, recording at least 40% abundance across treatment groups. Similar to the ileal microbiota, genus Lactobacillus and Romboutsia were also found in the ceca, although at lower relative abundance. Contrastingly, the genus Lachnospiraceae was higher in the ceca (22.39%) compared to the ileum (1.91%). Significant differences in the cumulative proportions of bacteria in the genera Christensenellaceae_R-7_group, Elsenbergiella, Lachnoclostridium, and Shuttleworthia were recorded between treatments in the ceca ( Figure 5). Compared to other treatments, the in-feed antibiotic treatment significantly (p < 0.05) increased the proportion of Eisenbergiella, Lachnoclostridium, and Shuttleworthia. Contrastingly, the proportion of bacteria Christensenellaceae_R-7_group (p < 0.01) in the ceca was reduced by the in-feed antibiotic treatment when compared to other treatments. No significant differences in the microbiota proportion between treatments were recorded in the ileum at both the phylum and genus levels. ; (E) in ovo essential oil treatment-eggs injec with 0.2 mL of a saline + essential oil blend mixture at a dilution ratio of 2:1; and (F) in ovo + water essential oil treatment-chicks offered the essential oil blend via the in ovo and in water rou successively. Boxes in the boxplots denote interquartile range, solid middle line in the boxes den the median, and ⊕ denote the means, Symbols ○ and * in (a,b) represent outliers.

Microbiota Composition
The relative abundance of the predominant bacteria phyla and genus in both the eum and ceca are presented in Figures 3 and 4, respectively. At the phyla level, ileu microbiota was dominated by Firmicutes (range of 99.5-99.8%), Proteobacteria (range 0.03-1.79%), and Actinobacteria (range of 0.03-0.12%) for all treatments. Ceca microbi phyla taxa showed a similar trend as the ileal microbiota, as the relative abundance Firmicutes (range of 98.3-99.6%) was found higher than Proteobacteria (range of 0.3 0.81%), which was also higher than Actinobacteria (range of 0.01-0.24%) across the tre ment groups. At the genus taxa, the ileal microbiota was 96% dominated by Lactobacill Clostridium sensu_stricto_1, Enterococcus, Romboutsia, and Lachnospiraceae_unclassified s cies, with Lactobacillus species being the prevalent species (occurring > 64% in all tre ments, except for the in-water EO treatment, which recorded a 46.2% Lactobacillus relat abundance). Faecalibacterium was the most abundant genus in the ceca, recording at le 40% abundance across treatment groups. Similar to the ileal microbiota, genus Lactoba lus and Romboutsia were also found in the ceca, although at lower relative abundan Contrastingly, the genus Lachnospiraceae was higher in the ceca (22.39%) compared to ileum (1.91%). Significant differences in the cumulative proportions of bacteria in the g era Christensenellaceae_R-7_group, Elsenbergiella, Lachnoclostridium, and Shuttleworthia w recorded between treatments in the ceca ( Figure 5). Compared to other treatments, the feed antibiotic treatment significantly (p < 0.05) increased the proportion of Eisenbergie Treatments include (A) negative control treatment-chicks fed a basal corn-soybean mealwheat-based diet; (B) in-feed antibiotics-chicks fed NC + 0.05% bacitracin methylene disalicylate and (C) in-water essential oil-chicks supplied the essential oil via the water route at the recommended dosage of 250 mL/1000 L of drinking water; (D) in ovo saline treatment-eggs injected with 0.2 mL of physiological saline (0.9% NaCl); (E) in ovo essential oil treatment-eggs injected with 0.2 mL of a saline + essential oil blend mixture at a dilution ratio of 2:1; and (F) in ovo + in-water essential oil treatment-chicks offered the essential oil blend via the in ovo and in water route, successively. Boxes in the boxplots denote interquartile range, solid middle line in the boxes denote the median, and ⊕ denote the means, Symbols and * in (a,b) represent outliers.
Lachnoclostridium, and Shuttleworthia. Contrastingly, the proportion of bacteria Chris-tensenellaceae_R-7_group (p < 0.01) in the ceca was reduced by the in-feed antibiotic treatment when compared to other treatments. No significant differences in the microbiota proportion between treatments were recorded in the ileum at both the phylum and genus levels.
(a) (b) Figure 3. Ileal microbiota bacteria composition at the (a) phylum and (b) genus levels of broiler chickens subjected to different treatments groups. Treatments include (A) negative control treatment-chicks fed a basal corn-soybean meal-wheat-based diet; (B) in-feed antibiotics-chicks fed NC + 0.05% bacitracin methylene disalicylate; (C) in-water essential oil-chicks supplied the essential oil via the water route at the recommended dosage of 250 mL/1000 L of drinking water; (D) in ovo saline treatment-eggs injected with 0.2 mL of physiological saline (0.9% NaCl); (E) in ovo essential oil treatment-eggs injected with 0.2 mL of a saline + essential oil blend mixture at a dilution

Ceca SCFA Concentration
The resulting concentrations of ceca SCFA are presented in Table 1. Only the conce tration of butyric acid recorded a statistical trend towards significance (p = 0.09) in the water EO treatment, compared to other treatments. All other acids that were quantifi recorded no statistical significance between treatment groups (p > 0.05). Nonetheless, t in-water essential oil treatment equally recorded numerically higher concentrations acetic, lactic, volatile, and total fatty acids. Total bacteria (copies/gram of sample) we also found to be higher (p > 0.05) in the in-water EO treatment when compared to oth treatments. Table 1. Effect of essential oil delivery route on ceca short-chain fatty acid concentration (SCFA) a total eubacteria (copies/gram of sample) in broiler chickens.

Ceca SCFA Concentration
The resulting concentrations of ceca SCFA are presented in Table 1. Only the concentration of butyric acid recorded a statistical trend towards significance (p = 0.09) in the in-water EO treatment, compared to other treatments. All other acids that were quantified recorded no statistical significance between treatment groups (p > 0.05). Nonetheless, the in-water essential oil treatment equally recorded numerically higher concentrations of acetic, lactic, volatile, and total fatty acids. Total bacteria (copies/gram of sample) were also found to be higher (p > 0.05) in the in-water EO treatment when compared to other treatments. Table 1. Effect of essential oil delivery route on ceca short-chain fatty acid concentration (SCFA) and total eubacteria (copies/gram of sample) in broiler chickens.

Transcriptome Analysis
In this study, to identify differentially expressed mRNAs in the liver of broiler chickens, a total of 6,360,427,350 raw reads were generated from 48 samples (Supplementary Table S2). After the trimming step, 6,357,176,660 clean reads were obtained and the clean reads were aligned to the whole genome of Gallus gallus domesticus (Fasta: Gallus_gallus. GRCg6a.fa, Annotation: Gallus_gallus.GRCg6a.Ensembl98.gtf, source: Ensembl98). An average of 95.56% of clean reads were mapped to the genome. To determine the statistical significance of differentially expressed genes (DEGs) between treatments, the genes with the parameter of false discovery rate (FDR < 0.05) were considered differentially expressed genes/transcripts. Only a limited number of DEGs were observed (6: 3 up-regulated, and 3 down-regulated) to be differentially influenced by the treatments ( Table 2). The in ovo + in-water EO recorded the highest number (2) of DEGs in this category, as it down-regulated the expression of both cubilin (CUBN) and aldehyde dehydrogenase 1 family member L2 (ALDH1L2) genes. A heatmap illustrating the top 100 most variable genes is presented in Supplementary Figure S4.

Treatments include (A) negative control treatment-chicks fed a basal corn-soybean meal-wheat-based diet;
(B) in-feed antibiotics-chicks fed NC + 0.05% bacitracin methylene disalicylate; (C) in-water essential oil-chicks supplied the essential oil via the water route at the recommended dosage of 250 mL/1000 L of drinking water; (D) in ovo saline treatment-eggs injected with 0.2 mL of physiological saline (0.9% NaCl); (E) in ovo essential oil treatment-eggs injected with 0.2 mL of a saline + essential oil blend mixture at a dilution ratio of 2:1; (F) in ovo + in-water essential oil treatment-chicks offered the essential oil blend via the in ovo and in water route, successively. Each comparison is specified in the format "B vs. A", where group B is compared to group A, with group A being the denominator for the comparison. Liver tissues (50-100 mg) were sampled from 8 replicate birds/treatment (independent of sex) using 1 mL TRIzol™ (Qiagen, Hilden, Germany).

Discussion
The supplementation of phytogenic feed additives, especially essential oil, is reported to promote lipid and cholesterol metabolism [39], as well as to enhance immunity [61], leading to improved poultry performance. These favorable effects are thought to be exerted through the modulation of gut microbiota and the expression of several unique genes [39][40][41]44,61]. However, in vivo results on the effect of EO on chicken microbiota are somewhat inconsistent. While EO blends have been reported to reduce the relative abundance of pathogenic bacteria like Escherichia coli [22,23], Salmonella [24], and Clostridium perfringens [25] in broiler chickens, a few studies have equally reported no effect of EO supplementation on gut commensal bacteria [24,26,27]. This study utilized 16S rRNA gene sequencing in combination with transcriptomic analysis to investigate the effect of essential oil and its delivery routes (in water, in ovo, and in ovo + in water) and 0.05% bacitracin on both the composition and diversity of ileal and ceca microbiota and liver transcriptomics. Additionally, ceca short-chain fatty acid concentration was also evaluated. Bacitracin, the positive control in this study, is an extensively used antibiotic growth promoter in the poultry industry [62]. There is no doubt that the detailed delineation of the effect of a classic AGP like bacitracin and an alternative to AGP and its delivery routes, as this study presents, is key to understanding the molecular mechanisms underlying growth promotion in poultry. Accordingly, this study provides insight into the microbiota-mediated mode of action of antibiotics growth promoters, as well as preliminary transcriptomic evidence suggesting sex-controlled hepatic differential gene expression in broiler chickens offered antibiotics and essential oil (via water, in ovo, and in ovo + in-water delivery routes).
The gut microbiota plays an important role in host health, immune modulation, nutrient absorption, and pathogen control [63,64]. Although no treatment effect was recorded, this study revealed higher alpha (Shannon index) diversity in broiler chicken ceca compared to the ileum. This agrees with other studies [65][66][67], which have also recorded higher microbial diversity in the ceca. Higher microbial richness and stability observed in the ceca compared to the ileum in broiler chicken have been correlated with the higher number of obligate anaerobic microbes present therein, compared to aerobes or facultative anaerobes [63]. Consistent with the results of this study, both Abdelli et al. [68] and Pham et al. [69] have equally reported no effect of EO on alpha diversity. Thibodeau et al. [70] have shown that only extreme events (dysbiosis and disease inclusive) which modify the number of ecological niches in different bacterial species can alter the alpha diversity. Furthermore, beta diversity analysis revealed no difference in microbial community structure between treatments at both the ileum and ceca; however, the bacteria communities clearly differed across both gut sections. Other studies involving AGP or EO supplementation in broiler chickens have also observed similar results [71][72][73]. Conversely, Pham et al. [69] have recently highlighted the potential of EO to modulate the gut bacterial community structure. Aside from differences in intestinal sections, time of sampling, and supplemented additives, other factors, including broiler chicken breed, age, environmental condition, and disease status, can potentially cause shifts in beta diversity [63,74]. Diseases accompanied by intestinal dysbiosis like necrotic enteritis and Eimeria infection are reported to cause a significant change in gut microbiota community structure [75,76], buttressing the healthy state of the flock in this study.
Furthermore, in-feed antibiotic treatment significantly (p < 0.05) increased the proportion of Eisenbergiella, Lachnoclostridium, and Shuttleworthia, while decreasing (p < 0.01) the proportion of Christensenellaceae_R-7_group in the ceca, as compared to other treatments in this study. Interestingly, all of the bacteria with an increased proportion belong to the family Lachnospiraceae. The abundance of bacteria in the family Lachnospiraceae has been associated with improved weight gain [77], feed conversion ratio [78], and butyrate production in broiler chickens [79,80]. The performance result from this study published in Oladokun et al. [47] shows the increased weight gain recorded by the in-feed antibiotic treatment in the early period (d [1][2][3][4][5][6][7][8][9][10][11][12][13][14], compared to the in ovo treatments group, supports this result. Consistent with our results, Zhong et al. [81] reported the increased abundance of bacteria in the genus Eisenbergiella in neonates offered probiotics and antibiotics concurrently. An increased abundance in bacteria of the genus Eisenbergiella has been associated with reduced incidence of gastrointestinal disorders linked to metabolic and microbiota changes (functional dyspepsia), resulting in improved nutrient metabolism in broiler chickens [82]. Nonetheless, a few studies have also associated the abundance of this genera with the incidences of subclinical enteritis and Eimeria infection in broiler chickens [83,84], emphasizing the cost-benefit effects of antibiotic use in poultry production and the need for more studies in this regard. Lachnoclostridium has been positively correlated with increased butyrate production with attendant gut health protection and pathogen control effects [85,86]. Probiotics [87], prebiotics (wheat bran) [88], and antibiotics, but not essential oil, have all been reported to enrich the abundance of Lachnoclostridium in broiler chicken ceca [89]. Similar to other enriched genera in the ceca in this study, the genus Shuttleworthia has also been associated with increased weight gain and growth performance resulting from a possible role in lipid and carbohydrate metabolic pathways [77]. Furthermore, disease conditions like avian leukosis virus [90], coccidia infection [91], and Salmonella infection have all been reported to decrease the abundance of bacteria in the genus Shuttleworthia in broiler chicken ceca [92]. Similar to the results presented here, Hung et al. [93] have also reported a reduced abundance of members of the genus Christensenellaceae in the feces of weaned piglets offered the antibiotic bacitracin. Although the functional role of bacteria in this genus in the chicken microbiota is not fully known, their abundance has been associated with the colonization of Campylobacter jejuni, a foodborne zoonotic pathogen [70]. While the results presented here suggest that antibiotic growth promoters might give the birds a better growth advantage than EOs under this experimental condition, such growth advantage might come with some metabolic costs, deducible from the metabolic functions of bacteria species enhanced by this treatment. Hence, more research is needed on potent alternatives to antibiotic growth promoters with no reported adverse effects on the poultry industry. Nonetheless, the results on gut microbiota presented here provide a critical perspective on microbiota-mediated mode of action of antibiotics growth promoters in broiler chickens.
The fermentation of dietary fibers to yield SCFA constitutes an important function of the ceca commensal microbiota. No significant effect of evaluated treatments on ceca SCFA concentration was recorded in this study. Only the in-water EO treatment showed a statistical trend (p = 0.09) of enhancing ceca butyric acid concentration relative to other treatments. Butyric acid serves as an important energy substrate for the maintenance and proliferation of gut colonic cells and structures [71,94]. Essential oils (in feed or in water) have been reported to increase the concentrations of acetic, butyric, propionic, and lactic acids and total SCFA in quail breeders [95] and broiler chickens [96][97][98]. The positive effect of EO on ceca SCFA concentration could be related to the capacity of their phytogenic formulations to enhance bacteria proliferation in the lower gut. Several variable factors that potentially influence SCFA concentrations in broiler chickens, including the microbiota composition, bird age, and the amount and type of available fermentable substrates, explain the observed results in this study [23,99,100].
Transcriptomic analysis in this study suggests unique sex-controlled gene expression in broiler chicken livers. To evaluate the similarities and dissimilarities between samples in an unsupervised manner, principal component exploratory analysis (PCA) was carried out. The PCA showed that samples were not segregated by treatments (Supplementary Figure  S5) but instead showed modest segregation based on an unknown variable, probably sex. To confirm the hypothesis that treatments indeed clustered based on sex, the expression of five highly expressed genes on the W chromosome was examined. The results showed that all five genes showed much higher expression in group two than group one, indicating that group two samples were probably female (group two were samples with PC1 score > 0, group one were samples with PC1 scores < 0) (Supplementary Figure S6). Based on these PCA gene expression plots, samples were thus assigned as male or female according to their PC1 score. A total of 14 DEGs were found to be influenced by the treatment and sex (Supplementary Table S3). Of these DEGs, six genes were up-regulated (fold change range from 0.4 to 1.1) and eight genes were down-regulated (fold change range from −0.5 to −0.9). Sex-based analysis revealed that in male transcripts, antibiotic treatments recorded the highest number of DEGs (seven genes: four upregulated and three downregulated) compared to other treatments. Similarly, in male transcripts, the BVES (blood vessel epicardial substance) gene was significantly downregulated in both antibiotics, in-water essential oil, and in ovo essential oil treatments. In female birds, only four DEGS (two upregulated, two downregulated) were recorded amongst treatment groups. To understand the functional roles of identified DEGs, GO and KEGG pathway analysis was carried out. The GO analysis showed that a total of 33 significant GO categories were enriched (p < 0.05) compared to the negative control treatment (Supplementary Figure S7a,b). The GO terms included both biological process (BP), cellular component (CC), and molecular function (MF). However, a vast majority (75.8%) of the significant GO terms in the male transcript were observed in the antibiotic treatment, with GO terms in the BP category including "noncanonical Wnt signaling pathway" and "snRNA 3 -end processing" being the principal terms. In the female transcript, the vast majority (51.5%) of the significant GO terms were observed in the in-water EO treatment, with GO terms in the CC category including "signal transduction" and "plasma membrane" being the principal terms. The KEGG pathway analysis results, also shown in Table 2, provides predictions of differentially regulated pathways across treatments and sex. The results revealed that the main enriched pathways were cell signaling-(acting in sodium-glucose transporter, ion channels, exosome, and inositol phosphate metabolism) and tight-junction-related pathways. Other highlighted en-riched pathways include genetic information processing, transcription machinery, spindle formation proteins, phosphoric diester hydrolases, and purine metabolism.
Several factors, including the threshold of significance levels, but certainly not the number of animals or sample number utilized in this study, could have contributed to the low number of DEGs observed in this study (n = 14). Compared to this study, other broiler chicken RNA-Seq experiments [101][102][103] have utilized lower animal or sample numbers to detect a higher number of DEGs. The BVES gene with the cellular GO term category "regulation of microtubule cytoskeleton organization and tight junction related pathways" was found to be ubiquitously downregulated in male transcripts in this study irrespective of treatments, suggesting that this gene plays a vital metabolic function in the cell. Since first identified in 2001, the BVES gene has mostly been functionally correlated with the maintenance of epithelial integrity and tight junctions [104][105][106][107]; this has been validated by decreased trans-epithelial resistance values (TER, a measure of tight junction integrity) [106]. Nonetheless, the exact mechanisms of its role in tight junction maintenance are yet to be fully elucidated [108]. Osler et al. [106] have proposed that BVES's role in tight junction maintenance might indeed be a secondary effect, with more primary roles likely related to cell signaling, and structural support, among others. Besides the intestine, where tight junction proteins are noted to ensure gut barrier integrity, the BVES gene is also reported to be highly expressed in cardiac and skeletal muscle [109][110][111]. More recently, Gu et al. [112] reported the detection of BVES following whole-genome resequencing of the autochthonous Niya chicken breed and associated its function to the regulation of heart rate and heart development [113,114]. Considering that ischemic hepatic necrosis, which is linked to heart failure, could occur in broiler chickens [115] and the healthy state of flocks in this study, it is hypothesized that the downregulation of the BVES gene in broiler chicken liver observed in male transcripts in this study might be functionally related to the regulation of heart rate. Moreover, male embryos and adult chickens are reported to exhibit slower heart rates as compared to females [116,117]. More studies are thus needed to validate the relationship between BVES expression in the liver and heart rate regulation in broiler chickens.
Furthermore, antibiotic treatment upregulated the expression of INTS2 (integrator complex subunit 2), SLC5A10 (solute carrier family 5-member 10), CEP70 (centrosomal protein 70), and MED13 (mediator complex subunit 13) genes, while downregulating the expression of PDE11A (phosphodiesterase 11A) and CLCA1 (chloride channel accessory 1) genes in male transcripts in this study. Only in female transcripts did the antibiotic treatment upregulate PANX2 (pannexin 2) gene expression. INTS2 is a subunit of the integrator complex, which interacts with the C-terminal domain (CTD) of RNA polymerase II (RNAP II) large subunit and modulates 3-prime end processing of small nuclear RNAs (snRNAs) U1 and U2 [118]. The snRNAs are components of the spliceosome involved with the processing of pre-mRNA while also modulating the expression of other genes [119]. The modulation of snRNAs has also been reported to impact the innate immune system [120]. Slc5a10 (encoding SGLT5) is a mannose, fructose, and to a less degree, a glucose and galactose transporter [121,122]. Although glucose transporter 2 (GLUT2) is considered the main sugar transporter relevant to liver function [123], Fukuzawa et al. [124] have reported exacerbated hepatic steatosis induced by diminished sodium-dependent fructose uptake in SGLT5-deficient mice, suggesting the potential use of this gene as an indirect biomarker of liver health. Moreover, while the liver is the main site of ingested fructose metabolism, the occurrence of excess fructose beyond the liver's metabolic capacity triggers GLUT5 transporter upregulation to ensure fructose absorption into the epithelial cells [125]. Similar to the BVES gene with predicted cardioprotective effect, the upregulation of the MED13 gene by the antibiotic treatment is also thought to exert a cardioprotective effect in the birds. The MED13 gene is a component of the mediator complex, working in synchronization with RNA polymerase II to direct transcription [126]. Its mutation has been implicated in lethal cardiac defects [127][128][129]. Similarly, upregulated CEP70 expression, as induced by antibiotic treatment in this study, has been implicated in the pathophysiology of numerous cancers [130]. It is a centrosomal-associated protein that has been linked with the regulation of microtubule nucleation in animal cells [131][132][133]. Centrosome dysfunction has been linked to the incidences of liver diseases and other non-apparent cell cycle defects in humans [134]. The downregulated PDE11A is a dual-specificity phosphodiesterase that catalyzes the breakdown of the cyclic nucleotides cyclic adenosine monophosphate (cAMP) and cyclic guanosine monophosphate (cGMP) [135]. Although mainly expressed in the prostate, it also finds expression to a lower degree in the pituitary gland, heart, and liver [136]. Although CLCA1 is noted for its role in the activation of calcium-activated chloride channels, its downregulation is reported to enhance pro-inflammatory cytokine release in both mice mucus cells [137] and human small airway epithelial cells [138]. Increased innate immune responses are usually associated with increased energy demands; this suggests that antibiotic use might be an energy-intensive means of growth promotion. Similar to this study, Farmahin et al. [139] have reported the differential expression of PANX2 in the liver of female, but not male, Fischer rats. PANX2 is functionally known for its potential to create gap junctions that facilitate ion exchange between cells and their role as a potential tumor suppressor in the human brain, skin, and liver tissues [140][141][142][143].
In a like manner, in-water EO treatment equally downregulated the expression of the MC5R (melanocortin 5 receptor) gene in female transcripts in this study. MC5R encodes a protein receptor for melanocyte-stimulating hormone and adrenocorticotropic hormone. It has been functionally designated as a candidate gene for obesity and fatness in humans and domestic animals [144]. Consistent with the results presented here, Ren et al. [145] have previously reported that its expression in the chicken liver might be estrogen activated, further buttressing its differential expression in female transcripts in this study. Similarly, Blankenship et al. [146] have reported that the downregulation of MC5R was critical to achieving feed efficiency phenotype in first-generation female, but not male, quails in their study. This is likely achieved directly by fatty acid metabolism or indirectly by glucose homeostasis. This result is not unexpected, considering that the broiler chickens utilized in this study have been bred for high feed efficiency. Conversely, the in-water EO treatment up-regulated the expression of the GUCY2C (guanylate cyclase 2C) gene, which encodes guanylate cyclase belonging to the membrane guanylyl cyclase family [147]. Mice deficient in GUCY2C have been reported to have reduced inflammatory response due to reduced expression of pro-inflammatory molecules [148]. In contrast, higher expression of GUCY2C in the liver of milk-restricted lambs has been associated with increased pro-inflammatory response [149]. This is likely the molecular basis of the antibacterial properties of essential oils, especially as it relates to pro-inflammatory hepatic stimulus. This also has an energy trade-off, as more energy might be directed towards countering systemic inflammation and not growth. Higher expression of GUCY2C in human females as compared to males has also been reported [150]. Furthermore, while the in ovo EO treatment downregulated the expression of the cubilin (CUBN) gene, the in ovo + in-water EO treatment downregulated the expression of the MTMR6 (myotubularin related protein 6) gene in male and female transcripts, respectively. Although the functional relevance of CUBN in chickens is not fully understood, CUBN is generally noted to play a role in the uptake of vitamin, iron, and lipoprotein endocytosis [151,152]. Lee et al. [153] have reported the downregulation of the CUBN gene in chicken lines with high residual feed intake, suggesting a possible role in amino acid metabolism and molecular transport network. Sun et al. [154] have also alleged that the downregulation of this gene could be induced by stressors, particularly heat stress. More research is thus needed to fully elucidate the functionality of this gene in chickens. Overexpression of the downregulated MTMR6 by the in ovo + in-water EO treatment has been reported to inhibit the Ca 2+ -activated potassium channel [155,156]. Given the physiological function of the Ca 2+ -activated potassium channel, which is to regulate cellular membrane potential and calcium signaling, this is considered to be a beneficial effect of essential oil delivery via this route. Moreover, hypoxia has been reported to increase mitoBKCa channel activity (big conductance potassium channel) of rat liver [157]. Furthermore, in this study, both in-feed antibiotic and in ovo + in-water EO treatments downregulated the expression of ALDH1L2 in the liver of broiler chickens. Similar to this study, Li et al. [158] have previously reported the downregulation of ALDH1L2 in a nonalcoholic steatohepatitis (NASH) rat model by analyzing the liver proteome, suggesting that ALDH1L2 may be involved in NASH progression. Contrary to the result presented here, Bajagai et al. [101] have reported upregulation of the ALDH1L2 gene with continuous EO (2% oregano powder) supplementation in the liver of male broiler chickens. Differences in routes of EO supplementation and length of study may have possibly influenced reported results. Although little is known about the BTN1A1 gene in chickens, which was upregulated by the in-water EO treatment in this study, Huang et al. [159] reported that this gene might play a role in immune response via inhibition of T-cell activation using lipopolysaccharide-challenged broiler chickens [160,161]. Low hepatic expression of this gene has also been reported in water buffalo [162]. In addition, while the ST8SIA6 gene has been reported to be downregulated in the liver of apolipoprotein E-knockout (apo E-KO) mice offered phytosterol treatment for 14 weeks [163], an upregulation of this gene in birds of the in ovo EO treatment was observed in this study. Phytosterols are of plant origin and have reported cholesterol-lowering effects [164,165]. In humans, high expression of the ST8SIA6 gene has been attributed an oncogenic function, including tumor cell proliferation, invasion, and migration [166]. As little is known of this gene in chickens, an overt attribution of a high expression of this gene to an increased likelihood of hepatic steatosis or liver cancer might be farfetched. More studies are needed in this regard to enable a more definite prognosis. On the other hand, the CUBN gene is reported to play an important role in the metabolism and transport of the active form of vitamin D (1,25-dihydroxy vitamin D) in the liver. This has been confirmed in transcriptomics studies involving mice supplemented with cholecalciferol [167,168]. Although this gene is reported to be downregulated by the in ovo EO treatment in this study, Collision et al. [169] have previously reported its upregulation in mice liver under a trans-fatty acid (TFA)-induced non-alcoholic fatty liver disease challenge. Overall, the results presented here provide transcriptomic evidence on the possibility of "natural" phytobiotics (including essential oil) having side effects depending on the length of use, dosage, and administration routes, an important concept to be considered in the development of potent human and animal pharmacotherapeutic strategies.

Conclusions
Summarily, while treatments yielded no difference in alpha and beta bacteria diversity in this study, clear differences in ileal and ceca microbiota distribution and structure were recorded. In-feed antibiotic treatment is also reported to significantly increase the proportion of specific beneficial bacteria in the family Lachnospiraceae while reducing the proportion of bacteria in the genus Christensenellaceae, all in the ceca. No significant effect of the evaluated treatments on ceca SCFA concentration was recorded in this study. Only the concentration of butyric acid recorded a statistical trend towards significance in the in-water essential oil treatment when compared to other treatments. The study also suggests unique sex-controlled gene expression in broiler chicken liver. Compared to the negative control treatment, the differential expression of the INTS2, SLC5A10, MED13, CEP70, PDE11A, and CLCA1 genes functionally associated with genetic information processing, glucose transport, mediator complex, spindle formation proteins, phosphoric-diester hydrolases, and ion channel activity, respectively, were all regulated by the antibiotic treatment in male transcripts. Only the BVES and CUBN gene sets, functionally associated with tight junctions and cholesterol homeostasis, were regulated by the in-water and in ovo EO treatments in male transcripts, respectively, compared to the negative control treatment. Conversely, in female transcripts, while the antibiotic treatment regulated the expression of the PANX2 gene functionally associated with ion exchange, the in-water and in ovo + in-water treatments regulated the differential expression of GUCY2C, MC5R, and MTMR6 genes functionally associated with peptide hormone binding, melanocortin receptor activity, and peptidyl-tyrosine dephosphorylation, respectively, all compared to the negative control treatment. Taken together, the results presented here provide mechanistic insights on the microbiota-mediated mode of action of antibiotics growth promoters by modulating the abundance of specific bacteria communities, as well as preliminary transcriptomic evidence suggesting sex-controlled hepatic differential gene expression in broiler chickens offered antibiotics and essential oil (via water, in ovo, and in ovo + in-water delivery routes). To our knowledge, this is the first study to suggest such sex-controlled hepatic differential gene expression in broiler chickens offered these treatments. There is thus a need for well-designed in vivo studies that take sex into consideration in order to fully validate the results presented herein. Nonetheless, the data presented here not only provide guidance on antibiotics and essential oil application in the poultry industry; they also provide a solid framework for further research in the field.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/microorganisms10050861/s1, Figure S1: Rarefaction curves of observed features obtained from 16S rRNA gene V4-V5 sequences on the basis of (a) Microbiota source-cecum and ileum and (b) Treatment; Figure S2: Rarefaction curves of Shannon's index obtained from 16S rRNA gene V4-V5 sequences on the basis of (a) Microbiota source-cecum and ileum and (b) Treatments; Figure S3: PCoA plots based on unweighted UniFrac metric; Figure S4: Heatmaps of the top 100 most variable genes; Figure S5: Principal component analysis (PCA) plot; Figure S6: Gene expression of the chrW gene; Figure S7: Gene ontology (GO) classifications of differentially expressed genes (DEGs) between treatment groups in (a) male and (b) female transcripts; Table S1: Ingredients, calculated, and analyzed compositions of experimental diets; Table S2: Sequencing data quality control metrics; Table S3: Differentially expressed genes in the liver of broiler chickens as influenced by sex and treatment groups.  Institutional Review Board Statement: All experimental procedures were approved by the Animal Care and Use Committee of Dalhousie University (Protocol number: 2020-035), in accordance with the guidelines of the Canadian Council on Animal Care [45].

Informed Consent Statement: Not applicable.
Data Availability Statement: All data relevant to this manuscript are available upon request from the corresponding author.