Genetic architecture of mouse skin inflammation and tumour susceptibility.

Germline polymorphisms in model organisms and humans influence susceptibility to complex trait diseases such as inflammation and cancer. Mice of the Mus spretus species are resistant to tumour development, and crosses between M. spretus and susceptible Mus musculus strains have been used to map locations of genetic variants that contribute to skin cancer susceptibility. We have integrated germline polymorphisms with gene expression in normal skin from a M. musculus x M. spretus backcross to generate a network view of the gene expression architecture of mouse skin. Here we demonstrate how this approach identifies expression motifs that contribute to tissue organization and biological functions related to inflammation, haematopoiesis, cell cycle control and tumour susceptibility. Motifs associated with inflammation, epidermal barrier function and proliferation are differentially regulated in backcross mice susceptible or resistant to tumour development. The intestinal stem cell marker Lgr5 is identified as a candidate master regulator of the hair follicle, and the vitamin D receptor (Vdr) is linked to coordinated control of epidermal barrier function, inflammation and tumour susceptibility.

Skin tumors were induced on the dorsal back skin of 71 spretus / musculus backcross mice ([SPRET/Ei X FVB/N] X FVB/N, hereafter FVBBX mice, see Methods). We performed a combined correlation and linkage analysis of mRNA expression in uninvolved tail skin to characterize the genetic architecture of gene expression regulation in FVBBX mice. This approach identifies genetic loci that influence gene expression (eQTL), with either cis- or trans-acting effects. Several loci were hot-spots for significant trans-eQTL, including chromosome 2 at 85 Mb., chromosome 1 at 29 Mb., and chromosome 15 at 73 Mb (Supplementary Fig. 1). However, not all transcripts influenced by a given locus are functionally related to each other or co-regulated. We used a network analysis approach to identify loci that control functionally related groups of genes. We started by representing FVBBX tail skin RNA expression as a network where significantly correlated transcripts are drawn as nodes connected by an edge (see Methods). We then identified fully connected gene sets (cliques) that were enriched for functions of interest (Fig. 1). These sub-networks were highly enriched for genes representing either structural or functional components of the skin. These included skin-resident cell types (hair follicles, muscle, melanocytes, hematopoietic cells) and physiological (e.g. inflammation) or cellular (e.g. cell cycle) functions (Supplementary Fig. 4, Supplementary Table 3). We next identified the genetic loci associated with control of these different expression motifs (see Methods).

FIGURE 1

A visual representation of the FVBBX gene expression network

Black edges connect significantly correlated genes (blue nodes) that are co-regulated in cliques with five or more members. Labels identify clusters significantly enriched for functional or structural properties. Thousands of genes have significant eQTL, which are not shown here to preserve clarity.

A network comprising 62 genes, the majority of which are involved in hair follicle biology, was identified (Fig. 2, Supplementary Table 4). Numerous Type I and Type II hair follicle keratin genes such as Krt32 and Krt73 were present (Fig. 2A) along with genes involved in keratin structure and differentiation (Fig. 2C) and many keratin-associated proteins (Fig. 2D). Also identified were Msx2, Dlx3, and Lhx2 (Fig. 2B); these homeobox genes have been implicated in transcriptional control of hair follicle morphogenesis. Notably, epidermal keratins such as Krt5 and Krt14 are under separate genetic control and their expression was not correlated with hair follicle keratins. We extracted epithelial tissue from the tails of spretus, musculus, and SPRET/Ei X FVB/N mice (hereafter SPRET/Ei, FVB/N, and F1) to investigate the hypothesis that expression of hair follicle keratins is under complex genetic control. Expression of Type I and Type II hair follicle keratins was significantly different in these parental mice: they were approximately 8-fold higher in FVB/N mice than in SPRET/Ei. Expression of these genes was significantly higher in F1 mice than in either of the parental strains (Supplementary Fig. 2A), suggesting that expression of these structural protein genes is under complex genetic control. The transcription factors Dlx3 and Msx2 were also significantly differentially expressed in the same directions.

FIGURE 2

The hair follicle gene expression and linkage network

Grey edges connect significantly correlated genes (blue nodes) from several overlapping highly correlated cliques. Red and Green edges connect eQTL loci (red and green nodes) with significant linkage to genes. Green arrow-terminated edges indicate the gene is up-regulated by the SPRET/Ei allele; red bar-terminated edges indicate the SPRET/Ei down-regulates the gene. See also Supplementary Table 4. (A) Type I and Type II Keratins. (B) Transcription factors Msx2 and Dlx3 and stem cell marker Lgr5. (C) Genes associated with keratin biology & other genes. (D) Keratin-associated proteins.

