Transcriptional Profiling of the Small Intestine and the Colon Reveals Modulation of Gut Infection with Citrobacter rodentium According to the Vitamin A Status

Vitamin A (VA) deficiency and diarrheal diseases are both serious public health issues worldwide. VA deficiency is associated with impaired intestinal barrier function and increased risk of mucosal infection-related mortality. The bioactive form of VA, retinoic acid, is a well-known regulator of mucosal integrity. Using Citrobacter rodentium-infected mice as a model for diarrheal diseases in humans, previous studies showed that VA-deficient (VAD) mice failed to clear C. rodentium as compared to their VA-sufficient (VAS) counterparts. However, the distinct intestinal gene responses that are dependent on the host’s VA status still need to be discovered. The mRNAs extracted from the small intestine (SI) and the colon were sequenced and analyzed on three levels: differential gene expression, enrichment, and co-expression. C. rodentium infection interacted differentially with VA status to alter colon gene expression. Novel functional categories downregulated by this pathogen were identified, highlighted by genes related to the metabolism of VA, vitamin D, and ion transport, including improper upregulation of Cl− secretion and disrupted HCO3− metabolism. Our results suggest that derangement of micronutrient metabolism and ion transport, together with the compromised immune responses in VAD hosts, may be responsible for the higher mortality to C. rodentium under conditions of inadequate VA.


Introduction
Diarrheal disease causes~0.8 million deaths per year in children under the age of five and ranks as the second leading cause of infection-related mortality in this demographic group [1,2]. Recurrent early childhood diarrhea, occurring during the period of rapid growth of the brain and other organ systems, can contribute to lasting impairment in fitness, growth, and cognition. According to a report of Petri et al. [3], repeated diarrhea in the first two years of life was associated with a loss of 10 IQ points and one year of incident with diarrhea prevalence. VA deficiency impairs the development, vision, immune functions, and is associated with the infectious disease severity [42]. It is estimated that VA deficiency impacts over 250 million preschool-aged children [45]. In the present study, we extended the current knowledge by seeking to identify the potential key mediators responsible for the resistance of VAS mice to the C. rodentium infection and the effects of interaction between the infection and the VA status on the global gene expression profiles in both the colon and the SI. Because the colon is a major site of the C. rodentium infection and the SI is a major site of VA metabolism, antigen presentation, and immune cell maturation, RNAseq was performed using both organs. We used the RNAseq analysis and a bioinformatics approach as a global and unbiased means to identify the genes and pathways affected by C. rodentium and as modified by the VA nutritional status. We hypothesized that the characterization of differential expression, co-expression networks, and enrichment analyses would identify the key genes and pathways involved in diarrheal disease that interact with the VA status.

Animals
C57BL/6 mice were bred and maintained according to the guideline of the Institutional Animal Care and Use Committee (IACUC) at Pennsylvania State University. The mice were exposed to a 12 h dark/light cycle with continuous access to water and food. VA-sufficient (VAS) and VA-deficient (VAD) mice were generated as described previously [23,25,46,47] by feeding VAD or VAS diets to pregnant dams and weanlings. The VAD diet contained no VA, whereas the VAS diet provided 25 µg of retinyl acetate per day. The diet composition table is provided in the supplementary Excel file [48]. After weaning, the mice were assigned to one of the four treatment groups (non-infected VAD, infected VAD, non-infected VAS, and infected VAS) until the end of the experiments. The VA status of the animals was confirmed by analyzing serum retinol using an ultra-performance liquid chromatography method [44].

Infection of C. rodentium
The origin, culture condition, and inoculation of C. rodentium have been described before [24]. Briefly, C. rodentium were grown overnight in a Difco Luria-Bertani (LB) broth containing 50 µg/mL nalidixic acid. The mice (8-10 weeks of age, individually housed) were fasted overnight prior to oral gavage with 5 × 10 9 CFU of C. rodentium. To quantify the bacterial burdens, fresh fecal samples were collected thrice a week, homogenized, and plated on LB agar plates containing nalidixic acid.

Tissue Collection and RNA Extraction
Two studies were conducted; one focused on the SI and one-on the colon, where the study was ended and tissues were collected on post-infection (p.i.) day 5 and day 10, respectively, the time of peak infection as shown in previous studies [23,25,38]. Details on tissue collection and RNA extraction were described previously [44]. Briefly, total RNA was extracted with a Qiagen RNeasy Midi Kit, genomics DNA was removed by a TURBO DNA-free kit, and RNA quality was determined on an Agilent Bioanalyzer to ensure RNAseq libraries were constructed with RNA samples of sufficient quality (RNA integrity number, RIN > 8).

RNAseq Library Preparation, Sequencing, and Mapping
Details on the methods of library preparation, sequencing, and mapping were provided before [44]. In brief, a TruSeq Stranded mRNA Library Prep kit and a HiSeq 2500 system (Illumina) were used to generate raw sequencing data. Quality trimming, mapping, coverage, and raw read counts were obtained with the following tools: trimmomatic [49], hisat2 alignment program [50], bedtools [51], and featureCounts [52].

