Gene Co-Expression Network Analysis Identifies Vitamin D-Associated Gene Modules in Adult Normal Rectal Epithelium Following Supplementation

Abstract

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 R² 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

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

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).

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).

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)₂D, 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.

References

• 1

  RossA. C.TaylorC. L.YaktineA. L.Del ValleH. B. (Editors) (2011). Dietary Reference Intakes for Calcium and Vitamin D (Washington (DC): The National Academies Collection: Reports funded by National Institutes of Health).

  • Google Scholar

• 2

  AlleyneD.WitonskyD. B.MapesB.NakagomeS.SommarsM.HongE.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. 10.1016/j.jsbmb.2017.02.001

  • CrossRef
  • Google Scholar

• 3

  BakhtiarizadehM. R.HosseinpourB.ShahhoseiniM.KorteA.GifaniP. (2018). Weighted Gene Co-expression Network Analysis of Endometriosis and Identification of Functional Modules Associated with its Main Hallmarks. Front. Genet.9, 453. 10.3389/fgene.2018.00453

  • CrossRef
  • Google Scholar

• 4

  BaronJ. A.BarryE. L.MottL. A.ReesJ. R.SandlerR. S.SnoverD. 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. 10.1056/nejmoa1500409

  • CrossRef
  • Google Scholar

• 5

  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.

  • Google Scholar

• 6

  CarlbergC. (2019). Vitamin D Signaling in the Context of Innate Immunity: Focus on Human Monocytes. Front. Immunol.10, 2211. 10.3389/fimmu.2019.02211

  • CrossRef
  • Google Scholar

• 7

  ChenE. Y.TanC. M.KouY.DuanQ.WangZ.MeirellesG. V.et al (2013). Enrichr: Interactive and Collaborative HTML5 Gene List Enrichment Analysis Tool. BMC Bioinformatics14, 128. 10.1186/1471-2105-14-128

  • CrossRef
  • Google Scholar

• 8

  CheriyathV.KuhnsM. A.JacobsB. S.EvangelistaP.ElsonP.Downs-KellyE.et al (2012). G1P3, an Interferon- and Estrogen-Induced Survival Protein Contributes to Hyperplasia, Tamoxifen Resistance and Poor Outcomes in Breast Cancer. Oncogene31 (17), 2222–2236. 10.1038/onc.2011.393

  • CrossRef
  • Google Scholar

• 9

  DahlS. L.WoodworthJ. S.LercheC. J.CramerE. P.NielsenP. R.MoserC.et al (2018). Lipocalin-2 Functions as Inhibitor of Innate Resistance to Mycobacterium tuberculosis. Front. Immunol.9, 2717. 10.3389/fimmu.2018.02717

  • CrossRef
  • Google Scholar

• 10

  DragM.Skinkyté-JuskienéR.DoD. N.KogelmanL. J. A.KadarmideenH. N. (2017). Differential Expression and Co-expression Gene Networks Reveal Candidate Biomarkers of Boar Taint in Non-castrated Pigs. Sci. Rep.7 (1), 12205. 10.1038/s41598-017-11928-0

  • CrossRef
  • Google Scholar

• 11

  DurinckS.MoreauY.KasprzykA.DavisS.De MoorB.BrazmaA.et al (2005). BioMart and Bioconductor: a Powerful Link between Biological Databases and Microarray Data Analysis. Bioinformatics21 (16), 3439–3440. 10.1093/bioinformatics/bti525

  • CrossRef
  • Google Scholar

• 12

  DurinckS.SpellmanP. T.BirneyE.HuberW. (2009). Mapping Identifiers for the Integration of Genomic Datasets with the R/Bioconductor Package biomaRt. Nat. Protoc.4 (8), 1184–1191. 10.1038/nprot.2009.97

  • CrossRef
  • Google Scholar

• 13

  FeldmanD.KrishnanA. V.SwamiS.GiovannucciE.FeldmanB. J. (2014). The Role of Vitamin D in Reducing Cancer Risk and Progression. Nat. Rev. Cancer14 (5), 342–357. 10.1038/nrc3691

  • CrossRef
  • Google Scholar

• 14

  Fernández-BarralA.Costales-CarreraA.BuiraS. P.JungP.Ferrer-MayorgaG.LarribaM. J.et al (2020). Vitamin D Differentially Regulates colon Stem Cells in Patient-Derived normal and Tumor Organoids. FEBS J.287 (1), 53–72. 10.1111/febs.14998

  • CrossRef
  • Google Scholar