The eQTLs in this hair follicle network controlling the largest number of genes were located on chromosome 10 at 118 Mb. (locus Chr10:118Mb) and on chromosome 3 at 19 Mb. (locus Chr3:19Mb). Inheritance of the SPRET/Ei allele at Chr3:19Mb is associated with higher expression levels of many genes in this network, while the SPRET/Ei allele on Chr10:118Mb decreases gene expression. Lgr5 (also known as Gpr49), a member of this network, is physically located on chromosome 10 at 115 Mb. and has a strong cis-eQTL at Chr10:118Mb (−Log10P = 4.3, q-value < 0.05, Fig. 2B). Lgr5 is a G-protein-coupled receptor recently identified as a stem cell marker in colon and small intestine[7]. Lgr5 is the only gene with a significant cis-eQTL in this sub-network (the other eQTL effects act in trans) and Lgr5 is plausibly the causative polymorphic gene on chromosome 10 that influences expression of these hair follicle genes. Lgr5 microarray results were validated by Taqman assay, and cis-regulation of Lgr5 expression in F1 mice was verified by a cis-trans test (Supplementary Fig. 5, 6). Sequencing of Lgr5 coding exons in SPRET/Ei and FVB/N identified four non-conservative differences between these two species (Supplementary Table 5). In agreement with this network view of the hair follicle, Lgr5 was recently identified as a hair follicle stem cell marker by lineage tracing[8]. Under normal homeostatic conditions, Lgr5 positive skin cells give rise to lineages of the hair follicle, but not to sebaceous glands or interfollicular epidermis[8].

Our approach also identified motifs associated with genetic control of hematopoietic cells and melanocytes (Fig. 3A, 3B, Supplementary Tables 6, 7). A hemoglobin production network including the hemoglobin isoform Hbb-b1, alpha-synuclein (Snca), and aminolevulinic acid synthase 2 (Alas2) was identified with eQTL that peak on chromosome 7 at 63 and 86 Mb. The observation of eQTL for expression of red blood cell-specific genes in the skin suggested that specific alleles control the number of skin-resident blood cells. This could be due to a physiological process that results in a specific homing of hematopoietic cells to the skins of particular animals, or may reflect increased systemic red blood cell production that is under genetic control. To test these alternative possibilities we generated a second FVBBX population (n = 83) and performed QTL analysis on Complete Blood Count (CBC) parameters (see Methods). Mean Corpuscular Volume showed a highly significant peak LOD score of 6.8 on chromosome 7 at 79 Mb and Mean Corpuscular Hemoglobin Content had a peak LOD score of 8.2 on chromosome 7 at 86 Mb (corrected p values < 0.001, Fig. 3C, 3D). This supports the hypothesis that genetic factors, rather than environmental factors such as skin inflammation, control blood cell phenotypes in the circulation and within the skin. It is also likely that variation in the numbers of melanocytes in the skin is under the control of genetic loci on chromosome 7 (Fig. 3B and Supplementary Table 7).

FIGURE 3

Hematopoiesis gene expression and linkage networks are confirmed by QTL results

Hematopoietic (A) and Melanosomal (B) eQTL networks in FVBBX epidermal RNA, drawn as in Fig. 2. See also Supplementary Tables 6, 7. Genome-wide plot of LOD scores from an 83 mouse FVBBX cohort for Mean Corpuscular Volume (C) and MCHC (D), where the horizontal line represents genome-wide significance.

We next identified networks associated with susceptibility to the development of skin tumors. The overall dorsal papilloma count in FVBBX mice after 20 weeks of DMBA / TPA treatment ranged from 0 to 35 (mean = 5.4, SD = 6.7). We assigned mice with zero papillomas at 20 weeks to a low susceptibility group and mice with eight or more papillomas to a high susceptibility group and performed differential expression analysis of tail epidermal RNA. 371 putative susceptibility genes were identified at a False Discovery Rate cut-off of 10%. This list was significantly enriched for epidermal keratinization and barrier function, the mitotic cycle and DNA replication, and interleukin 1 antagonist activity (Supplementary Table 8). Many genes involved in mitosis and cell proliferation (including Aurkb, Bub1, Cdca1, Cdca8, Ect2) were significantly co-regulated and associated with susceptibility, but do not share a statistically significant eQTL in normal skin. This motif includes Aurora-A, which has previously been identified as a genetic modifier of papilloma number[5].