Differential Expression
We removed low-expression transcripts from the raw count data matrix (Figure 1a). Then, differential expression (DE) analyses were performed using the DESeq2 package [53]. For the SI study, the transcript was removed if it had less than 10 counts in more than 8 of the 12 samples or if its total count was less than 200. For the colon study, the transcript was removed if it had less than 10 counts in more than 10 of the 16 samples or if its total count was less than 220 [44]. Differential expression analyses were performed and the 2 × 2 experimental design was resolved into three effects related to infection (Figure 1b): (I) the "VAS-Inf" effect, which corresponded to the comparison between the non-infected VAS verses the infected VAS groups; (II) the "VAD-Inf" effect, which corresponded to the comparison between the non-infected VAD and the infected VAD groups; and (III) "VA effect under infection", which compared the infected VAD and the infected VAS groups. Criteria for differentially expressed genes (DEGs) were set at adjusted p (padj) < 0.05 and |FoldChange| > 2. The stringent criteria of DEG may lead to a large number of unchanged genes, i.e., genes that are neither upregulated (Log2FoldChange > 1, padj < 0.05) nor downregulated DEG (Log2FoldChange < −1, padj < 0.05). Visualizations such as heatmaps and volcano plots were created via R packages pheatmap, ggplot, and ggrepel.
As a criterion for the DEGs corresponding to the interaction effect, padj < 0.05 was used, corresponding to the interaction effect. By definition, the interaction effect is the "difference on top of difference". In the current study, the interaction effect examined how the infection effects differed between the different VA statuses. In other words, the interaction effect reflected the differences between the "VAS-Inf" and "VAD-Inf" groups defined earlier (Figure 1b), which would suggest the mechanisms underlying the compromised host resistance to C. rodentium that has been demonstrated in VAD mice [23,25].

WGCNA
For better comparability and consistency, the same data matrix in differential expression was used to build co-expression networks using the weighted correlation network analysis (WGCNA) package in R (version 3.4.4, Vienna, Austria) [54]. In other words, the dataset went through the same preprocessing steps (screening and normalization) as the differential expression analysis (Figure 1a). The application of the WGCNA tool is consistent with our earlier work, where the approaches were detailed [44]. To obtain distinct modules with moderate sizes, we set the minimum height for merging modules at 0.25 and the minimum module size to 30. According to the WGCNA standard usage, modules are henceforth referred to by their color labels, e.g., Colon("Color"). Genes belonging to no specific modules were assigned to the Grey module. Heatmaps were depicted to visualize the expression pattern of each module using R package named pheatmap. Eigengenes are the first principal component of each set of module transcripts, describing most of the variance in the module gene expression [55]. As representations of individual modules, eigengenes were computed. Module-trait relationships were analyzed through correlating module eigengenes with the trait measurements, VA status, and infection status. The significance and correlation coefficient of each correlation between the traits and module eigengenes were visualized in a labeled heatmap, color-coded using red (positive correlation) and blue (negative correlation) shades.
Nutrients 2022, 14, 1563 5 of 26 enrichment background was set as the filtered gene list of the colon (15,340 genes) or the SI (14,368 genes) ( Figure 1a). As the statistical significance threshold for all functional enrichment analyses, padj (Benjamini and Hochberg adjusted for multiple comparisons) < 0.05 was used.

Functional Enrichment Analysis
Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichment were assessed via clusterProfiler (v3.6.0) in R (version 3.4.4, Vienna, Austria) through one-tailed Fisher's exact test, also known as the hypergeometric test [56]. The enrichment background was set as the filtered gene list of the colon (15,340 genes) or the SI (14,368 genes) ( Figure 1a). As the statistical significance threshold for all functional enrichment analyses, padj (Benjamini and Hochberg adjusted for multiple comparisons) < 0.05 was used.

Model Validation: Infection Status and VA Status
Infection with C. rodentium was confirmed in all the mice used for transcriptomic analysis by measuring fecal shedding specific to this bacterium thrice weekly. In agreement with our previous reports [23,25], shedding increased after oral inoculation and plateaued during the period of peak infection (p.i. day 10). Bacterial counts increased with time (p < 0.0001) but did not differ between the infected VAD and VAS mice during either the 5-day course of the SI study ( Figure S1a) or the 10-day course of the colon study ( Figure S1b). Therefore, based on these kinetics and our previous kinetic studies of the passage of fluorescently tagged C. rodentium through the mouse intestinal tract [23,57], p.i. day 5 was selected for collection of the SI and p.i. day 10-for collection of the colon for further analysis.
The VA status of the mice was confirmed by measuring retinol concentrations in serum, which were significantly higher in the VAS mice than in their VAD counterparts (p < 0.001 [44]). These findings for both fecal shedding curves and serum retinol concentrations validated that our animal model is appropriate for the transcriptomic analysis that followed.

Mapping and Overview of the Differential Expression Results in the SI and the Colon
The flow of analysis and mapping results are illustrated in Figure 1a. Mapping revealed 24,421 genes, which were subsequently filtered for low-expression genes and normalized using DESeq2. Among the 15,340 genes expressed in the colon, the number of DEGs (Figure 1b) was as follows: 4524 genes corresponding to the infection effect under the VAS status (VAS-Inf, list I in the supplementary Excel file); 4278 genes corresponding to the infection effect under the VAD status (VAD-Inf, list II in the supplementary Excel file); 959 genes corresponding to the VA effect under infection (list III in the supplementary Excel file); and 1329 genes corresponding to the interaction effect (list IV in the supplementary Excel file). Since no DEG was discovered for the VAS-Inf or interaction effect in the SI, the results in the colon will be the primary focus of this paper.

