Burn injury alters the intestinal microbiome’s taxonomic composition and functional gene expression

Burn patients have a high risk of sepsis-related mortality even after surviving the initial injury. Immunosuppression increases the risk of sepsis after burn injury, as does the disruption of the intestinal epithelial barrier, which allows the translocation of bacteria and bacterial products into the circulation. The integrity of the intestinal epithelial barrier is largely maintained by the intestinal microbiota. Burn injury has been reported to result in significant changes in the intestinal microbiome composition. In this mouse study, we confirm these taxonomic differences in a full-thickness scald injury model using CF-1 mice. For the first time, we also address alterations in functional gene expression of the intestinal microbiota after burn injury to assess the microbiome’s physiological capabilities for overgrowth and pathogenic invasion: 38 pathways were differentially abundant between the sham and burn injury mice, including bacterial invasion of epithelial cells and gap- and adherens junction pathways.


Introduction
Burn patients, who survive the initial trauma, develop a state of relative immune compromise, making them highly susceptible to infections [1,2]. Burn patients also suffer from a loss of endothelial and mucosal integrity [3]. In combination with the relative immune suppression the outcome is detrimental (massive tissue edema, bacterial translocation, multi-organ failure) and in some cases lethal [4,5].
In healthy individuals, the gut microbiome protects against pathogenic microbes and helps maintain the intestinal epithelial barrier [6]. This commensal system, consisting of over 100 trillion microbes, also provides a number of other benefits to the host, i.e. metabolism and de novo synthesis of nutrients and assistance in immune development and function [7]. The gut microbiome is composed of between 300 and 1,000 species of bacteria. However, 99% of the total mass consists of only 40 species, among which three highly abundant phyla contribute the most to pathogen control and gut function: Firmicute, Proteobacteria and Bacterioidete [8,9]. Changes to the microbiome result in a dysbiosis, that disrupts this regulatory environment, enabling low abundant and pathogenic populations (i.e. Clostridium difficile) to expand [10,11]. Dysbiosis also contributes to the disruption of the intestinal barrier, enabling the translocation of bacteria and their products from the gut lumen to the mesenteric lymph nodes and on to the circulation [12]. Thus, the gut becomes a potential source of bacterial infections and sepsis after burn injury.
Several factors have been reported to affect the composition and biodiversity of the microbiome, including diet, environment, medication, infection/inflammation and burn injury itself [13][14][15]. On day one after burn injury, the gut is overgrown with gram-negative aerobic bacteria, including opportunistic pathogens [13]. We confirm this finding in a mouse model and expand the existing knowledge by analyzing gene expression of the whole microbiome to analyze their physiological capabilities for overgrowth and pathogenic invasion.

Mice
Male CF-1 mice were obtained from Charles River Laboratories (Wilmington, MA) at five weeks of age and allowed to acclimate for one to two weeks prior to conducting experiments. All mice were housed in standard environmental conditions and had ad libitum access to pellet diet and water. All experiments are approved by the Institution Animal Care and Use Committee (IACUC number 08-09-19-01) of the University of Cincinnati.

