RNAseq studies reveal distinct transcriptional response to vitamin A deficiency in small intestine versus colon, uncovering novel vitamin A-regulated genes

 Mbl2, Cxcl14, and Nr0b2 may be new candidate genes regulated by Vitamin A (VA) in small intestine..Mast cells and Tuft cells may be targets of VA in small intestinal..Small intestine and colon, although as consecutive organs, were impacted by VA distinctively

Vitamin A (VA) deficiency remains prevalent in resource limited areas

• Previous reports showed reduced pathogen clearance and survival due to vitamin A deficient (VAD) status
• To characterize the impact of preexisting VA deficiency on gene expression patterns in the intestines, and to discover novel target genes in VA-related biological pathways, VA deficiency in mice was induced by diet.

Vitamin A (VA), retinol and its metabolites

• In humans, a nutritional deficiency of VA is both a major cause of xerophthalmia and a risk factor for infectious diseases
• VA deficiency is estimated to affect nearly 20% of preschool aged children worldwide
• Several lines of evidence indicate that VA is a pleiotropic regulator of gut immunity
• Preexisting VA status has been shown to affect the immune response in several animal models of gut infection, including infection with Citrobacter rodentium (C. rodentium), a Gram-negative natural mouse intestinal pathogen that mimics enteropathogenic Escherichia coli infection in humans
• The initial VA status of mice at the time of infection could be an important factor in their ability to respond to gut infection

VAD and VAS mice were generated as previously described

• Briefly, VAS or VAD diets were fed to pregnant mothers and the weanlings.
• VAS diet contained 25 μg of retinyl acetate per day whereas VAD diet did not contain any VA. At weaning, mice were maintained on their respective diets until the end of the experiments.

Quantification of serum retinol by UPLC

• serum aliquots (30-100 μL) were incubated for 1 h in ethanol then 5% potassium hydroxide and 1% pyrogallol were added.
• After saponification in a 55°C water bath, hexanes (containing 0.1% butylated hydroxytoluene, BHT) and deionized water were added to each sample for phase separation.- The total lipid extract was transferred and mixed with a known amount of internal standard, trimethylmethoxyphenyl-retinol.

Statistical Analysis

• Two-way analysis of variance (ANOVA) with Bonferroni’s post-hoc tests was used to compare serum retinol levels.

Tissue collection and RNA extraction

• Tissues were collected at the end of each study. During tissue collection, the dissected colons were first cut open to remove fecal contents.
• Both lower SI and colon tissues were separated into segments and stored in individual tubes, and snap frozen in liquid nitrogen upon tissue harvest.

RNAseq library preparation, sequencing and mapping

• 200ng of RNA per sample were used as the input of TruSeq Stranded mRNA Library Prep kit
• The library concentration was determined by qPCR using the Library Quantification Kit Illumina Platforms
• An equimolar pool of the barcoded libraries was made and the pool was sequencing on a HiSeq 2500 system (Illumina) in Rapid Run mode using single ended 150 base-pair reads
• 30-37 million reads were acquired for each library
• Mapping were done in the Penn State Bioinformatics Consulting Center
• Quality trimming and adapter removal was performed using trimmomatic
• Hisat2 alignment was performed by specifying the 'rna-strandness' parameter
• Coverages were obtained using bedtools
• Reads mapped to each gene were counted using featureCounts

Data pre-processing

• Transcripts with low expression levels were first removed to ensure that each analysis would focus on tissue-specific mRNAs
• SI study: if the expression level of a gene was lower than 10 in more than eight samples, or if the sum of the expression in all 12 samples was higher than 200, this transcript was removed from the later analysis
• Colon study: low expression level in 10 samples or higher than 220, the transcript was regarded as lowly expressed gene and removed

Differential expression

• Analysis was performed using the DESeq2 package
• adjusted P<.05 to ensure statistical significance
• |Fold change|> 2 (to ensure biological significance) as the criteria for differentially expressed genes (DEGs).
• To visualize the result and prioritize genes of interest, heatmaps and volcano plots were depicted using the R packages pheatmap, ggplot, and ggrepel.

WGCNA