Functional Categories of Colon DEGs in the VAS Infected vs. Non-Infected States (VAS-Inf Effect)
Of the 4524 DEGs corresponding to the VAS-Inf effect in the colon (list I in the supplementary Excel file), 1424 were upregulated while 3100 were downregulated ( Figure S2 showing the volcano plot, the heat map, and enrichment categories). The genes belonging to the immune-related and epithelial cell (EC)-related categories are listed alphabetically by gene symbol in Table S1, grouped first by the genes upregulated in VAS-Inf, followed by the genes downregulated in VAS-Inf, each compared to VAS alone, and subdivided into the genes likely originating in immune cells and the genes likely from ECs: the first included caspases, chemokine ligands, interleukin receptors, matrix metallopeptidases, antimicrobial proteins, and tumor necrosis factor-related genes, while those related to ECs included villin, actin, claudins, ROS production, endoplasmic reticulum (ER) stress signals, and goblet cell-specific transcripts. Some examples of upregulated DEGs that are also in the category of the top 30 most abundant genes are as follows: actin beta (Actb), mucin 2 (Muc2), chloride channel accessory 3B (Clca3b), Clca4b, and anterior gradient 2 (Agr2). Functional enrichment analyses for the 1424 upregulated DEGs in the colon in the VAS-Inf conditions suggested significant enrichment in the categories primarily rele-vant to inflammation, host-pathogen interactions, innate and adaptive immune functions ( Figure S2c), as well as cell proliferation, apoptosis, ER stress, and ROS production ( Table 1). Several enrichment categories for biological pathways and processes previously shown to be induced by the C. rodentium infection were confirmed in the colon study (Table 1), adding to the validity of the model. Table 1. Enrichment categories for the biological pathways and processes induced by the C. rodentium infection confirmed in the colon study.

Biological Process
Ref.
Downregulated DEGs (Table S1) contained fewer genes that represent immune signals while including several categories that are detailed further in the Discussion, including the genes involved in RA metabolic processes (Table 2); the genes involved in transcellular calcium absorption (Table 3); and the genes for the transporters that affect the absorption of water, Cl − , and Na + (Table 4). Some most abundant downregulated DEGs were carbonic anhydrase 4 (Car4), carbonic anhydrase 2 (Car2), solute carrier family 26 member 3 (Slc26a3), and vitamin D (1,25-dihydroxyvitamin D3) receptor (Vdr). The GO and KEGG enrichment analyses for the

The Effect of Infection under VA Deficiency (VAD-Inf) and the Effect of the VA Status under Infection
We identified 4278 DEGs corresponding to the effect of infection under the VA deficiency conditions (VAD-Inf, list II in the supplementary Excel file), which were visualized ( Figure S3a), among which 984 upregulated DEGs were enriched for functions such as defense response to the bacterium, IL-17 signaling pathway, acute-phase response, keratinization, water homeostasis, and neutrophil apoptosis ( Figure S3c, blue box on the right). The 3294 downregulated DEGs were enriched for categories like ion channel activity, muscle contraction, enteric nervous system, chemotaxis, solute sodium symporter activity, leukocyte activation, and Th17 cell differentiation ( Figure S3d, blue box on the right). The enriched functional categories above were further compared with those derived from the VAS-Inf comparison (left pink boxes of Figure S3c,d). Whereas several GO and  Figure S3c,d, middle overlapping section in purple), it is clear that the VAS-Inf comparison, as compared to the VAD-Inf comparison, contained more biological functions in the upregulated category ( Figure S3c) and fewer in the downregulated category ( Figure S3d). These results showed that the C. rodentium infection interacted differentially with the VA status to alter the colon gene expression. The effect of the VA status under infection, reflected by 959 DEGs, was also identified (list III in the supplementary Excel file) and visualized ( Figure S3b).

The Effect of Infection under VA Deficiency (VAD-Inf) and the Effect of the VA Status under Infection
We identified 4278 DEGs corresponding to the effect of infection under the VA deficiency conditions (VAD-Inf, list II in the supplementary Excel file), which were visualized ( Figure S3a), among which 984 upregulated DEGs were enriched for functions such as defense response to the bacterium, IL-17 signaling pathway, acute-phase response, keratinization, water homeostasis, and neutrophil apoptosis ( Figure S3c, blue box on the right). The 3294 downregulated DEGs were enriched for categories like ion channel activity,