• 15

  FletcherJ.CooperS. C.GhoshS.HewisonM. (2019). The Role of Vitamin D in Inflammatory Bowel Disease: Mechanism to Management. Nutrients11 (5), 1019. 10.3390/nu11051019

  • CrossRef
  • Google Scholar

• 16

  GalloneG.HaertyW.DisantoG.RamagopalanS. V.PontingC. P.Berlanga-TaylorA. J. (2017). Identification of Genetic Variants Affecting Vitamin D Receptor Binding and Associations with Autoimmune Disease. Hum. Mol. Genet.26 (11), 2164–2176. 10.1093/hmg/ddx092

  • CrossRef
  • Google Scholar

• 17

  GerkeA. K.PezzuloA. A.TangF.CavanaughJ. E.BairT. B.PhillipsE.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. 10.1186/2049-6958-9-18

  • CrossRef
  • Google Scholar

• 18

  HimesB. E.Koziol-WhiteC.JohnsonM.NikolosC.JesterW.KlandermanB.et al (2015). Vitamin D Modulates Expression of the Airway Smooth Muscle Transcriptome in Fatal Asthma. PLoS One10 (7), e0134057. 10.1371/journal.pone.0134057

  • CrossRef
  • Google Scholar

• 19

  HorvathS. (2011). Weighted Network Analysis: Applications in Genomics and Systems Biology. Berlin: Springer Science & Business Media.

  • Google Scholar

• 20

  Hossein-nezhadA.HolickM. F. (2013). Vitamin D for Health: a Global Perspective. Mayo Clinic Proc.88 (7), 720–755. 10.1016/j.mayocp.2013.05.011

  • CrossRef
  • Google Scholar

• 21

  Hossein-nezhadA.SpiraA.HolickM. 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 one8 (3), e58725. 10.1371/journal.pone.0058725

  • CrossRef
  • Google Scholar

• 22

  HuangZ.ZhangY.LiH.ZhouY.ZhangQ.ChenR.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 Dis10 (12), 936. 10.1038/s41419-019-2177-x

  • CrossRef
  • Google Scholar

• 23

  JenabM.Bueno-de-MesquitaH. B.FerrariP.van DuijnhovenF. J. B.NoratT.PischonT.et al (2010). Association between Pre-diagnostic Circulating Vitamin D Concentration and Risk of Colorectal Cancer in European Populations:a Nested Case-Control Study. BMJ340, b5500. 10.1136/bmj.b5500

  • CrossRef
  • Google Scholar

• 24

  KeumN.LeeD. H.GreenwoodD. C.MansonJ. E.GiovannucciE. (2019). Vitamin D Supplementation and Total Cancer Incidence and Mortality: a Meta-Analysis of Randomized Controlled Trials. Ann. Oncol.30 (5), 733–743. 10.1093/annonc/mdz059

  • CrossRef
  • Google Scholar

• 25

  KliewerS. A.UmesonoK.MangelsdorfD. J.EvansR. M. (1992). Retinoid X Receptor Interacts with Nuclear Receptors in Retinoic Acid, Thyroid Hormone and Vitamin D3 Signalling. Nature355 (6359), 446–449. 10.1038/355446a0

  • CrossRef
  • Google Scholar

• 26

  KogelmanL. J. A.CireraS.ZhernakovaD. V.FredholmM.FrankeL.KadarmideenH. 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. Genomics7, 57. 10.1186/1755-8794-7-57

  • CrossRef
  • Google Scholar

• 27

  LagunovaZ.PorojnicuA. C.LindbergF.HexebergS.MoanJ. (2009). The Dependency of Vitamin D Status on Body Mass index, Gender, Age and Season. Anticancer Res.29 (9), 3713–3720. 10.14341/2071-8713-4886

  • CrossRef
  • Google Scholar

• 28

  LamprechtS. A.LipkinM. (2003). Chemoprevention of colon Cancer by Calcium, Vitamin D and Folate: Molecular Mechanisms. Nat. Rev. Cancer3 (8), 601–614. 10.1038/nrc1144

  • CrossRef
  • Google Scholar

• 29

  LangfelderP.HorvathS. (2017). WGCNA Package FAQ. Available from: https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/faq.html.

  • Google Scholar

• 30

  LangfelderP.LuoR.OldhamM. C.HorvathS. (2011). Is My Network Module Preserved and Reproducible. Plos Comput. Biol.7 (1), e1001057. 10.1371/journal.pcbi.1001057

  • CrossRef
  • Google Scholar