Scald burn injury
We used a scald burn injury model as previously described [16] Briefly, male CF-1 mice 6-8 weeks of age were randomized into two groups: scald and sham. Scald mice were anesthetized with 4.5% inhaled isoflurane in oxygen. Mice were placed in a template exposing their previously shaven dorsal surface and immersed in in 90.0˚C water for 9 s. The scald procedure yields a 28% total body surface area burn (total surface area calculated based on the Meeh formula [17]. It is a full thickness, third degree, insensate lesion, as the destruction of the entire thickness of the dermis and its peripheral sensory endings was confirmed histologically [18]. The procedure typically has a mortality of less than 10%. Mice were subsequently resuscitated intraperitoneally with 1.5 mL sterile saline. Sham treated mice underwent the same procedure expect for the immersion in water. All murine experiments were performed between 8 AM and 1 PM. After injury, the mice are monitored to ensure that they wake up from the anesthesia and subsequently allowed to recover on a 42.0˚C heating pad for 3 h. Mice are then returned to their home cage and monitored for any complications twice a day until the experiment is completed. No animal showed overt signs of morbidity that required euthanasia prior to the experimental endpoint.

Sample collection and isolation
On post burn day 1 (PBD1), mice were euthanized with CO 2 at a 30% fill rate and cecal fecal samples (100-200mg) were collected. Nucleic acid isolation was performed with the MoBio PowerMag Microbiome kit (Carlsbad, CA) according to the manufacturer's guidelines and optimized for high-throughput processing. All samples were quantified via the Qubit Quant-iT dsDNA High Sensitivity Kit (Invitrogen, Life Technologies, Grand Island, NY) to ensure that they met minimum concentration and mass of DNA.

Library preparation
Samples were prepared for sequencing with the Illumina Nextera XT kit and quantified with Quant-iT dsDNA High Sensitivity assays. Libraries were pooled and run with 100 bp pairedend sequencing protocols on the Illumina NextSeq 500 platform.

Data analysis
Raw data processing. Host sequences were removed with Kraken [19], which uses exact alignments of raw shotgun sequences to k-mers derived from the resus macaque reference genome. Remaining reads were processed with Trimmomatic [20] to trim adapter sequences and low-quality ends (<Q20). Reads shorter than 35bp after trimming were discarded. rRNA sequences from all three domains of life were identified and removed with SortMeRNA 2.0 [21]. Contaminant sequences, like PhiX174 and sequencing primers, were removed with Bowtie2 [22].
Taxonomic profiling using MetaPhIAn. MetaPhIAn2 (Metagenomic Phylogenetic Analysis, version 2.0, [23]) was used for the taxonomic profiling of the metagenomic samples. Raw non-host reads were used directly, because low quality reads are ignored in addition to human, 16S rRNA and tRNA reads. MetaPhIAn2 includes and expanded set of~1 million markers (184 ± 45 for each bacterial species) from > 16,000 reference genomes and > 7,500 unique species. Marker genes were identified for Bacteria, Archaea, viruses and eukaryotic microbes (Fungi and Protozoa) that are all crucial components of microbial communities.
Functional analysis. Filtered DNA sequences were mapped against a reference database of all proteins within the KEGG databases (version 75.0). The database is composed of 8,973,418 protein sequences extracted from 2,854 genomes. The search of translated DNA sequences was executed using Diamond [24] and hits that spanned ! 20 amino acids with ! 80% similarity were collected. In cases where one read matched these criteria against multiple proteins, only the protein or proteins (in the event of a tie) with the maximum bit score were collected. For pathway analysis, all proportional counts of KOs belonging to the pathway were summed for each pathway. Since one KO can be a member of multiple pathways (one-to-many relationship), this step increases the total proportional count in the pathway abundance table.
Alpha-diversity (within sample diversity) metrics. Richness is the sum of unique species found in each sample. Shannon diversity utilizes the richness of a sample along with the relative abundance of the present OTUs to calculate a diversity index.
Beta-diversity (sample-to-sample dissimilarity) metrics. All profiles are inter-compared in pair-wise fashion to determine a dissimilarity score and store it in a distance dissimilarity matrix. Distance functions produce low dissimilarity scores when comparing similar samples. Abundance-weighted sample pair-wise differences were calculated using the Bray-Curtis dissimilarity. Bray-Curtis dissimilarity is calculated by the ratio of the summed absolute differences in counts to the sum of abundances in the two samples [25].
Ordination, clustering and classification methods. Two-dimensional ordinations and hierarchical clustering maps of the samples in the form of dendrograms were created to graphically summarize the inter-sample relationships. Principal Coordinate Analysis (PCoA) is a method of two-dimensional ordination plotting that is used to help visualize complex relationships between samples. PCoA uses the sample-to-sample dissimilarity values to position the points relative to each other by maximizing the linear correlation between the dissimilarity values and the plot distances.
Whole microbiome significance testing. Permutational Analysis of Variance (PERMA-NOVA) is utilized for finding significant differences among discrete categorical or continuous variables. In this randomization/Monte Carlo permutation test, the samples are randomly reassigned to the various sample categories, and the between-category differences are compared to the true between-category differences. PERMANOVA utilizes the sample-to-sample distance matrix directly, not a derived ordination or clustering outcome.
Taxon, gene and pathway significance testing. The Wilcoxon Rank Sum test was employed to identify differentially abundant taxa, genes or pathways. Where samples could be paired across categories, a paired Wilcoxon Signed Rank test was employed. P-values were adjusted by Benjamini-Hochberg procedure to control for false discovery rates from multiple testing.

Burn injury severely impacts overall health on post burn day 1
Mice were subjected to a 28% full-thickness scald injury. To assess their overall health status after burn injury, we first assessed changes in body weight. On post burn day 1 (PBD1) mice had lost 3.3 g body weight on average, corresponding to approx. 10% of their total body mass ( Fig 1A). We also analyzed spleen-to-body weight ratios, as the spleen is the biggest single secondary lymphoid organ in mice and of great immunological importance. Spleen-to-body weight ratios dropped significantly in burn mice compared to sham-treated controls (Fig 1B). To further assess the systemic immunological response, we quantified serum G-CSF and IL-6 levels. Both were significantly increased upon burn injury (Fig 1C and 1D).

Taxonomic composition of the microbiome is altered upon burn injury
In order to analyze potential changes to the microbiome upon burn injury, cecal feces samples were collected on PBD1. The taxonomic richness (alpha diversity) was statistically similar Body weight (A) of mice was determined before and one day after burn or sham injury. Data are presented as weight change in g. Spleen-to-body weight ratios (B) were also determined. Serum G-CSF (C) and serum IL-6 (D) levels of the same mice were quantified by cytometric bead assay. For all experiments, each replicate is plotted (n = 7 mice/group), as well as mean ± SD. Asterisks indicate significant differences as assessed by unpaired, two-tailed Student's t test: ÃÃ p < 0.01; ÃÃ p < 0.001. between both groups: Similar sums of unique species were found in each sample (Fig 2A) and the relative abundance was similar as well (Fig 2B), indicating that the microbial community structure was not affected by burn injury. Nevertheless, beta diversity (sample-to-sample dissimilarity) was significantly altered between the two groups: Samples clustered together by injury type in the Bray-Curtis comparison, although one sham and two burn injury mice also clustered together (Fig 2C; p = 0.001, PERMANOVA).
Regarding taxonomic composition, fecal microbial communities of control mice were dominated by Firmicutes (Fig 3A) on the phylum level. In burn mice Firmicutes and Actinobacteria were reduced, whereas Proteobacteria, Deferribacteres, Bacteroidetes, Viruses and Apicomplexa were enriched (Fig 3A). Out of the top eight phyla, only Verrucomicrobia was not significantly altered (Fig 3A). On the family level, Enterobacteriaceae, Deferribacteraceae and Porphyromonadaceae were enriched in the burn injury mice samples, while Lactobacilaceae, Coriobacteriaceae, Clostridiaceae, Desulfovibrionaceae were reduced (Fig 3B). Ruminococcaceae were not altered (Fig 3B). Regarding individual strains, Adlercreutzia equolifaciens and Enterorhabdus caecimuris were enriched in the sham injury mice samples, whereas Parabacteroides goldsteinii, Mucispiillum schaedleri, Escherichia coli, Enterococcus faecalis, Oscillibacter unclassified, Mouse mammary tumor virus (genus Betaretrovirus) and Spleen focus forming virus (genus Gammaretrovirus) were enriched in the burn injury mice samples (Fig 3C).

signaling pathways are differentially abundant upon burn injury
We determined pathway richness upon burn injury and observed a significant increase in alpha diversity metrics in burn mice (Fig 4A). We also determined the abundance of 331 signaling pathways: The abundance of 4 pathways was significantly reduced in burn injured mice, 34 were significantly enriched (Fig 4B). Altered pathways regulate cellular processes, environmental or genetic information processing and metabolism, or are associated with Each replicate is presented, as well as mean ± SD. Asterisks indicate significant differences as assessed by Kruskal-Wallis rank sum test. Differentially abundant pathways in the cecal microbiome (B) were considered significant if their FDR-corrected p-value was less than or equal to 0.05 and the absolute value of the log-2 fold change was greater than or equal to 1. Infinite log-2 fold changes indicate that the other groups mean abundance is 0. Out of the 331 pathways tested, only the 38 pathways, whose abundance was significantly altered, are reported. human diseases, infection or specific organ systems. Among the 38 altered pathways, 10 were not expressed in the microbiome of sham injured mice.

Discussion
Previous work indicates that burn injury significantly alters the taxonomic composition of the microbiome [13] and that in mice, these changes were most pronounced one day after injury. Three days after burn, the dysbiosis had already decreased or normalized. We thus focused our analysis on PBD1. While the previous analysis addressed both the microbiome of the small and of the large intestine [13], we did not see increased permeability of the small intestine following burn injury [26]. Given that the cecum is the main site of fermentation and microbialderived metabolite production in the mouse, changes to the cecal microbiome have the greatest potential to impact murine health. We thus focused our analyses on the cecal microbiome.
Similar to the previous report [13], we observed an overgrowth with normally low-abundant gram-negative bacteria of the Proteobacteria, Deferribacteres and Bacteroidetes phyla, compared to the normal predominance of gram-positive Firmicutes. The increase in Proteobacteria was the most pronounced change, particularly regarding the subfamily Enterobacteriaceae. This corresponds to the previously published work [13]. The Enterobacteriaeae family includes many opportunistic pathogens which are commonly found in septic patients, i.e. Escherichia, Klebsiella, Proteus and Citrobacter [27], but it not yet known which individual strains of these bacteria elicit systemic inflammation after burn injury. We observed an increase in the abundance of Enterococcus faecalis, Oscilibacter, Escherichia coli, Parabacteroides goldsteinii and Mucispirillum schaedleri upon burn injury in our study. Enterococci are the third leading cause of nosocomical bloodstream infections, among which 60% are caused by E. faecalis [28]. Next to Escherichia, which is commonly found in septic patients [27], this identifies E. faecalis as a strong candidate which may induce systemic inflammation after burn. Individual reports of bacteraemia also exist for Oscilibacter [29] and P. goldsteinii [30]. Further studies will be necessary to determine which strains are the culprits in septic burn patients.
We also noted decreases in potentially protective bacteria. The majority of butyrate producing bacteria belong to the Firmicutes phylum, which was significantly less abundant in the burn injury mice. Regarding individual strains, the butyrate producer Enterorhabdus caecimuris was less abundant in the burn injury mice as well. Butyrate is well known for regulating T-cell function and the systemic immune response [31]. We have previously shown that colonic butyrate levels decrease following burn injury [26] and that butyrate reduces T cell death in an acid sphingomyelinase-dependent manner [32]. T cell depletion is a key component in burn-induced immunosuppression and enhanced susceptibility to opportunistic infections. Thus, butyrate restoration is a potential new therapy to ameliorate immunosuppression in burn. Indeed, we have previously shown that fecal microbiota transplants that reconstituted butyrate producers to burn injured mice ameliorate burn-induced colon permeability [26]. Taken together, our data motivates further research into probiotic supplementation with butyrate producing bacteria as a novel treatment for burn patients to hopefully prevent sepsis.
To our knowledge, our study is the first to assess functional gene expression of the microbiome upon burn injury, to determine the microbiome's physiological capabilities for overgrowth and pathogenic invasion. Pathway richness was significantly increased upon burn injury and the abundance of 38 pathways were significantly different between scald and sham mice. Among these, 9 metabolic pathways were significantly upregulated and 2 significantly suppressed.
Tetracycline biosynthesis was one of the suppressed pathways. Tetracyclines are antibiotics produced by Streptomyces (Actinobacteriae). A reduced abundance of this pathway corresponds to the reduced abundance of Actinobacteriae in burn injured mice and constitutes a loss of a self-regulation of the microbiome. The second downregulated pathway, N-glycan biosynthesis, can result in immune evasion, as non-glycosylated bacterial proteins are often no longer recognized by the host immune system [33].
Many of the metabolic pathways that were more abundant in the microbiome of burn injury mice have potential immunomodulatory consequences: Increased production of prostaglandin E2 -either through increased biosynthesis of its precursor alpha-linoleic acid, through increased arachidonic-acid biosynthesis or through increased ether lipid metabolism-suppresses the type 1 immune response [34].
LPS biosynthesis was also more abundant upon burn injury. This matches the expansion of Bacteriodetes found in burn injury microbiome samples. Bacteriodetes are major producers of LPS. Bacteriodetes-derived LPS has been reported to be immunoinhibitory, which may contribute to disease susceptibility upon burn injury [35]. Additionally, the metabolome of the gut microbiome has been previously shown to be predictive of host dysbiosis [36]. LPS biosynthesis was among the most predictive pathways, as was Chlorocyclohexane and chlorobenzene degradation, which was also significantly more abundant upon burn injury [36].
A large number of significantly affected pathways upon burn injury are linked to infection and pathogenic invasion. Particularly the pathways regulating the actin cytoskeleton, gap-and adherens junctions and bacterial invasion of epithelial cells would contribute to the disruption of the intestinal barrier and enable the translocation of bacteria and subsequent sepsis.
Our data confirms taxonomic differences to the intestinal microbiome upon burn injury, particularly the overgrowth of Enterobacteriaceae. We also identify functional gene expression changes upon burn injury, linking particular pathways to dysbiosis, lack of barrier integrity and pathogenic invasion for the first time. The disruption of the intestinal barrier following dysbiosis and the resulting translocation of bacteria and their products from the gut lumen to the circulation can lead to (Gram-negative) bacteremia, resulting in sepsis, multi-organ failure and even death of burn patients.
Supporting information S1 Table. Burn injury severely impacts overall health on post burn day 1-Supplement to  Table 1 includes individual data points for body weight before burn and on PBD1, spleen weight, spleen-to-body weight ratio, serum G-CSF and serum IL-6 levels for each sample. (XLS) S2  Table 4 includes individual data points on pathway alpha diversity and count tables for pathways and individual genes, as well as the mean log2 fold change values for each gene and pathway. (XLS)