• The same data matrix in differential expression was used to build signed co-expression networks
• All 16 samples were taken into consideration during the construction of the SI and colon studies
• For each set of genes, a pair-wise correlation matrix was computed and an adjacency matrix was calculated by raising the correlation matrix to a power
• To obtain moderately large and distinct modules, we set the minimum module size to 30 genes and the minimum height for merging modules at 0.25 meters
• Eigengenes were computed and used as a representation of individual modules

Functional enrichment analyses

• Gene ontology (GO) and KEGG enrichment were assessed using clusterProfiler (v3.6.9.0) in R

Cell type marker enrichment

• The significance of cell type enrichment was assessed for each module in a given network with a cell type marker list of the same organ, using a one-sided hypergeometric test (as we were interested in testing for over-representation).
• To account for multiple comparisons, a Bonferroni correction was applied on the basis of the number of cell types in the organ-specific cell type list.

Module visualization

• igraph was used to visualize the co-expression network of genes in each individual modules based on the adjacency matrix computed by WGCNA
• The connectivity of a gene is defined as the sum of connection strength (adjacency) with the other module members
• In this way, the connectivity measures how correlated the gene is with all other genes in the same module.

Module preservation

• Genes with low expression (based on the criteria stated in Differential Expression session) in the SI or Colon datasets were excluded from this analysis.
• Reconstruction of co-expression networks were done in both datasets separately, with minimum module size = 30 and the minimum height for merging modules = 0.25
• All the modules from the analysis were marked as mpSI(Color) or mpColon(Color), and color assignment was remapped so that the corresponding modules in the re-constructed mpSI network sharing the similar expression pattern were assigned with the same color.

Module overlap

• The overlap between all possible pairs of modules was calculated along with the probability of observing such overlaps by chance.
• To quantify module comparison, a Bonferroni correction was applied on the basis of the number of modules in the comparison network (mpColon network). A module pair with corrected hypergeometric P value <.05 was considered a significant overlap.

RAR/RXR binding site prediction in promoter regions of genes of interest

• CiiiDER was used to predict transcription factor binding sites (TFBSs) across the promoter regions.

Model validation: Serum retinol is lower in VAD mice than VAS mice

• The VA status of the animal model was validated by determining serum retinal concentrations at the end of the SI study.
• Body weights of VAD compared to VAS animals were significantly lower in both groups.

Distribution of DEGs in SI and colon and their co-expression

• Whole tissue RNAseq analysis yielded 24,421 genes, of which 14,059 genes were shared between the two datasets and were used for the module preservation analysis.
• Among the 14,368 SI-expressed genes, 49 DEGs differed with VA status, and 9 co- expression modules were identified (Fig. 2).
• Similarly, among the 15,340 genes mapped in the colon dataset, 94 DEG differed with CA status and 13 co-expressions were identified.

Identification of DEGs

• Of the 49 DEGs in the SI that differed significantly by VA status, 18 were upregulated and 31 were downregulated (higher in VAS than VAD, Table S1)
• A volcano plot of the SI DEGs (Fig. 3A) shows that genes with the largest negative or positive fold changes as well as highest statistical significance included Cd207, Mcpt4, and Fcgr2b.
• Several genes with known regulatory functions were among the most strongly regulated in our datasets.

GO and KEGG enrichment of DEGs

• 49 SI DEGs were significantly enriched in "retinoid metabolic process" (GO term) and downregulated DEGs in VAS vs VAD mice.
• 18 of the 31 downregulated SI deGs (n=34) had a higher expression of Bco1, a beta-carotene cleavage enzyme compared to VAS (Fig. 3C).
• In the colon, downregulated VAS compared to upregulated VAD DEGs (N=60) had significantly enriched levels of cell division and cell cycle.

WGCNA modules significantly correlated with VA status

• Analysis suggests that innate immune activity is stronger, while adaptive immune response is weaker in the VAS SI compared to the VAD group
• To visualize the co-expression relationships within the modules of interest, the top 50 genes selected by connectivity are depicted as networks in Figure 5.5.

Associations with cell markers

• To query whether the expression pattern of each module was mainly driven by a regulatory effect of VA on cell differentiation, a cell type marker enrichment was performed.
• The SI(Turquoise) module, the largest module with MS=5,187 genes (Table S5), was significantly enriched for markers of goblet cells (CP=1 × 10−6) and Paneth cells.