GO and KEGG Enrichment of DEGs: Four Interaction Scenarios Comparing Up-and Downregulated Gene Patterns in the Colon of VAS-Infected and VAD-Infected Mice
To further compare the gene expression patterns that distinguished the VAS-Inf group from the VAD-Inf group, we constructed four interaction "scenarios," i.e., four subgroups of DEGs that were distinct based on their expression patterns, for the reason that these patterns may provide clues to mechanisms of interaction between the nutritional status and infection. DEGs corresponding to each of the four interaction scenarios in the colon are visualized in Figure 3 and listed in the supplementary Excel file. Scenario 1 comprised the genes that were upregulated in the VAS-Inf comparison while being lower (six DEGs) or remaining unchanged (241 DEGs) in the VAD-Inf comparison; the higher expression of these genes was expected to be positively associated with host resistance and/or pathogen clearance in the VAS mice, while an absence of a response might suggest a reason why VAD hosts were less likely to effectively clear the pathogen and survive. Scenario 1 included 247 DEGs, with functional enrichment categories related to MHC class I, T cells, interferons, mitotic spindle, cytokine signaling, protein catabolism, apoptosis, and transcriptional regulation ( Figure S4l). Scenario 2 comprised the genes that were downregulated in the VAS-Inf comparison while being either higher (six DEGs) or unchanged (110 DEGs) in the VAD-Inf comparison; these transcripts may help to explain the greater severity of infection and a higher rate of mortality in the VAD mice. Scenario 2 included 116 DEGs (example genes see in Figure 4), but was, however, enriched for only two functional categories: "apical plasma membrane" and "apical part of cell," both of which are cellular component categories in the GO database and likely related to EC functions. Scenario 3 comprised the genes that were unchanged in the VAS-Inf comparison but upregulated in the VAD-Inf comparison; however, although there were 90 DEGs in Scenario 3 (example genes see in Figure S5), there was no enrichment in any functional category. Finally, Scenario 4 comprised the genes that were unchanged in the VAS-Inf comparison while being downregulated in the VAD-Inf comparison; these genes might be likely to be associated with host resistance and pathogen clearance, i.e., the downregulation of those genes under VAD conditions could be a reason for the observed higher rate of lethality. Scenario 4 included 134 DEGs and was enriched for functional categories relevant to immune functions such as antigen processing and presentation, natural killer cell-mediated immunity, as well as cellular processes of proliferation, differentiation, migration, activation of T cells, and leukocytes ( Figure 5).

Identifying WGCNA Modules That Were Significantly Correlated with Infection Effects and the Interaction of the VA Status and Infection in the Colon
The construction of co-expression gene networks via the WGCNA provides an independent systems biology perspective (Figure 1a). For this analysis, we set a soft threshold power β = 26; this identified 13 different modules in the colon. To determine if any of the modules of the co-expressed colon genes were correlated with the infection effect and the interaction effect, we tested the correlations between the traits (VA status, and infection status) and module eigengenes. Eigengenes are the first principal component of the given modules and may be considered as a representative of the gene expression profiles of the module. Six modules (Colon(Brown, Blue, Turquoise, Yellow, Purple, and Green-yellow)) were found to be significantly correlated with the infection status (p-values < 0.05, Figure 6a), among which the Colon (Blue, Yellow, Purple, and Green-yellow) modules were negatively correlated with the trait "infection status" (Figure 6a). Among those four modules, the Colon(Blue) module was the most significant (p-values = 9 × 10 −12 , correlation coefficient = −0.98). The correlation coefficient was negative, meaning the module contained genes with overall lower expression levels in the infected colons (Figure 6a). Not surprisingly, among the 5674 module members, more than half of the genes (n = 2967) were identified as the downregulated DEGs. The Colon(Blue) module mainly exhibited functional enrichment in neurological functions, transport, and extracellular matrix (Table S2), meaning those activities were reduced in the colon of infected mice, comparing with their non-infected counterparts. The Colon(Turquoise) module was positively correlated with the infection status and enriched for functional categories involved in mRNA processing, chromosome segregation, protein catabolic process, etc. (Table S2). The Colon(Brown) and Colon(Purple) modules were simultaneously associated with VA and the infection status (Figure 6b,c). Colon(Purple) exhibited enrichment of the "anion transmembrane transporter activity," whereas Colon(Brown) was not enriched for any functional categories (Table S2). For each module, a module member gene list (supplementary Excel file) and an eigengene bar graph ( Figure S6) are provided.      enrichment in neurological functions, transport, and extracellular matrix (Table S2), meaning those activities were reduced in the colon of infected mice, comparing with their non-infected counterparts. The Colon(Turquoise) module was positively correlated with the infection status and enriched for functional categories involved in mRNA processing, chromosome segregation, protein catabolic process, etc. (Table S2). The Colon(Brown) and Colon(Purple) modules were simultaneously associated with VA and the infection status (Figure 6b,c). Colon(Purple) exhibited enrichment of the "anion transmembrane transporter activity," whereas Colon(Brown) was not enriched for any functional categories (Table S2). For each module, a module member gene list (supplementary Excel file) and an eigengene bar graph ( Figure S6) are provided.