We identified a locus on distal chromosome 15 (Chr15:101Mb) that influences expression of many of these genes (Fig. 4, Supplementary Table 9). This susceptibility network includes IL-1f6 and IL-1f5, ligands for the interleukin-1 receptor that can act positively or negatively in the inflammatory response, depending on the context[9]. The interleukin 1 receptor is the major mediator of inflammation in the skin[10]. Mice with a SPRET/Ei allele at Chr15:101Mb had higher expression of inflammation antagonists IL-1f5 and IL-1f6 and lower expression of pro-inflammatory genes such as Pde4b[. This pattern is compatible with the overall resistance of Mus spretus to inflammation noted by others[12]. We observed that higher expression of inflammation antagonists and lower expression of Pde4b was associated with high papilloma count, suggesting that in this context the spretus allele increases susceptibility, counter to the idea that inflammation is invariably associated with increased tumor susceptibility[13,14]. However, anti-inflammatory drugs can have contradictory effects on skin tumor development[15] and over-expression of pro-inflammatory cytokines such as IL-1 can prevent skin tumor formation in mouse models of chemically induced skin cancer[10]. The role of inflammation in skin cancer is therefore highly complex, with possibly different consequences associated with acute or chronic inflammatory conditions.

FIGURE 4

The Inflammation / Barrier Function gene expression and linkage network

Drawn as in Fig. 2. These genes are associated with susceptibility to papillomas, and many share a common significant eQTL at Chr. 15, 101 Mb. See also Supplementary table 9.

Several genes related to skin barrier function were highly correlated with the inflammation network, although genetic linkage to Chr15:101Mb was only detected at a lower level of stringency than that used for this analysis. These included Rptn and small proline-rich proteins Sprr2d and Sprr2i, major constituents of the epidermal barrier and cornified cell envelope. Disruption of skin barrier function leads to inflammation in epithelial tissue[16,17,18] but the exact mechanisms linking these processes are unclear. This network also includes the microbial peptide gene defensin beta 3 (Defb3), implicated in inflammation and psoriasis[19] and the anti-microbial calcium-binding gene S100a8, which complexes with S100a9 to form a pro-inflammatory, anti-bacterial cytokine present in keratinocytes that are inflamed or infected[20].

The best candidate master regulator of this network is the gene encoding the Vitamin D receptor (Vdr). This gene is physically located on chromosome 15 at 97.7 Mb. and shows a strong cis-eQTL (−Log10P = 5.8, q-value < 0.001). Sequencing Vdr identified several non-conservative sequence changes between SPRET/Ei and FVB/N in a phosphorylation domain and the hormone binding domain (Supplementary Fig. 7, 8). Vdr microarray results were validated by Taqman assay, and cis-regulation of Vdr in F1 mice was verified by a cis-trans test (Supplementary Fig. 5,6). Low levels of Vdr in FVBBX mice were associated with higher tumor susceptibility. Epidemiological studies have demonstrated an association between low serum levels of Vitamin D and human cancer susceptibility, emphasizing an important role for vitamin D in human cancer prevention[21]. The vitamin D signaling pathway has been linked directly to control of skin barrier function[22], inflammation and microbial defense[23], and pathways involved in stem cell growth[24]. Importantly, loss of the Vdr gene in knockout mice causes the phenotype that would be predicted from the network shown in Fig. 4, namely increased skin tumor susceptibility[25]. This suggests that Vdr is a master regulator of tissue damage and acute inflammation. These results emphasize the value of integrated systems analysis of polymorphisms and gene expression to identify groups of interacting genes and candidate master regulators that contribute to individual cancer susceptibility. Applying this approach to human samples may reveal the major genetic components of cancer susceptibility that remain to be discovered.

METHODS

Male SPRET/Ei and female FVB/N mice from the Jackson Laboratory were crossed. Female F1 hybrids were mated to male FVB/N mice. Seventy-one backcross mice (8-12 weeks old) received a dose of DMBA (25 μg per mouse in 200 μl acetone applied to shaved dorsal back skin). One week after initiation tumors were promoted with TPA (200 μl of 10−4 M solution in acetone) twice weekly for 20 weeks. Susceptibility groups were defined by papilloma count at 20 weeks (Low: zero at 20 weeks, n = 20, High: ≥ 8 at 20 weeks, n=17). Normal tail skin from FVBBX, SPRET/Ei, FVB/N and F1 mice (4 replicates per group) was snap frozen at sacrifice. mRNA expression profiles were generated with the Affymetrix M430 2.0 platform. Mice were genotyped using the ABI platform at 223 SNPs. Correlation analysis (Wilcoxon rank-sum) was performed with 24,357 transcripts expressed above background levels in FVBBX samples. Eighty-three additional FVBBX mice were generated and CBC was analyzed using a blood cell counter (HEMAVET; CDC Technologies, CT, USA). Linkage testing was performed by linear regression. An FDR method that accounts for signal dependence between loci on the same chromosome[26] identified 5% and 10% FDR p value cut-offs for linkage. CBC QTL were calculated and plotted using R/QTL[27], using 1000 permutations to calculate the 5% genome-wide error rate p-value cut-off. Gene Ontology analysis was performed with BiNGO[28]. The relevance network used an r2 cut-off of 0.64 ( < 0.01% alpha level, 1000 permutations). Cliques with at least five members were identified using a modification of the Bron-Kerbosch algorithm. Network figures were generated using Cytoscape version 2.5.1 (). Differentially expressed genes were identified using the Significance Analysis of Microarrays[29] with a maximum FDR of 10%.