• 31

  LimketkaiB. N.MullinG. E.LimsuiD.ParianA. M. (2017). Role of Vitamin D in Inflammatory Bowel Disease. Nutr. Clin. Pract.32 (3), 337–345. 10.1177/0884533616674492

  • CrossRef
  • Google Scholar

• 32

  MacdonaldH. M.ReidI. R.GambleG. D.FraserW. D.TangJ. C.WoodA. 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. 10.1002/jbmr.3442

  • CrossRef
  • Google Scholar

• 33

  MansonJ. E.CookN. R.LeeI.-M.ChristenW.BassukS. S.MoraS.et al (2019). Vitamin D Supplements and Prevention of Cancer and Cardiovascular Disease. N. Engl. J. Med.380 (1), 33–44. 10.1056/nejmoa1809944

  • CrossRef
  • Google Scholar

• 34

  MapesB.ChaseM.HongE.LudvikA.CeryesK.HuangY.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. genomics46 (8), 302–308. 10.1152/physiolgenomics.00194.2013

  • CrossRef
  • Google Scholar

• 35

  MartinaityteI.KamychevaE.DidriksenA.JakobsenJ.JordeR. (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. 10.1210/jc.2017-01187

  • CrossRef
  • Google Scholar

• 36

  MartineauA. R.ForouhiN. G. (2020). Vitamin D for COVID-19: a Case to Answer. Lancet Diabetes Endocrinol.8 (9), 735–736. 10.1016/s2213-8587(20)30268-0

  • CrossRef
  • Google Scholar

• 37

  MizunashiK.FurukawaY.AbeK.TakayaK.YoshinagaK. (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. 10.1007/bf00298993

  • CrossRef
  • Google Scholar

• 38

  MungerK. L.LevinL. I.HollisB. W.HowardN. S.AscherioA. (2006). Serum 25-hydroxyvitamin D Levels and Risk of Multiple Sclerosis. JAMA296 (23), 2832–2838. 10.1001/jama.296.23.2832

  • CrossRef
  • Google Scholar

• 39

  MunroF. 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.

  • Google Scholar

• 40

  MuscogiuriG.BarreaL.SommaC. D.LaudisioD.SalzanoC.PuglieseG.et al (2019). Sex Differences of Vitamin D Status across BMI Classes: An Observational Prospective Cohort Study. Nutrients11 (12), 3034. 10.3390/nu11123034

  • CrossRef
  • Google Scholar

• 41

  NemeA.SeuterS.MalinenM.NurmiT.TuomainenT.-P.VirtanenJ. 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. 10.1016/j.jsbmb.2018.11.019

  • CrossRef
  • Google Scholar

• 42

  PasingY.FentonC. G.JordeR.PaulssenR. H. (2017). Changes in the Human Transcriptome upon Vitamin D Supplementation. J. Steroid Biochem. Mol. Biol.173, 93–99. 10.1016/j.jsbmb.2017.03.016

  • CrossRef
  • Google Scholar

• 43

  PatroR.DuggalG.LoveM. I.IrizarryR. A.KingsfordC. (2017). Salmon Provides Fast and Bias-Aware Quantification of Transcript Expression. Nat. Methods14 (4), 417–419. 10.1038/nmeth.4197

  • CrossRef
  • Google Scholar

• 44

  ProtivaP.PendyalaS.NelsonC.AugenlichtL. H.LipkinM.HoltP. 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. 10.3945/ajcn.114.105304

  • CrossRef
  • Google Scholar

• 45

  RitchieM. E.PhipsonB.WuD.HuY.LawC. W.ShiW.et al (2015). Limma powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. Nucleic Acids Res.43 (7), e47. 10.1093/nar/gkv007

  • CrossRef
  • Google Scholar

• 46

  RobinsonM. D.McCarthyD. J.SmythG. K. (2010). edgeR: a Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data. Bioinformatics26 (1), 139–140. 10.1093/bioinformatics/btp616

  • CrossRef
  • Google Scholar

• 47

  RyynänenJ.NemeA.TuomainenT.-P.VirtanenJ. K.VoutilainenS.NurmiT.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. 10.1002/mnfr.201400291

  • CrossRef
  • Google Scholar

• 48

  SausE.Iraola-GuzmánS.WillisJ. R.Brunet-VegaA.GabaldónT. (2019). Microbiome and Colorectal Cancer: Roles in Carcinogenesis and Clinical Potential. Mol. Aspects Med.69, 93–106. 10.1016/j.mam.2019.05.001

  • CrossRef
  • Google Scholar

• 49

  ScraggR. K. R. (2019). Overview of Results from the Vitamin D Assessment (ViDA) Study. J. Endocrinol. Invest.42 (12), 1391–1399. 10.1007/s40618-019-01056-z

  • CrossRef
  • Google Scholar

• 50

  SonesonC.LoveM. I.RobinsonM. D. (2015). Differential Analyses for RNA-Seq: Transcript-Level Estimates Improve Gene-Level Inferences. F1000Res4, 1521. 10.12688/f1000research.7563.1

  • CrossRef
  • Google Scholar

• 51

  SongL.LangfelderP.HorvathS. (2012). Comparison of Co-expression Measures: Mutual Information, Correlation, and Model Based Indices. BMC Bioinformatics13, 328. 10.1186/1471-2105-13-328

  • CrossRef
  • Google Scholar

• 52

  SunQ.ZhaoH.ZhangC.HuT.WuJ.LinX.et al (2017). Gene Co-expression Network Reveals Shared Modules Predictive of Stage and Grade in Serous Ovarian Cancers. Oncotarget8 (26), 42983–42996. 10.18632/oncotarget.17785

  • CrossRef
  • Google Scholar

• 53

  SzklarczykD.GableA. L.LyonD.JungeA.WyderS.Huerta-CepasJ.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. 10.1093/nar/gky1131

  • CrossRef
  • Google Scholar

• 54

  TheodoratouE.PalmerT.ZgagaL.FarringtonS. M.McKeigueP.DinF. 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 one7 (6), e37662. 10.1371/journal.pone.0037662

  • CrossRef
  • Google Scholar

• 55

  TheodoratouE.TzoulakiI.ZgagaL.IoannidisJ. P. A. (2014). Vitamin D and Multiple Health Outcomes: Umbrella Review of Systematic Reviews and Meta-Analyses of Observational Studies and Randomised Trials. BMJ348, g2035. 10.1136/bmj.g2035

  • CrossRef
  • Google Scholar

• 56

  Vaughan-ShawP. G.BuijsL. F.BlackmurJ. P.TheodoratouE.ZgagaL.DinF. 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. Cancer123 (11), 1705–1712. 10.1038/s41416-020-01060-8

  • CrossRef
  • Google Scholar

• 57

  Vaughan-ShawP. G.GrimesG.BlackmurJ. P.TimofeevaM.WalkerM.OoiL. 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. 10.1186/s12916-021-02044-y

  • CrossRef
  • Google Scholar

• 58

  Wactawski-WendeJ.KotchenJ. M.AndersonG. L.AssafA. R.BrunnerR. L.O'SullivanM. J.et al (2006). Calcium Plus Vitamin D Supplementation and the Risk of Colorectal Cancer. N. Engl. J. Med.354 (7), 684–696. 10.1056/NEJMoa055222

  • CrossRef
  • Google Scholar

• 59

  YuG.WangL.-G.HanY.HeQ.-Y. (2012). clusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters. OMICS: A J. Integr. Biol.16 (5), 284–287. 10.1089/omi.2011.0118

  • CrossRef
  • Google Scholar

• 60

  ZhangB.HorvathS. (2005). A General Framework for Weighted Gene Co-expression Network Analysis. Stat. Appl. Genet. Mol. Biol.4, Article17. 10.2202/1544-6115.1128

  • CrossRef
  • Google Scholar

• 61

  ZhangJ.ChalmersM. J.StayrookK. R.BurrisL. L.WangY.BusbyS. A.et al (2011). DNA Binding Alters Coactivator Interaction Surfaces of the Intact VDR-RXR Complex. Nat. Struct. Mol. Biol.18 (5), 556–563. 10.1038/nsmb.2046

  • CrossRef
  • Google Scholar

• 62

  ZhangX.FengH.LiZ.LiD.LiuS.HuangH.et al (2018). Application of Weighted Gene Co-expression Network Analysis to Identify Key Modules and Hub Genes in Oral Squamous Cell Carcinoma Tumorigenesis. Ott11, 6001–6021. 10.2147/ott.s171791

  • CrossRef
  • Google Scholar

• 63

  ZhuZ.JinZ.DengY.WeiL.YuanX.ZhangM.et al (2019). Co-expression Network Analysis Identifies Four Hub Genes Associated with Prognosis in Soft Tissue Sarcoma. Front. Genet.10, 37. 10.3389/fgene.2019.00037

  • CrossRef
  • Google Scholar