Approach and Validation
In this study, we used transcriptional profiling and bioinformatics in an unbiased, exploratory manner to better understand the ways in which the VA status and gut infection modulate the expression of genes in the SI and the colon, which may in turn suggest mechanisms of host defense and thus help to explain the observations that VA-deficient hosts are more likely than their VA-adequate counterparts to die from the C. rodentium infection [23,25]. An admitted limitation of the RNAseq analysis is that it only provides information on transcript levels, while other stages of regulation, such as protein expression and activity, are not detected. Nevertheless, RNAseq provides a comprehensive analysis that can be subjected to further mining using bioinformatics approaches. Moreover, RNAseq also does not rely on individual reference genes, an advantage over traditional PCR methods. Because there has been no previous comprehensive analysis of the gene expression patterns comparing the VAS and VAD nutritional conditions under both noninfected and infected states, the current study was designed to address these gaps. Firstly, we found that infection affected gene expression in the colon, but not the SI (which, however, was affected by the VA status, see [44]), and, therefore, our discussion herein is focused on the colon. As was illustrated in Figure 1a, the design used in this study allowed for several comparisons: the effect of the VA status (VAS vs. VAD) in the absence of infection, as previously reported [44], and the effects of infection in the VAS condition (VAS-Inf), which revealed several interesting patterns of DEGs that concern (i) VA metabolic genes, (ii) calcium transport and VD-related genes, and (iii) ion and water transport, and which, therefore, may be of interest in the relationship to diarrheal disease. We then discuss the main differences observed when comparing infection in the VAS and VAD conditions (interaction effect, which we divided into four scenarios). Finally, we note that the complete gene lists are provided as the supplementary Excel file, and the datasets are available at NCBI GEO (accession No. GSE143290), providing additional information for further analysis.
We first validated our model by measuring the serum retinol, which was higher in the VAS animals than in the VAD animals [44], and by monitoring fecal shedding, as in previous studies, to ensure the infection status ( Figure S1). As we began our RNAseq analysis, we looked for the DEGs that are known as C. rodentium-responsive in the colon to assure replication under our VAS-Inf conditions. Such genes included processes of epithelial hyperplasia and apoptosis [58], goblet cell depletion, as well as innate and adaptive immune responses [6], which were replicated in the present dataset (Table 1 and Figures S2 and S3). This is in line with the fact that Colon(Turquoise), a WGCNA module positively correlated with the infection status, was enriched for biological functions such as chromosome segregation, and protein catabolic process (Table S2). With respect to epithelial hyperplasia and apoptosis, a hallmark pathological feature of the C. rodentium infection is the transmissible murine crypt hyperplasia (TMCH) [58], defined by the thickening of colonic mucosa and caused by the excessive induction of epithelial repair and regeneration mechanisms. During peak infection, the expedited colonic hyperplasia and the apoptosis ( Table 1) can altogether serve protective roles for the host against C. rodentium because both processes can speed up the life cycle of the ECs, a proportion of which may have been attached by C. rodentium. When the rate of differentiation, migration, detachment, and programmed cell death can exceed that of the proliferation of the attached C. rodentium, gradual clearance of the pathogen will take place. Regarding goblet cell depletion, our data summarized in Table 1 agree with a previously published work that goblet cell activities are induced by and are protective against A/E pathogens, as increased mucin mRNA expression and mucus thickness during the C. rodentium infection were observed [18]. Mucin can displace the pathogen to the outer mucus layer and aid in flushing away the pathogen [73]. Some genes such as Il17a, Il22, and Ifng that would be expected to increase were not detected as DEGs in our study [74,75], perhaps due to the filtering and the expression-level criteria used. However, inspection of the raw data prior to filtering showed that there were changes in the expected direction in these genes. Moreover, the more abundant expressions of the receptors of those cytokines (Il-17 receptors a-e, IL-22 receptors a1 and a2, IFNg receptors Ifngr1 and Ifngr2), combined with other DEGs, helped to drive the enrichment categories "IL-17 signaling pathway," "Th17 cell differentiation," and "response to interferon-gamma," which can indirectly serve as a proof of the principle that cytokines were induced during peak infection.

