- 1MRC Human Genetics Unit, Institute of Genetics and Cancer, University of Edinburgh, Edinburgh, United Kingdom
- 2Cancer Research UK Edinburgh Centre, Institute of Genetics and Cancer, University of Edinburgh, Edinburgh, United Kingdom
- 3Department of Public Health, Danish Institute for Advanced Study, University of Southern Denmark, Odense, Denmark
Colorectal cancer (CRC) is a common, multifactorial disease. While observational studies have identified an association between lower vitamin D and higher CRC risk, supplementation trials have been inconclusive and the mechanisms by which vitamin D may modulate CRC risk are not well understood. We sought to perform a weighted gene co-expression network analysis (WGCNA) to identify modules present after vitamin D supplementation (when plasma vitamin D level was sufficient) which were absent before supplementation, and then to identify influential genes in those modules. The transcriptome from normal rectal mucosa biopsies of 49 individuals free from CRC were assessed before and after 12 weeks of 3200IU/day vitamin D (Fultium-D3) supplementation using paired-end total RNAseq. While the effects on expression patterns following vitamin D supplementation were subtle, WGCNA identified highly correlated genes forming gene modules. Four of the 17 modules identified in the post-vitamin D network were not preserved in the pre-vitamin D network, shedding new light on the biochemical impact of supplementation. These modules were enriched for GO terms related to the immune system, hormone metabolism, cell growth and RNA metabolism. Across the four treatment-associated modules, 51 hub genes were identified, with enrichment of 40 different transcription factor motifs in promoter regions of those genes, including VDR:RXR. Six of the hub genes were nominally differentially expressed in studies of vitamin D effects on adult normal mucosa organoids: LCN2, HLA-C, AIF1L, PTPRU, PDE4B and IFI6. By taking a gene-correlation network approach, we have described vitamin D induced changes to gene modules in normal human rectal epithelium in vivo, the target tissue from which CRC develops.
Introduction
Colorectal cancer (CRC) is common, with over 40,000 incident cases and over 15,000 deaths associated with the disease per year in the United Kingdom (Cancer Research UK, 2019). In case-control and prospective cohort studies, higher plasma 25(OH)D level and higher dietary intake of vitamin D are associated with lower CRC risk (Jenab et al., 2010; Theodoratou et al., 2012; Theodoratou et al., 2014). Potential co-causality of CRC risk and vitamin D status (e.g. socioeconomic status, diet, physical exercise) and reverse causation (CRC or its treatment affecting serum vitamin D concentration) mean these observational studies may be confounded. Supplementation trials have been inconclusive, with randomised-controlled trials (RCTs) failing to show effect on CRC or adenoma (precursor lesion) incidence (Wactawski-Wende et al., 2006; Baron et al., 2015; Manson et al., 2019; Scragg, 2019). However, these supplementation studies are themselves confounded by short duration of follow up and vagaries in genetic, environmental, ethnic, dietary and ecological factors such as latitude, weather and sunlight exposure, with many participants in both control and intervention arms starting trials replete in vitamin D, and/or taking low-dose vitamin D supplementation. Whether vitamin D supplementation reduces risk of CRC therefore remains an open question. In addition, recent studies have suggested a beneficial effect for vitamin D supplementation on CRC mortality (Keum et al., 2019; Vaughan-Shaw et al., 2020).
The mechanisms by which vitamin D may modulate CRC risk and survival are not well understood (Alleyne et al., 2017). Potential mechanisms include induction of cell differentiation and apoptosis or inhibition of cell growth and proliferation (Lamprecht and Lipkin, 2003; Feldman et al., 2014). Understanding of vitamin D transcriptional responses in the colon and rectum has primarily come from studies in transformed cancer cell lines, which may not represent responses in normal healthy tissue (Mapes et al., 2014), hence it is advantageous and timely to prioritise human studies to explore potential mechanisms in the normal target tissue.
Gene co-expression networks are used to describe the pairwise relationships of a large number of gene expression variables (Zhang and Horvath, 2005). Genes related by high correlation coefficients are thought to be functionally related, members of the same pathway and/or co-regulated. Weighted gene co-expression network analysis (WGCNA) is one of the most widely used network methods, whereby gene correlations are raised by a soft-thresholding power to form a scale-free network. The advantage of this method is that all correlations are included for analysis. Methods that apply a hard-threshold (e.g. r > 0.7 being arbitrarily biologically relevant) lose connections below that threshold for further analysis (Zhang and Horvath, 2005). Highly connected genes within each module identified by WGCNA are defined as hub genes, with such genes thought to play key biological roles in that particular module or in regulation of a particular trait (Zhang and Horvath, 2005; Kogelman et al., 2014; Drag et al., 2017; Bakhtiarizadeh et al., 2018). By comparing differences in network structures (e.g. whether gene modules are preserved between conditions), it is possible to assess how groups of genes are perturbed by a certain condition (such as vitamin D supplementation) (Langfelder et al., 2011). This method has been successfully used to discover genes involved in endometriosis (Bakhtiarizadeh et al., 2018) and multiple cancer types (Sun et al., 2017; Zhang et al., 2018; Zhu et al., 2019).
In this study, we aimed to determine effects of vitamin D on normal rectal epithelium by assessing gene expression before and after 12 weeks of vitamin D supplementation in human subjects. We first assessed vitamin D effects on the expression of single genes, and then constructed a weighted correlation network to assess effects of vitamin D on groups of genes, and to identify hub genes in modules emerging after supplementation. We sought to functionally annotate genes and gene modules using Gene Ontology (GO) pathway analysis and to determine transcription factor binding sites common to hub genes in emergent modules. Finally, we sought to validate expression changes of hub genes in adult normal colorectal mucosa organoids treated with vitamin D. These organoid models are isolated from many of the vagaries in heterogeneity of the human population, yet maintain the genetic architecture and 3D-cell arrangement present in the parent tissue.
Materials and Methods
Human Vitamin D Supplementation Study
In the Scottish Vitamin D (SCOVID) study, rectal normal mucosa biopsies (via rigid sigmoidoscopy) and blood were collected after informed consent from a cohort of human subjects free from colorectal cancer (n = 50, 49 whose samples passed QC). Demographic information is provided in Table 1, while the study protocol has been described elsewhere (Vaughan-Shaw et al., 2021). The study had approval from NHS Research Ethics Committee (REC No 13/SS/0248) and local Research and Development Committee (R&D Project ID 2014/0058). Participants received 3200IU daily oral vitamin D (Fultium-D3), with resampling at 12 weeks. Biopsy samples were stored immediately in RNA Later (Invitrogen) and kept for 48 h at 4°C before RNA extraction.
Plasma 25(OH)D was assayed from blood by mass spectrometry. Plasma extracted from blood taken in lithium heparin tubes was immediately frozen at −40°C and subsequently submitted to the Clinical Biochemistry department, Glasgow Royal Infirmary, United Kingdom for measurement of 25(OH)D.
To extract RNA, human biopsy samples transferred to 2 ml Eppendorf tubes and homogenised in Trizol. RNA was then extracted by Ribopure Kit (Invitrogen) according to the manufacturer’s protocol. RNA samples from 49 subjects passed QC and were submitted to the Edinburgh Genomics facility, with sequencing on the Illumina HiSeq 2,500 in “rapid mode” with 150 bp paired-end reads as described in Supp Methods.
Analysis
Transcript quantification from RNAseq was conducted using Salmon v0.11 (Patro et al., 2017) using Ensembl version GRCh38, March 2017, Ensembl 88. Gene level counts were generated by R packages txiimport (Soneson et al., 2015) and annotated using biomaRt (Durinck et al., 2005; Durinck et al., 2009). Expression was normalised using a TMM algorithm based on gene expression thresholds of >0.1 Transcripts Per Million (TPM) and ≥6 reads in ≥20% of samples. Trimmed mean of M-values (TMM) between-sample normalisation was applied to the counts, as per the GTEx (v7) protocol. Effects of vitamin D on gene expression were assessed using edgeR v3.32.1 (Robinson et al., 2010). Following dispersion estimation, a quasi-likelihood negative binomial generalized log-linear model was fitted to the count data. Differential expression was assessed by quasi-likelihood F-test, with the parent-tissue anonymous identifier included as a covariate in the model as per a paired-design. Significance was determined as FDR corrected p value < 0.05. GO pathway analysis was carried out by clusterProfiler (Yu et al., 2012).
Normalised counts were first log transformed before proceeding to correlation analysis by WGCNA. To minimise the biological noise from genes not functionally related to vitamin D and to limit the dataset for computational analysis, the top 25% most variable genes after vitamin D supplementation (determined by logFC) were taken forward for analysis. WGCNA analysis is sensitive to the presence of outliers (Horvath, 2011), therefore samples with a standardized connectivity score of less than −5 were removed. The goodSamplesGenes function was then used to remove samples and genes with missing entries (more than 50% missing entries) and genes with zero variance.
Vitamin D supplementation achieved a 25(OH)D concentration which might be termed “healthy” whereas the pre-supplementation samples were depleted in 25(OH)D, and could be termed the “disease” state. Based on the assumption that modules present after vitamin D supplementation but not present before would be of interest in discerning vitamin D activity, the post-vitamin D samples were considered as the reference set for module derivation. Module preservation analysis was then undertaken in the pre-vitamin D samples. Using WGCNA, we created a signed weighted gene co-expression network based on normal gene expression data. A weighted network was created from the pairwise biweight midcorrelation coefficients between genes using the blockwiseModules function, with module merge cut height of 0.25 and a minimum module size of 30 genes (Song et al., 2012; Bakhtiarizadeh et al., 2018). A weighted adjacency matrix was formed by raising correlations to the power of 7, which was chosen using the scale-free topology criterion (Zhang and Horvath, 2005; Horvath, 2011). The relationship between the power (β) and R2 for a scale-free network is demonstrated along with sample dendrograms and their trait relationships in Supplementary Figure S1.
Age, gender and body mass index (BMI) have been reported to be associated with vitamin D concentration (Lagunova et al., 2009; Muscogiuri et al., 2019) and hence to assess association of those traits with identified gene modules, we assessed trait correlations to module eigenvectors (first principal component of each module). Each module eigenvector represents the expression profiles of all genes within that module.
To assess the preservation of post-vitamin D network modules in the pre-vitamin D dataset, the modulePreservation function in the WGCNA package was applied. We then applied Zsummary and medianRank (with 200 permutations) to detect module preservation. A module was considered as non-preserved if it had Zsummary<5 or medianRank≥8 (Langfelder et al., 2011).
To identify hub genes, intramodular connectivity (kIM) and module membership (kME) measures were used. Intramodular connectivity measures the degree of co-expression of a given gene with respect to the genes of a particular module. This was determined for both post- and pre-supplementation networks from the respective matrices of log transformed normalised counts using the intramodularConnectivity.fromExpr function (using pairwise biweight midcorrelation coefficients, power 7, a signed network and using the module labels identified above). Module membership was determined similarly by the signedKME function, which determines the correlation between the expression profile of a gene and the module eigengene (first principal component of a particular module). Genes were kME ≥0.7 or kIM ≥0.7 were considered as hub genes to the respective module (Horvath, 2011).
Validation of Gene Modules in the STRING Dataset
We sought to validate gene modules identified by WGCNA in the STRING (Search Tool for the Retrieval of Interacting Genes/Proteins) curated database of protein-protein interactions (https://string-db.org/). The STRING database collects, scores and integrates publicly available sources of protein–protein interaction information, and can be used to assess if groups of genes identified by the user are enriched for protein-protein interactions in those sources (Szklarczyk et al., 2019). HGNC gene names within each module were uploaded to the STRING web-browser, and interactions assessed across all seven of the STRING “interaction sources”.
Investigation of Common Transcription Factor Binding Sites
We were interested to assess if hub genes in non-preserved modules may be regulated by common transcription factors. As per Bakhtiarizadeh et al. (2018), the “TRANSFAC_and_JASPAR_PWMs” section of the Enrichr tool (Chen et al., 2013) (https://amp.pharm.mssm.edu/Enrichr/) was applied to determine common transcription factor binding sites in promoter regions of such genes.
Validation of Hub Gene Expression Changes in Adult Normal Colorectal Organoids
Fernandez-Barral et al. (2020) undertook differential expression analysis of adult normal mucosa organoids (FB-ANMO) derived from six individuals treated for 96 h with 100 nM calcitriol or 1% ethanol. Normalised counts and anonymised meta-data were downloaded from GEO (GSE100785; https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc = GSE100785). The same edgeR pipeline as described above was then used to determine differentially expressed genes in the FB-ANMO. Gene set testing was carried out using the geneSetTest function in limma (Ritchie et al., 2015).
Results
No Genes Were Differentially Expressed at a Genome-Wide Level Following Vitamin D Supplementation
Fultium-D3 supplementation increased plasma 25(OH)D in all study participants [mean (SD) nmol/l baseline 39.4 (20.4), 12 weeks 92.2 (27.0), p < 0.001; mean increase 52.8 (28.3)]. 21,650 genes passed expression filters and were taken forward for differential expression analysis. 2,492 genes were nominally differentially expressed (p < 0.05), but none remained significant after genome-wide FDR correction (FDR p < 0.05, Supplementary Table S1). 150 GO Biological Process terms were enriched in nominally differentially expressed genes including “protein-containing complex localization”, “DNA conformation change”, and “RNA splicing” (Supplementary Table S2).
Gene Modules Identified by WGCNA After Vitamin D Supplementation
Having identified no genes differentially expressed at a genome-wide significance level, we were interested to examine if vitamin D changes occurred across groups of genes by network connectivity analysis. From the total pool of 21,650 genes, the top 25% (5,412) genes by logFC (irrespective of direction of effect) were taken forward for further analysis. This limited analysis to the genes potentially most responsive to vitamin D supplementation for computational reasons. No samples had outlying standardized connectivity score, and no genes were lost by the goodSamplesGenes function.
17 modules were identified in the post-vitamin D network including the unclassified module (grey). No module eigenvector was correlated with plasma 25(OH)D or gender (Figure 1, Supplementary Table S3). The tan and greenyellow module eigenvectors were nominally significantly correlated with age, though did not remain significantly correlated after correction for multiple testing (r = 0.39, p = 0.005, FDR p = 0.09; r = -0.29, p = 0.04, FDR p = 0.38 respectively). The purple, salmon, red and turquoise module eigenvectors were nominally correlated with BMI, however again were not significant after accounting for multiple testing (r = 0.39, p = 0.005, FDR p = 0.09; r = -0.35, p = 0.01, FDR p = 0.13; r = -0.30, p = 0.04, FDR p = 0.16; r = -0.30, p = 0.04, FDR p = 0.16 respectively).
FIGURE 1. Module-trait relationships in the post-vitamin D network. Pearson correlation of module eigenvector (first principal component) and trait (along with nominal p-value).
Identification of Treatment-Associated Modules From the Pre-Vitamin D Transcriptome
One sample with outlying standardized connectivity score was removed from the pre-vitamin D network prior to module preservation analysis. 12 modules from the post-vitamin D network were strongly preserved in the pre-vitamin D network (defined as Zsummary statistic >10), i.e. were not associated with vitamin D treatment, while four modules from the post-vitamin D network were not preserved in the pre-vitamin D network, and hence were considered to be treatment-associated modules (Zsummary statistic <5 and median rank >8 for salmon, midnightblue and tan modules. Lightcyan module Zsummary statistic 5.2, Zdensity. pres 4.1 and Zconnectivity 6.3 with median rank 12). Module preservation statistics are described in Supplementary Table S4 and their distribution in Figure 2. GO terms for each module are described in Table 2 and Supplementary Table S5.
FIGURE 2. Median rank and Z-summary statistics for preservation of modules from the post-vitamin D network in the pre-vitamin D network. Z-statistic >10 strong evidence of preservation, 5–10 moderate evidence of preservation, 2–5 weak evidence of preservation and <2 no evidence of preservation (Langfelder et al., 2011).
TABLE 2. Select GO biological process terms in each of the modules present in the post-vitamin D network, with associated statistics for preservation in the pre-vitamin D network. The top three GO biological process terms for each module (determined by FDR p-value) are shown. Where highly similar terms related to overlapping genes existed (e.g., purine ribonucleoside binding and purine nucleoside binding) only one term is shown.
As a means of testing if WGCNA had identified gene modules which had also been found to be interacting in other datasets, we next sought to test if these modules were enriched in the STRING curated database of protein-protein interactions (PPI). Three of the four modules which were not preserved in the pre-vitamin D network were enriched for protein-protein interactions (Table 3). In addition 11 of the 12 preserved modules were enriched for protein-protein interactions (Supplementary Table S6).
TABLE 3. Protein-protein interaction enrichment of genes in non-preserved modules identified in the STRING database. FDR correction for 16 modules tested (grey unclassified module excluded).
Identification of Hub Genes in Treatment-Associated Modules
The genes with the highest degree of connectivity within a particular module are considered as hub genes. We were interested to identify hub genes in treatment-associated (non-preserved) modules, and in particular those genes which gained or lost hub-status. Hub genes were defined as those with kME ≥0.7 or with kIM≥0.7. By this definition, eight hub genes were identified in the salmon module, 13 in the midnightblue, 21 in the tan and nine in the lightcyan module (Supplementary Table S7).
Enrichment of Transcription Factor Binding Sites in Hub Gene Promoter Regions
Given hub genes are thought to be functionally important within their respective modules, we next sought to identify transcription factors common to hub genes in the non-preserved modules. A total of 40 different transcription factor motifs were enriched in promoter regions of hub genes in treatment-associated modules (Supplementary Table S8). Each of the four treatment-associated modules included hub genes which contained either VDR:RXR or RXR binding motifs in their respective promoter regions.
Crossover With Differentially Expressed Genes in NM Organoids
Finally, we reviewed differential expression of hub genes in a study of adult normal mucosa organoids treated with vitamin D. 15 of the 51 hub genes in treatment-associated modules identified above were present in the FB-ANMO dataset after gene filters. Those 15 genes were significantly downregulated on gene set testing in FB-ANMO (p = 0.02). Six of the 15 hub genes were also nominally differentially expressed in FB-ANMO, with five remaining significant after correction for multiple testing (Genome-wide FDR <0.05, Table 4, Supplementary Table S9). WGCNA is not recommended on datasets with fewer than 15 samples (Langfelder and Horvath, 2017), and hence we did not proceed to network analysis in this organoid dataset (6 samples per condition).
TABLE 4. Hub genes in non-preserved modules which are also differentially expressed in adult normal mucosa organoids from re-analysis of Fernandez-Barral et al. (FB-ANMO).
Discussion
In this study, we have demonstrated that while vitamin D supplementation does not affect the expression of any single gene at a genome-wide level of significance, it does induce changes in the inter-correlated expression patterns across genes, reflected in modules. 16 gene modules were identified after vitamin D supplementation (a state which reflects those individuals having a sufficient level of vitamin D for health). 14 of the 16 identified modules were also noted to be significantly enriched for protein-protein interactions in the STRING curated database, suggesting that WGCNA is able to identify modules which may have functional relevance.
Setting these results in context, other studies have similarly found no or few differentially expressed single genes following vitamin D supplementation (Hossein-nezhad et al., 2013; Munro, 2016; Pasing et al., 2017), however pathway analyses have identified effects on fatty acid metabolism and PPAR signaling (Munro, 2016), MAPK signaling, NF-kappa B signaling, T cell receptor signaling and prostate cancer (Hossein-nezhad et al., 2013). No study has previously investigated effects of vitamin D supplementation by the gene-correlation network approach. We were particularly interested to identify four treatment-associated modules that were not observed before vitamin D supplementation (a state of vitamin D insufficiency). Functional terms enriched in those modules included multiple GO terms related to the immune system, along with terms related to hormone metabolism, cell growth and RNA metabolism. Enrichment of GO terms associated with immunity is of particular note given the general interest in vitamin D and the immune system, in particular with risk of conditions such as multiple sclerosis and inflammatory bowel disease being associated with lower serum vitamin D (Munger et al., 2006; Limketkai et al., 2017; Fletcher et al., 2019). Vitamin D has also been postulated to favourably benefit the immune response in severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) (Martineau and Forouhi, 2020). Mechanistically, an excess of VDR binding variants identified by ChIP-exo has been reported to overlap with genomic variants associated with autoimmune disorders such as inflammatory bowel disease, Crohn’s disease and rheumatoid arthritis (Gallone et al., 2017).
Hub genes represent the most connected genes within a module, and are thought to be functionally important. Interestingly each of the non-preserved modules contained hub genes which contained VDR:RXR or RXR motifs in their respective promoter regions. Vitamin D signaling occurs principally following the binding of the active form of vitamin D, 1,25(OH)2D, to the Vitamin D Receptor (VDR) (Kliewer et al., 1992). VDR forms a heterodimer complex with Retinoid X Receptor (RXR) to bind DNA via VDR-responsive elements (VDRE) largely characterized by the VDR-RXR motif (Kliewer et al., 1992; Zhang et al., 2011). While the presence of a motif does not indicate actual transcription factor binding in-vivo, it does add further support to the notion that these genes are central to the vitamin D response in normal rectal epithelium.
Finally, when the list of hub-genes in non-preserved modules was cross-referenced with a separate study of vitamin D response in adult normal mucosa organoids, a number of interesting candidate genes were identified including LCN2, HLA-C and IFI6. LCN2, found to have a RXR motif in its promoter region, is expressed by macrophages and epithelia in response to inflammation (Dahl et al., 2018) and inhibition of LCN2-modulated NF-kB pathway activation by vitamin D has been noted to promote cisplatin sensitivity of oral squamous cell carcinomas (Huang et al., 2019). LCN2 acts in a bacteriostatic fashion (Dahl et al., 2018), which is noteworthy given the potential role of the gut microbiota in development of CRC (Saus et al., 2019). Human leukocyte antigens have been reported to be a major target of vitamin D physiological activity (Carlberg, 2019), with HLA-C being differentially expressed in peripheral blood mononuclear cells (PBMCs) following vitamin D supplementation of adult humans (Neme et al., 2019). Interferon alpha-inducible protein 6 (IFI6), also known as G1P3 has been shown to contribute to hyperplasia, tamoxifen resistance and poor outcomes in breast cancer (Cheriyath et al., 2012). It was one of the top differentially expressed genes (log FC -3.04) following vitamin D treatment of airway smooth muscle cells derived from individuals following a fatal asthma episode (Himes et al., 2015).
It is worthy of note that no module eigenvector was correlated with plasma 25(OH)D. The relationships between vitamin D and gene expression may not be linear, and instead sigmoidal or U-shaped relationships may exist (Mizunashi et al., 1995; Ross et al., 2011; Macdonald et al., 2018). Assessing the relationships between plasma 25(OH)D and gene modules may therefore fail to show effects in individuals sufficient in vitamin D, as was the case following vitamin D supplementation in this study. In addition, plasma vitamin D may not be an accurate measure of vitamin D in the target tissue of interest (in this case the rectal epithelium), with previous studies reporting a marked discrepancy between serum and tissue concentrations (Martinaityte et al., 2017).
This study represents a novel approach to assessing vitamin D effects. Unlike many of the large randomised-controlled trials of vitamin D effects (Wactawski-Wende et al., 2006; Baron et al., 2015; Manson et al., 2019; Scragg, 2019), the majority of participants were deficient in vitamin D at the start of the study period (n = 39 baseline plasma 25(OH)D < 50 nmol/l). The effects of vitamin D were also assessed directly in the target tissue of interest, the rectum, as opposed to assessing effects in blood which may have been technically easier to sample. This study is larger than many of the published studies assessing effects of vitamin D supplementation on gene expression (Hossein-nezhad and Holick, 2013; Gerke et al., 2014; Ryynänen et al., 2014; Protiva et al., 2016). Limitations of this study have been discussed elsewhere (Vaughan-Shaw et al., 2021). This study may have been too small to achieve sufficient power to assess individual gene significance. Individuals taking part in the study were not selected on the basis of initial plasma 25(OH)D; supplementation may have a sigmoidal or U-shaped relationship with gene expression (Mizunashi et al., 1995; Ross et al., 2011; Macdonald et al., 2018) and hence failing to select participants based on initial 25(OH)D could blunt the observed effect of supplementation. Finally sampling after 12 weeks of supplementation may not adequately capture early or later gene expression changes, however more frequent or delayed sampling would provide additional practical and ethical challenges.
Summary
By taking a gene-correlation network approach, we have described vitamin D-induced changes to groups of genes in normal human rectal epithelium. By reviewing treatment-associated modules before and after vitamin D supplementation, we have identified hub genes which may play a key role in modulating vitamin D actions in normal rectal epithelium. This provides novel understanding of the mechanisms by which vitamin D may have beneficial effects on CRC risk and survival.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, GSE157982.
Ethics Statement
The studies involving human participants were reviewed and approved by the NHS Research Ethics Committee (REC No 13/SS/0248). The patients/participants provided their written informed consent to participate in this study.
Author Contributions
Conceptualization, JPB, SMF, PVS, MGD; Methodology, JPB, PVS, KD, BH, VS, MT, TG, CS, FVND, SMF, MGD; Investigation, PVS, JPB, AMOF, PF, MT, MW, TG, SR, VS, KD; Writing - Original Draft JPB; Writing - Review & Editing, JPB, PVS, SMF, MGD, MT, KD, BH, FVND, CS, VS, AMOF, PF, MW, TG, SR; Funding Acquisition, SMF, MGD, PVS, JPB; Resources, SMF, MGD; Supervision, SMF, MGD.
Funding
This work was supported by funding for the infrastructure and staffing of the Edinburgh CRUK Cancer Research Centre; CRUK programme grant C348/A18927 (MGD/SMF). JB was supported by an ECAT-linked CRUK ECRC Clinical training award (C157/A23218). PV-S was supported by a NES SCREDS clinical lectureship, MRC Clinical Research Training Fellowship (MR/M004007/1), a Research Fellowship from the Harold Bridges bequest and by the Melville Trust for the Care and Cure of Cancer. The work received support from COST Action BM1206. FVND is supported by a CSO Senior Clinical Fellowship. This work was also funded by a grant to MGD as Programme Leader with the MRC Human Genetics Unit Centre Grant (U127527202 and U127527198 from 1/4/18).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Acknowledgments
We are grateful to Donna Markie and Fiona McIntosh, and all those who continue to contribute to recruitment, data collection, and data curation for the Scottish Vitamin D Study. We acknowledge the expert support on sample preparation from the Genetics Core of the Edinburgh Wellcome Trust Clinical Research Facility in addition to the nursing and study facilities provided by the Clinical Research Facility.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.783970/full#supplementary-material
References
A. C. Ross, C. L. Taylor, A. L. Yaktine, and H. B. Del Valle (Editors) (2011). Dietary Reference Intakes for Calcium and Vitamin D (Washington (DC): The National Academies Collection: Reports funded by National Institutes of Health).
Alleyne, D., Witonsky, D. B., Mapes, B., Nakagome, S., Sommars, M., Hong, E., et al. (2017). Colonic Transcriptional Response to 1α,25(OH) 2 Vitamin D 3 in African- and European-Americans. J. Steroid Biochem. Mol. Biol. 168, 49–59. doi:10.1016/j.jsbmb.2017.02.001
Bakhtiarizadeh, M. R., Hosseinpour, B., Shahhoseini, M., Korte, A., and Gifani, P. (2018). Weighted Gene Co-expression Network Analysis of Endometriosis and Identification of Functional Modules Associated with its Main Hallmarks. Front. Genet. 9, 453. doi:10.3389/fgene.2018.00453
Baron, J. A., Barry, E. L., Mott, L. A., Rees, J. R., Sandler, R. S., Snover, D. C., et al. (2015). A Trial of Calcium and Vitamin D for the Prevention of Colorectal Adenomas. N. Engl. J. Med. 373 (16), 1519–1530. doi:10.1056/nejmoa1500409
Cancer Research UK (2019). Colorectal Cancer Statistics 2019. Available from: https://www.cancerresearchuk.org/health-professional/cancer-statistics/statistics-by-cancer-type/bowel-cancer#heading-Three.
Carlberg, C. (2019). Vitamin D Signaling in the Context of Innate Immunity: Focus on Human Monocytes. Front. Immunol. 10, 2211. doi:10.3389/fimmu.2019.02211
Chen, E. Y., Tan, C. M., Kou, Y., Duan, Q., Wang, Z., Meirelles, G. V., et al. (2013). Enrichr: Interactive and Collaborative HTML5 Gene List Enrichment Analysis Tool. BMC Bioinformatics 14, 128. doi:10.1186/1471-2105-14-128
Cheriyath, V., Kuhns, M. A., Jacobs, B. S., Evangelista, P., Elson, P., Downs-Kelly, E., et al. (2012). G1P3, an Interferon- and Estrogen-Induced Survival Protein Contributes to Hyperplasia, Tamoxifen Resistance and Poor Outcomes in Breast Cancer. Oncogene 31 (17), 2222–2236. doi:10.1038/onc.2011.393
Dahl, S. L., Woodworth, J. S., Lerche, C. J., Cramer, E. P., Nielsen, P. R., Moser, C., et al. (2018). Lipocalin-2 Functions as Inhibitor of Innate Resistance to Mycobacterium tuberculosis. Front. Immunol. 9, 2717. doi:10.3389/fimmu.2018.02717
Drag, M., Skinkyté-Juskiené, R., Do, D. N., Kogelman, L. J. A., and Kadarmideen, H. N. (2017). Differential Expression and Co-expression Gene Networks Reveal Candidate Biomarkers of Boar Taint in Non-castrated Pigs. Sci. Rep. 7 (1), 12205. doi:10.1038/s41598-017-11928-0
Durinck, S., Moreau, Y., Kasprzyk, A., Davis, S., De Moor, B., Brazma, A., et al. (2005). BioMart and Bioconductor: a Powerful Link between Biological Databases and Microarray Data Analysis. Bioinformatics 21 (16), 3439–3440. doi:10.1093/bioinformatics/bti525
Durinck, S., Spellman, P. T., Birney, E., and Huber, W. (2009). Mapping Identifiers for the Integration of Genomic Datasets with the R/Bioconductor Package biomaRt. Nat. Protoc. 4 (8), 1184–1191. doi:10.1038/nprot.2009.97
Feldman, D., Krishnan, A. V., Swami, S., Giovannucci, E., and Feldman, B. J. (2014). The Role of Vitamin D in Reducing Cancer Risk and Progression. Nat. Rev. Cancer 14 (5), 342–357. doi:10.1038/nrc3691
Fernández-Barral, A., Costales-Carrera, A., Buira, S. P., Jung, P., Ferrer-Mayorga, G., Larriba, M. J., et al. (2020). Vitamin D Differentially Regulates colon Stem Cells in Patient-Derived normal and Tumor Organoids. FEBS J. 287 (1), 53–72. doi:10.1111/febs.14998
Fletcher, J., Cooper, S. C., Ghosh, S., and Hewison, M. (2019). The Role of Vitamin D in Inflammatory Bowel Disease: Mechanism to Management. Nutrients 11 (5), 1019. doi:10.3390/nu11051019
Gallone, G., Haerty, W., Disanto, G., Ramagopalan, S. V., Ponting, C. P., and Berlanga-Taylor, A. J. (2017). Identification of Genetic Variants Affecting Vitamin D Receptor Binding and Associations with Autoimmune Disease. Hum. Mol. Genet. 26 (11), 2164–2176. doi:10.1093/hmg/ddx092
Gerke, A. K., Pezzulo, A. A., Tang, F., Cavanaugh, J. E., Bair, T. B., Phillips, E., et al. (2014). Effects of Vitamin D Supplementation on Alveolar Macrophage Gene Expression: Preliminary Results of a Randomized, Controlled Trial. Multidiscip Respir. Med. 9 (1), 18. doi:10.1186/2049-6958-9-18
Himes, B. E., Koziol-White, C., Johnson, M., Nikolos, C., Jester, W., Klanderman, B., et al. (2015). Vitamin D Modulates Expression of the Airway Smooth Muscle Transcriptome in Fatal Asthma. PLoS One 10 (7), e0134057. doi:10.1371/journal.pone.0134057
Horvath, S. (2011). Weighted Network Analysis: Applications in Genomics and Systems Biology. Berlin: Springer Science & Business Media.
Hossein-nezhad, A., and Holick, M. F. (2013). Vitamin D for Health: a Global Perspective. Mayo Clinic Proc. 88 (7), 720–755. doi:10.1016/j.mayocp.2013.05.011
Hossein-nezhad, A., Spira, A., and Holick, M. F. (2013). Influence of Vitamin D Status and Vitamin D3 Supplementation on Genome Wide Expression of white Blood Cells: a Randomized Double-Blind Clinical Trial. PloS one 8 (3), e58725. doi:10.1371/journal.pone.0058725
Huang, Z., Zhang, Y., Li, H., Zhou, Y., Zhang, Q., Chen, R., et al. (2019). Vitamin D Promotes the Cisplatin Sensitivity of Oral Squamous Cell Carcinoma by Inhibiting LCN2-Modulated NF-Κb Pathway Activation through RPS3. Cell Death Dis 10 (12), 936. doi:10.1038/s41419-019-2177-x
Jenab, M., Bueno-de-Mesquita, H. B., Ferrari, P., van Duijnhoven, F. J. B., Norat, T., Pischon, T., et al. (2010). Association between Pre-diagnostic Circulating Vitamin D Concentration and Risk of Colorectal Cancer in European Populations:a Nested Case-Control Study. BMJ 340, b5500. doi:10.1136/bmj.b5500
Keum, N., Lee, D. H., Greenwood, D. C., Manson, J. E., and Giovannucci, E. (2019). Vitamin D Supplementation and Total Cancer Incidence and Mortality: a Meta-Analysis of Randomized Controlled Trials. Ann. Oncol. 30 (5), 733–743. doi:10.1093/annonc/mdz059
Kliewer, S. A., Umesono, K., Mangelsdorf, D. J., and Evans, R. M. (1992). Retinoid X Receptor Interacts with Nuclear Receptors in Retinoic Acid, Thyroid Hormone and Vitamin D3 Signalling. Nature 355 (6359), 446–449. doi:10.1038/355446a0
Kogelman, L. J. A., Cirera, S., Zhernakova, D. V., Fredholm, M., Franke, L., and Kadarmideen, H. N. (2014). Identification of Co-expression Gene Networks, Regulatory Genes and Pathways for Obesity Based on Adipose Tissue RNA Sequencing in a Porcine Model. BMC Med. Genomics 7, 57. doi:10.1186/1755-8794-7-57
Lagunova, Z., Porojnicu, A. C., Lindberg, F., Hexeberg, S., and Moan, J. (2009). The Dependency of Vitamin D Status on Body Mass index, Gender, Age and Season. Anticancer Res. 29 (9), 3713–3720. doi:10.14341/2071-8713-4886
Lamprecht, S. A., and Lipkin, M. (2003). Chemoprevention of colon Cancer by Calcium, Vitamin D and Folate: Molecular Mechanisms. Nat. Rev. Cancer 3 (8), 601–614. doi:10.1038/nrc1144
Langfelder, P., and Horvath, S. (2017). WGCNA Package FAQ. Available from: https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/faq.html.
Langfelder, P., Luo, R., Oldham, M. C., and Horvath, S. (2011). Is My Network Module Preserved and Reproducible. Plos Comput. Biol. 7 (1), e1001057. doi:10.1371/journal.pcbi.1001057
Limketkai, B. N., Mullin, G. E., Limsui, D., and Parian, A. M. (2017). Role of Vitamin D in Inflammatory Bowel Disease. Nutr. Clin. Pract. 32 (3), 337–345. doi:10.1177/0884533616674492
Macdonald, H. M., Reid, I. R., Gamble, G. D., Fraser, W. D., Tang, J. C., and Wood, A. D. (2018). 25-Hydroxyvitamin D Threshold for the Effects of Vitamin D Supplements on Bone Density: Secondary Analysis of a Randomized Controlled Trial. J. Bone Miner Res. 33 (8), 1464–1469. doi:10.1002/jbmr.3442
Manson, J. E., Cook, N. R., Lee, I.-M., Christen, W., Bassuk, S. S., Mora, S., et al. (2019). Vitamin D Supplements and Prevention of Cancer and Cardiovascular Disease. N. Engl. J. Med. 380 (1), 33–44. doi:10.1056/nejmoa1809944
Mapes, B., Chase, M., Hong, E., Ludvik, A., Ceryes, K., Huang, Y., et al. (2014). Ex Vivo culture of Primary Human Colonic Tissue for Studying Transcriptional Responses to 1α,25(OH)2 and 25(OH) Vitamin D. Physiol. genomics 46 (8), 302–308. doi:10.1152/physiolgenomics.00194.2013
Martinaityte, I., Kamycheva, E., Didriksen, A., Jakobsen, J., and Jorde, R. (2017). Vitamin D Stored in Fat Tissue during a 5-Year Intervention Affects Serum 25-Hydroxyvitamin D Levels the Following Year. J. Clin. Endocrinol. Metab. 102 (10), 3731–3738. doi:10.1210/jc.2017-01187
Martineau, A. R., and Forouhi, N. G. (2020). Vitamin D for COVID-19: a Case to Answer. Lancet Diabetes Endocrinol. 8 (9), 735–736. doi:10.1016/s2213-8587(20)30268-0
Mizunashi, K., Furukawa, Y., Abe, K., Takaya, K., and Yoshinaga, K. (1995). Resetting of Parathyroid Hormone Secretion after Vitamin D3 Treatment in Hypoparathyroidism and after Parathyroid Adenectomy in Primary Hyperparathyroidism. Calcif Tissue Int. 57 (1), 30–34. doi:10.1007/bf00298993
Munger, K. L., Levin, L. I., Hollis, B. W., Howard, N. S., and Ascherio, A. (2006). Serum 25-hydroxyvitamin D Levels and Risk of Multiple Sclerosis. JAMA 296 (23), 2832–2838. doi:10.1001/jama.296.23.2832
Munro, F. M. (2016). The Effect of Vitamin D on Gene Expression in Colorectal Tumours and normal colon. Otago: University of Otago. (Thesis, Master of Science) Retrieved from http://hdl.handle.net/10523/6170.
Muscogiuri, G., Barrea, L., Somma, C. D., Laudisio, D., Salzano, C., Pugliese, G., et al. (2019). Sex Differences of Vitamin D Status across BMI Classes: An Observational Prospective Cohort Study. Nutrients 11 (12), 3034. doi:10.3390/nu11123034
Neme, A., Seuter, S., Malinen, M., Nurmi, T., Tuomainen, T.-P., Virtanen, J. K., et al. (2019). In Vivo transcriptome Changes of Human white Blood Cells in Response to Vitamin D. J. Steroid Biochem. Mol. Biol. 188, 71–76. doi:10.1016/j.jsbmb.2018.11.019
Pasing, Y., Fenton, C. G., Jorde, R., and Paulssen, R. H. (2017). Changes in the Human Transcriptome upon Vitamin D Supplementation. J. Steroid Biochem. Mol. Biol. 173, 93–99. doi:10.1016/j.jsbmb.2017.03.016
Patro, R., Duggal, G., Love, M. I., Irizarry, R. A., and Kingsford, C. (2017). Salmon Provides Fast and Bias-Aware Quantification of Transcript Expression. Nat. Methods 14 (4), 417–419. doi:10.1038/nmeth.4197
Protiva, P., Pendyala, S., Nelson, C., Augenlicht, L. H., Lipkin, M., and Holt, P. R. (2016). Calcium and 1,25-dihydroxyvitamin D 3 Modulate Genes of Immune and Inflammatory Pathways in the Human colon: a Human Crossover Trial. Am. J. Clin. Nutr. 103 (5), 1224–1231. doi:10.3945/ajcn.114.105304
Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. Nucleic Acids Res. 43 (7), e47. doi:10.1093/nar/gkv007
Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2010). edgeR: a Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data. Bioinformatics 26 (1), 139–140. doi:10.1093/bioinformatics/btp616
Ryynänen, J., Neme, A., Tuomainen, T.-P., Virtanen, J. K., Voutilainen, S., Nurmi, T., et al. (2014). Changes in Vitamin D Target Gene Expression in Adipose Tissue Monitor the Vitamin D Response of Human Individuals. Mol. Nutr. Food Res. 58 (10), 2036–2045. doi:10.1002/mnfr.201400291
Saus, E., Iraola-Guzmán, S., Willis, J. R., Brunet-Vega, A., and Gabaldón, T. (2019). Microbiome and Colorectal Cancer: Roles in Carcinogenesis and Clinical Potential. Mol. Aspects Med. 69, 93–106. doi:10.1016/j.mam.2019.05.001
Scragg, R. K. R. (2019). Overview of Results from the Vitamin D Assessment (ViDA) Study. J. Endocrinol. Invest. 42 (12), 1391–1399. doi:10.1007/s40618-019-01056-z
Soneson, C., Love, M. I., and Robinson, M. D. (2015). Differential Analyses for RNA-Seq: Transcript-Level Estimates Improve Gene-Level Inferences. F1000Res 4, 1521. doi:10.12688/f1000research.7563.1
Song, L., Langfelder, P., and Horvath, S. (2012). Comparison of Co-expression Measures: Mutual Information, Correlation, and Model Based Indices. BMC Bioinformatics 13, 328. doi:10.1186/1471-2105-13-328
Sun, Q., Zhao, H., Zhang, C., Hu, T., Wu, J., Lin, X., et al. (2017). Gene Co-expression Network Reveals Shared Modules Predictive of Stage and Grade in Serous Ovarian Cancers. Oncotarget 8 (26), 42983–42996. doi:10.18632/oncotarget.17785
Szklarczyk, D., Gable, A. L., Lyon, D., Junge, A., Wyder, S., Huerta-Cepas, J., et al. (2019). STRING V11: Protein-Protein Association Networks with Increased Coverage, Supporting Functional Discovery in Genome-wide Experimental Datasets. Nucleic Acids Res. 47 (D1), D607–D613. doi:10.1093/nar/gky1131
Theodoratou, E., Palmer, T., Zgaga, L., Farrington, S. M., McKeigue, P., Din, F. V. N., et al. (2012). Instrumental Variable Estimation of the Causal Effect of Plasma 25-Hydroxy-Vitamin D on Colorectal Cancer Risk: a Mendelian Randomization Analysis. PloS one 7 (6), e37662. doi:10.1371/journal.pone.0037662
Theodoratou, E., Tzoulaki, I., Zgaga, L., and Ioannidis, J. P. A. (2014). Vitamin D and Multiple Health Outcomes: Umbrella Review of Systematic Reviews and Meta-Analyses of Observational Studies and Randomised Trials. BMJ 348, g2035. doi:10.1136/bmj.g2035
Vaughan-Shaw, P. G., Buijs, L. F., Blackmur, J. P., Theodoratou, E., Zgaga, L., Din, F. V. N., et al. (2020). The Effect of Vitamin D Supplementation on Survival in Patients with Colorectal Cancer: Systematic Review and Meta-Analysis of Randomised Controlled Trials. Br. J. Cancer 123 (11), 1705–1712. doi:10.1038/s41416-020-01060-8
Vaughan-Shaw, P. G., Grimes, G., Blackmur, J. P., Timofeeva, M., Walker, M., Ooi, L. Y., et al. (2021). Oral Vitamin D Supplementation Induces Transcriptomic Changes in Rectal Mucosa that Are Linked to Anti-tumour Effects. BMC Med. 19 (1), 174. doi:10.1186/s12916-021-02044-y
Wactawski-Wende, J., Kotchen, J. M., Anderson, G. L., Assaf, A. R., Brunner, R. L., O'Sullivan, M. J., et al. (2006). Calcium Plus Vitamin D Supplementation and the Risk of Colorectal Cancer. N. Engl. J. Med. 354 (7), 684–696. doi:10.1056/NEJMoa055222
Yu, G., Wang, L.-G., Han, Y., and He, Q.-Y. (2012). clusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters. OMICS: A J. Integr. Biol. 16 (5), 284–287. doi:10.1089/omi.2011.0118
Zhang, B., and Horvath, S. (2005). A General Framework for Weighted Gene Co-expression Network Analysis. Stat. Appl. Genet. Mol. Biol. 4, Article17. doi:10.2202/1544-6115.1128
Zhang, J., Chalmers, M. J., Stayrook, K. R., Burris, L. L., Wang, Y., Busby, S. A., et al. (2011). DNA Binding Alters Coactivator Interaction Surfaces of the Intact VDR-RXR Complex. Nat. Struct. Mol. Biol. 18 (5), 556–563. doi:10.1038/nsmb.2046
Zhang, X., Feng, H., Li, Z., Li, D., Liu, S., Huang, H., et al. (2018). Application of Weighted Gene Co-expression Network Analysis to Identify Key Modules and Hub Genes in Oral Squamous Cell Carcinoma Tumorigenesis. Ott 11, 6001–6021. doi:10.2147/ott.s171791
Keywords: colorectal cancer, vitamin D, vitamin D supplementation, gene correlation network, WGCNA (weighted gene co-expression network analysis)
Citation: Blackmur JP, Vaughan-Shaw PG, Donnelly K, Harris BT, Svinti V, Ochocka-Fox A-M, Freile P, Walker M, Gurran T, Reid S, Semple CA, Din FVN, Timofeeva M, Dunlop MG and Farrington SM (2022) Gene Co-Expression Network Analysis Identifies Vitamin D-Associated Gene Modules in Adult Normal Rectal Epithelium Following Supplementation. Front. Genet. 12:783970. doi: 10.3389/fgene.2021.783970
Received: 27 September 2021; Accepted: 02 December 2021;
Published: 04 January 2022.
Edited by:
Stephen J. Bush, University of Oxford, United KingdomReviewed by:
Wei Liu, Fujian Agriculture and Forestry University, ChinaMohammad Farhadian, University of Tabriz, Iran
Copyright © 2022 Blackmur, Vaughan-Shaw, Donnelly, Harris, Svinti, Ochocka-Fox, Freile, Walker, Gurran, Reid, Semple, Din, Timofeeva, Dunlop and Farrington. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: James P. Blackmur, jblackm2@ed.ac.uk; Susan M. Farrington, susan.farrington@ed.ac.uk