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.

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

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

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.

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

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

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

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.

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

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

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.

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