Module preservation analysis suggests that gene networks controlled by VA are poorly preserved between SI and colon

• To determine whether the modules identified as being significantly correlated with VA in SI have a similar module structure in the colon, we performed module preservation analysis using 14,059 genes that were robustly expressed in both datasets.
• The modules from this analysis are referred to as mpSI(Color) and mpColon(Color).
• For all three of the mpSI modules that we found, we found only “weak to moderate” evidence of preservation in the MPColon network according to Zsummary.

RAR/RXR binding sites were predicted in the promoter regions of genes of interest

• Here are the genes that we focused on:
• Scarb1, Bco1, Isx, Cd207, Mbl2, Clca1, Mmp9, and Nr0b2.

Differential expression analysis confirmed known targets under the control of VA in SI and validated the nutritional model

• Of the 49 DEGs differing by VA status in SI, there was a strong enrichment in "retinoid metabolic process" driven by three DEGs: Cyp26b1, Isx, and Bco1
• RA induces the expression of ISX, a transcriptional repressor for both BCO1, which catalyzes ß-carotene cleavage, and SCARB1, a receptor for carotene uptake, in the intestine
• A series of immunological genes previously reported to be regulated by VA, albeit in a variety of tissues, were also identified in the SI dataset, including Mmp9, Cd207, and enzymes produced by mast cells
• Mast cell proteases (MCPT) 1, 2, 4, and 6 (also known as TPSB2 or tryptase beta 2) were uniformly downregulated in the VAS group

Potential new target genes regulated by VA in SI

• Mannose-binding lectin (MBL) is an important pattern-recognition receptor in the innate immune system
• SI is a predominant site of extrahepatic expression of MBL
• Allelic variants of Mbl2 have been associated with deficiencies in MBL protein production, which may increase the risk of mother-to-child HIV transmission
• VA plus β-carotene supplementation partially counteracted the increased risk of transmission

RA treatment study and in silico prediction validated VA targets in SI

• Transcriptomes of lower SI of mice that were either VA deficient or VA deficient treated four times with an oral dose RA, the bioactive form of VA, were compared.
• With the stringent cutoff for DEG (adjusted P<.05 and |Fold change|> 2), the discoveries of the Scarb1, Bco1, Isx, Cd207, and Mbl2 corresponding to the VA effect in our SI study were replicated (Fig. 7A, Tables S8 and S9).
• The observed consistency between the genes discovered in the situation of acute RA treatment and with long-term VA status, taken together with the fact that RAR/RXR binding sites were predicted in the promoter regions of those genes, suggests that those DEGs we have identified could be direct targets under the transcriptional regulation of RA in SI.

VA as a potential regulator of the SI Tuft cell population

• Tuft cells normally represent only about 0.5% of the epithelial cells in the murine small and large intestines, with a slightly more prevalent distribution in the distal SI.
• One of the advantages of the co-expression network analysis used in our study is that cell type-specific information may be derived through the analysis of whole intact tissues, without the isolation of homogeneous populations of cells.

VA differentially regulates SI and colon transcriptomes

• DEGs in the SI were fewer in number but were condensed in more biological activities than those in the colon
• The primary transcriptional role of VA in colon is related to inhibition of cell proliferation and induction of epithelial differentiation, two widely applicable functions of RA
• Colon was enriched for fewer categories that differed by VA status
• Although most genes were commonly expressed in both organs, the overlap of VA-correlated modules between the two organs were likely to be insignificant

Limitations and future directions

• The current study focused on SI-only but lack the colon counterpart, and is potentially involved with a discrepancy between the effect of acute oral RA supplementation and the long-term accumulative effect of VA status
• In the future, larger sample sizes should be considered to increase power
• Single cell techniques may be employed to tease out the cell-type specific phenomena

Conclusion

• DEGs corresponding to VA effect are present in both the SI and colon. In SI, DEGs altered by VA are mainly involved in the retinoid metabolic pathway and immunity-related pathways.
• With regard to the colon, VA appears to play an overall suppressive role on cell division while promoting cell differentiation, consistent with known functions of VA in other organ systems.

Reference

10.1016/j.jnutbio.2021.108814