DEGs Positively Associated with Host Resistance under the VAS Status
The analysis of DEGs expressed under VAS-Inf conditions also revealed upregulated genes that are likely attributable to immune cells in the lamina propria, such as genes for caspases, chemokine ligands, interleukin receptors, matrix metallopeptidases, antimicrobial proteins and tumor necrosis factor-related genes, etc., as well as genes likely derived from ECs, including villin, actin, claudins, ROS production, ER stress signals, and goblet cellspecific transcripts (Table S1), and the processes related to inflammation, host-pathogen interactions, innate and adaptive immune functions were increased ( Figure S2c), as well as cell proliferation, apoptosis, ER stress, and ROS production ( Figure S3c).
The host response to C. rodentium depends on the VA status, wherein the VAD mice showed a lower survival rate and a slower clearance rate [23,25]. To study the key pathways mediating the mitigated symptoms in the VAS hosts, the interaction effect was analyzed. Because the transcriptomic interaction between the VA status and the infection response could be either positive or negative, the analysis resulted in four different scenarios, of which Scenarios 1 and 4 were positively associated with a more robust host resistance. Scenario 1 contained the DEGs that were higher during infection in the VAS condition while being lower or unchanged with infection in the VAD condition. Therefore, this group of genes might be essential for the host response against the C. rodentium infection, and, conversely, the absence of their response might explain why the VAD hosts were unable to effectively clear the pathogen [23,25]. The Muc4 DEG belonged to Scenario 1. Muc4 mRNA was not influenced by the VA status (VAS vs. VAD) in a naïve mouse colon [44]. Herein, Muc4 was significantly upregulated during peak C. rodentium infection only in the VAS mice (VAS-Inf), but not in the VAD mice (VAD-Inf), suggesting it might play a protective role against the C. rodentium infection. The expression pattern of Muc4 ( Figure S4) is not only a good example as the interaction effect, but also suggests that VA signaling in the colon is a prerequisite for the induction of Muc4 during colitis.
C. rodentium-infected mice have been widely used to investigate innate and adaptive mucosal immunity. Lymphocytes, especially B cells and CD4 + T cells, are crucial for the resistance to the C. rodentium infection [66,67]. The cytokine IL-17 is required for the clearance of the pathogen [24,26,70,71]. Th17 cells are a distinct branch of mucosal immunity that deals with A/E lesions, and therefore can be induced by the C. rodentium infection [76]. Atarashi et al. [12] demonstrated that EC adhesion by C. rodentium and the subsequent ROS production mediated the induction of Th17 cells. In the current RNAseq dataset, upregulated DEGs nitric oxide synthase 2 (Nos2), dual oxidase 2 (Duox2), and dual oxidase maturation factor 2 (Duoxa2) ( Table 1 and Table S1), together with other DEGs, drove the enrichment of the "reactive oxygen species metabolic process" (Table 1 and Figure S3c). Furthermore, "Th17 cell differentiation" was a KEGG enrichment term for Scenario 4 of the interaction effect, that is, for the genes that were unchanged in VAS-Inf while being downregulated in VAD-Inf and which, therefore, might likely be associated with impaired host resistance and pathogen clearance under VAD conditions ( Figure 5a and Table 1).
Infection by several enteropathogens, including C. rodentium, leads to a drastic reduction in the number of the phenotypically distinct goblet cells, a process defined as the "goblet cell depletion" [58]. This may seemingly contradict the increased transcriptional activities in goblet cells during the C. rodentium infection (Muc, Retnlb, Clcas, and Fut2) observed in our study (Table S1), but a closer examination indicates that this evidence may indeed line up well with each other. Goblet cell depletion, reflected by a reduction of periodic acid-Schiff (PAS) staining on histological colon sections, was actually caused by a reduction of the mucin glycoprotein content of goblet cells rather than an actual lineage loss of goblet cells.
Whereas the majority of previous research was focused on the genes induced by C. rodentium, the downregulated DEGs were more than double the number of the upregulated DEG under the VAS-Inf conditions (3100 vs. 1424, Figure S2b and list I of the supplementary Excel file); those downregulated DEGs included categories of genes that may be most closely associated with the nutritional status and/or with diarrheal disease, as discussed next.

Three Categories of Downregulated DEGs Related to Nutritient Utilization and Ion Transport Relevant to Diarrheal Disease
Three groups of genes were identified from this study that could be of particular importance for understanding nutrient processes during infection and infection-related diarrhea. Firstly, our results showed that the C. rodentium infection concordantly downregulated VA metabolic genes in the VAS colon. Several genes involved in RA metabolic processes, which are known to be inducible by RA, were significantly downregulated in the colon of the VAS mice during peak infection (Table 2), including Aldh1a1, Aldh1a2, Cyp26b1, and Lrat [77,78]. Enzymes RALDH1 and RALDH2 (encoded by Aldh1a1 and Aldh1a2, respectively) convert retinal to RA [78], whereas cytochrome P450 enzymes (CYP26), present in a variety of tissues, catabolize RA to its less bioactive form and thus may be critical for preventing the buildup of high concentrations of RA intracellularly [79]. Similarly, lecithin retinol acyltransferase (LRAT) is also considered among the most important genes for regulating VA metabolism by esterifying retinol to its storage form, retinol ester. Previously, the expression levels of several of these genes were identified as being reduced in VAD vs. VAS SI [44]; however, here, differences were observed due to infection in the VAS host. Besides the reduction of these RA-responsive genes, a uniform downregulation was observed for other VA-metabolic genes (Bco2 (a beta-carotene cleavage enzyme that provides a precursor for RA synthesis [80]); Rarb, Rxrg (nuclear receptors mediating RA signaling), Rbp2 (binds retinol for its intracellular transport and conversion), Adh1, Adh7, Rdh5 (enzymes that catalyze the conversion from retinol to retinal), and Crabp1 (an intracellular RA-binding protein [81]), as well as upregulation of the two rate-limiting enzymes (Rdh1 and Dhrs9) for the production of retinal as a precursor of RA [82] (Table 2). These changes suggest that infection reduces the capacity to produce RA, while it may also interfere with VA storage and RA catabolism. Overall, the pattern of changes in the expression of genes suggests a reduced ability to utilize VA and/or produce RA for retinoid signaling. Although on a systemic level, the vast majority of VA absorption occurs higher up in the SI, a small amount may be directly absorbed by the colon, considering the recent findings that watersoluble vitamins, known to be uptaken in bulk in the SI, now found efficiently absorbed in colon [83]. Rodent colons do contain VA [84], and the colonic gene expression profile is influenced by the VA status [44], indicating there is either direct absorption in the colon and/or some VA is absorbed in SI then allocated to colon via blood circulation. Since VA needs to be converted to its bioactive form to take effect, the reduced local expression of RA metabolic genes in colon can compromise the health of colonocytes, as well as the education and differentiation of colonic immune cells (e.g., DCs, lymphocytes, and ILCs) [43]. Not necessarily expected, the C. rodentium infection in the colon recapitulated some of the effects previously noted for VA deficiency in the SI [44], suggesting there could be regulated VA metabolism, even in the colon, under the control of the VA status [44] or bacterial infection.
Secondly, we also observed that the C. rodentium infection concordantly downregulated the VD metabolic and calcium transport-related genes in the VAS colon. An interesting finding was the almost uniform downregulation of Vdr and the genes involved in calcium absorption during infection in the VAS-Inf group compared to the VAS control mice (Table 3). VDR is a nuclear receptor that induces the expressions of several genes in the transcellular calcium absorption process. We observed reductions during peak infection in Cacna1d (which encodes Ca v 1.3, a calcium channel allowing for Ca 2+ entry through the apical membrane of enterocytes), S100g (which encodes calbindin CB 9k , a binding protein that chaperones Ca 2+ to move from apical to basolateral membrane in enterocytes), Atp2b2, Atp2b3, Atp2b4 (which encode the plasma membrane Ca 2+ ATPase family that extrude Ca 2+ into the blood stream across the basolateral membrane), and Slc8a1 (which encodes sodium calcium exchanger, or NCX1, which also plays a role in the extrusion step) ( Table 3) [85]. From these results, it appears that infection causes an apparently coordinate reduction in processes necessary for calcium uptake.
In addition to being a regulator of the calcium absorption process, VD is also known as a regulator of mucosal immunity in the gut [86]. VD plays dual roles on activating/suppressing cells in the mucosal immune system, which seem to be dependent on the phase of the infection, where at the beginning of the infection VD is required for the proper initiation of the immune response, while the immunosuppressive role of VD emerges during the peak of infection and becomes essential during the resolution phase of the infection [87]. Macrophages, DCs, T cells, and B cells are all VD targets since they all express VDR [87]. Immune cells generally take 2-3 days after infection to become activated and maximize their VDR expression [87]. During the early phase of the C. rodentium infection, VD is essential for the differentiation of ILC3 and the production of IL-22, which may be required for the normal expansion of Th17 cells during the later phase of infection [88]. Previous studies using VD-deficient mice showed some similar changes, such as significantly fewer ILC3 cells and lower IL-22 production, failure to expand the Th17 population, clearance of the pathogens at a lower rate, more severe infection, and more rapid mortality comparing with their VD-sufficient counterparts [88]. Those results resemble our observation during VAS-Inf, when Vdr and the vast majority of VD target genes were downregulated (Table 3); meanwhile, enrichment terms for "interleukin-6 production," "TNF signaling pathway," "response to interferon-gamma," and "IL-17 signaling pathway" were increased ( Figure S2c). To sum up, these results suggest the idea that the C. rodentium infection may result in an overall reduction of the effective concentration of calcitriol in the colon, which may aid the host resistance/survival mechanisms through partially releasing the immunosuppressive effect of VD. Additionally, considering that IL-17 is required for the clearance of this A/E pathogen [24,26,70,71], reduced VA and vitamin D signaling locally in the colon under the state of immune activation status during peak C. rodentium infection might have been utilized by the host (either passively or actively) as an "on-switch" for the Th17 immune response. However, further experiments are necessary to test this possible connection.
Thirdly, our study provides evidence that C. rodentium infection concordantly downregulated ion transport pathways, which may underlie C. rodentium-induced diarrhea. Water loss during diarrhea is a major cause of mortality in experimental models and in humans [14,29,89,90]. Although the secretion of fluid and mucus are important mucosal defense mechanisms that can dilute and wash away pathogens and their products from the epithelial surface, diarrhea can be fatal due to the fluid losses, hypovolemia, and organ failure [91]. In the current study, downregulation of genes involved in Cl − absorption (Slc26a3), Cl − secretion (Cftr and Ano1), Na + absorption (Slc9a2 and Chp2), and HCO 3 − homeostasis (Car 2, 3, 4, 11, 14, 15 and Best2) were observed in VAS-Inf (Table 4), resulting in enriched functional categories of "sodium ion transport," "regulation of cytosolic calcium ion concentration," "regulation of cation transmembrane transport," "potassium ion transport," and "regulation of ion transmembrane transport" (Figure 2). In addition, Colon (Blue) and Colon (Purple) were two WGCNA modules suggesting lower transmembrane transport activities in infected colon (Table S2, Figure 4c, and Figure S6).
Regarding chloride transepithelial transport, DEGs were identified for both absorption and secretion-related processes. Slc26a3 encodes for the Downregulated in Adenoma (DRA), a Cl − /HCO 3 − exchanger, which mediates the apical chloride absorption into the ECs and the secretion of bicarbonate into the lumen. Mutation of DRA is associated with an intestinal disorder characterized by watery diarrhea, severe dehydration, high levels of fecal chloride, hypochloremia, and hyponatremia [92]. In previous studies, the downregulation of Slc26a3 have been accused for the fatal diarrhea, characterized by the loss of water in the stool [28]. The reduction of the Slc26a3 transcription may be resulted from the elevated TNFα signaling during peak C. rodentium infection [93]. Regarding chloride secretion, Cftr and Ano1, two transporters located on the apical membrane of colonic ECs that secrete chloride into the lumen [94] were found drastically reduced during VAS-Inf (Table 4), suggesting downregulation of Cl − loss through feces in VAS mice during diarrhea. This may have compensated the reduction of Cl − absorption, which may be essential for the survival of the VAS mice. Interestingly, Cftr and Ano1 were two genes belonging to Scenario 2 of the Interaction effect, which were not effectively downregulated in VAD-Inf (Figure 4a,c). This may suggest that excessive Cl − secretion into the lumen may have exacerbated the water loss through diarrhea.
With respect to sodium absorption, previous studies have showed suppression of Slc9a2, Slc9a3, and Chp2 expression by C. rodentium infection [28,29], wherein the reduction of Slc9a2 and Chp2 were replicated in the VAS-Inf comparison of the current study (Table 4). NHE2 (encoded by Slc9a2) and NHE3 (encoded by Slc9a3) are both major Na + /H + exchangers responsible for the apical sodium absorption, pH regulation, and fluid balance in the intestine [95,96]. Calcineurin-like EF hand protein 2 (CHP2), an essential cofactor for the NHE family members, supports the Na + /H + exchange activity [97]. The differential expression of these genes may also contribute to a disturbed ion balance.
With respect to bicarbonate metabolism, carbonic anhydrases (CAs) catalyze the reversible dehydration/hydration of CO 2 and water [98], thus maintaining the bicarbonate pool required for the ion-exchange activity of chloride and bicarbonate (e.g., DRA) and supplying protons for the apical membrane ion exchangers (e.g., NHEs and cHKA), and therefore, inhibition of carbonic anhydrases is associated with marked reductions in bicarbonate secretion as well as chloride, sodium, and water absorption [29,96,99]. Significant suppression of Car2 and Car4 by C. rodentium infection were observed previously [18,28,29] and in our study (Table 4). Car4 and Car14 encode for CA isoforms that are membranebound, whereas Car2, Car3, and Car13 encode cytosolic isoforms. Bestrophin-2 (Best2), a gene encoding a HCO 3 channel in the basolateral membrane, was also found to be down-regulated by the C. rodentium infection in previous works [18] and in our study under the VAS-Inf condition (Table 4). However, as a DEG belonging to Scenario 2 of the Interaction effect, it was not significantly downregulated in the VAD-Inf condition (Figure 4b). Slc4a7, encoding a Na + -HCO 3 − -cotransporter named NBCn1 in the basolateral membrane of the goblet cells in the colonic crypts [100,101], were found significantly upregulated during VAD-Inf, but was not significantly changed during the VAS-Inf (Scenario 3 of Interaction effect, Figure S5). HCO 3 − is of the utmost importance for the build-up of the mucus layer [102]. Since the upregulation of Slc4a7 was coincident with the increased mortality rate in the VAD mice, one possible explanation is that higher NBCn1 channel facilitated the HCO 3 − transportation from the blood stream into the colonocytes, which supports a higher mucus production. Even though extra mucus can help isolate, dilute, and flush away the pathogens, it also requires more water to hydrate the mucin molecules, resulting in more water loss and less water absorption in the VAD hosts. Dehydration, rather than sepsis, is the major cause of lethality in the C. rodentium infection within the susceptible mouse strains [29]. Therefore, our study supports an overall downregulation of transcriptional machinery to maintain the intracellular HCO 3 − level, and decreased chloride/sodium absorption during VAS-Inf, which is concordant with a mild diarrhea and high survival rate in the VAS mice. In contrast, the diminished compensatory effect reflected by the expression patterns of Cftr, Ano1, Best2, Slc4a7 and possibly other key genes in Scenario 2 and 3 may have exacerbated the water loss during the VAD-Inf, contributing to the higher mortality in the VAD hosts.

Conclusions
The RNAseq approach used in this study provides an abundance of data for these and further analyses. Finding that the response of VAD mice to infection differed from that of VAS mice with similar infection, as shown in terms of heat maps and several gene enrichment pathways, helps to generate new hypotheses for future experiments. Given that VAS mice have been shown to survive C. rodentium infection whereas a significant proportion of VAD mice succumb [23,25], it is of interest that several genes related to chloride, sodium, and bicarbonate ions, which responded in the VAS-Inf condition, did not respond in the VAD-Inf condition (Scenario 2). These observations lend support to the idea that that the host's response related to ion and water loss is crucial with respect to preventing mortality due to this gut infection, and that a deficiency of VA handicaps the host's ability to make an appropriate ion transport and fluid response.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/nu14081563/s1, Figure S1: C. rodentium clearance shows no significant difference between the infected VAD and VAS mice during the course of the SI and colon studies. Figure S2: Transcriptomic changes corresponding to VAS-Inf in colon study. Figure S3: Transcriptomic changes corresponding to VAD-Inf and the VA effect under infection in the colon study. Figure S4: Transcriptomic changes in Scenario 1. Figure S5: Representative DEGs corresponding to Scenario 3 in the colon study. Figure S6: Module eigengenes in the colon WGCNA network. Table S1: Immune-and epithelial-related DEGs upregulated and downregulated in VAS-Inf in the colon study.