Dissecting genetics of cutaneous miRNA in a mouse model of an autoimmune blistering disease

MicroRNAs (miRNAs) are small endogenous non-coding RNAs that control genes at post-transcriptional level. They are essential for development and tissue differentiation, and such altered miRNA expression patterns are linked to the pathogenesis of inflammation and cancer. There is evidence that miRNA expression is genetically controlled similar to the transcription of protein-coding genes and previous studies identified quantitative trait loci (QTL) for miRNA expression in the liver. So far, little attention has been paid to miRNA expression in the skin. Moreover, epistatic control of miRNA expression remains unknown. In this study, we characterize genetic regulation of cutaneous miRNA and their correlation with skin inflammation using a previously established murine autoimmune-prone advanced intercross line. We identified in silico 42 eQTL controlling the expression of 38 cutaneous miRNAs and furthermore found two chromosomal hot-spots on chromosomes 2 and 8 that control the expression of multiple miRNAs. Moreover, for 8 miRNAs an interacting effect from pairs of SNPs was observed. Combining the constraints on genes from the statistical interaction of their loci and further using curated protein interaction networks, the number of candidate genes for association of miRNAs was reduced to a set of several genes. A cluster analysis identified miR-379 and miR-223 to be associated with EBA severity/onset, where miR-379 was observed to be associated to loci on chromosome 6. The murine advanced intercross line allowed us to identify the genetic loci regulating multiple miRNA in skin. The recurrence of trans-eQTL and epistasis suggest that cutaneous miRNAs are regulated by yet an unexplored complex gene networks. Further, using co-expression analysis of miRNA expression levels we showed that multiple miRNA contribute to multiple pathways that might be involved in pathogenesis of autoimmune skin blistering disease. Specifically, we provide evidence that miRNA such as miR-223 and miR-379 may play critical role in disease progression and severity.


Background
The discovery of microRNAs (miRNAs), a class of small endogenous non-coding molecules ranging from 18-24 nt brought a new level of complexity for understanding the mechanisms that constitute various biological processes [1,2]. The involvement of miRNAs in the control of gene expression has been thoroughly defined for the cell cycle, metabolism, and immune system and in cancer [3][4][5][6]. Binding to the 3′ or 5′ UTR region of genes, miRNAs may yield increased or decreased gene expression levels and have been described to affect various molecular pathways [7].
Approximately half the miRNAs are intergenic with few also located in intronic regions [8,9]. These are understood to have their own enhancers and promoters and are transcribed by RNA polymerase II [10]. However, it remains unclear whether these are produced as by-products of protein-coding gene transcription or whether their biogenesis has its own machinery [11]. Recent studies in murine and human fibroblasts of liver tissues revealed that miRNA expression can either be regulated by their transcriptional genomic location or by other regions in the genome [12]. Furthermore, mutations in genes involved in miRNA processing, such as AGO1, DGCR8, and DICER, can cause significant changes in the expression of miRNAs, resulting in altered disease susceptibility [13]. Regarding the skin, in vivo studies show that miRNA biogenesis is dependent on both DICER and DGCR8. A lack of these enzymes causes severe phenotypes [14], underscoring the importance of miRNAs in the regulation of morphogenesis and homeostasis of the skin [15]. Differentially expressed miRNAs are associated with different physiological and pathological processes in the skin such as melanoma, Sézary syndrome, psoriasis and atopic dermatitis [16][17][18][19]. Advancements have been achieved in understanding various processes regulated by miRNAs. On the other hand, the mechanisms that regulate their own expression have so far remained unknown. One possible approach towards understanding their regulation is constraining on genetic loci whose variations statistically link the variable phenotypic effect. We thus performed an in silico expression QTL-based analysis, which has been widely used for deciphering genetic loci regulating gene expression in various biological processes [20]. Su et al. have further shown that expression QTL-based analysis can provide insights into miRNA regulation [12]. However, no study has yet been performed to determine the impact of genetic variation on miRNAs associated to skin and autoimmune disorders.
This study (i) uses expression QTL-based analysis to define genetic loci that control miRNA expression in the skin and (ii) provides insights into interacting genes to control a single transcript (epistasis). The relevance of these findings is illustrated by employing a mouse model for the autoimmune blistering disease EBA (Epidermolysis Bullosa Acquisita), providing a direct link between miRNA expression and disease phenotype. Taken together, the data provide insights into the complexity of miRNA regulation and possible means to understand various gene interactions altering the expression of miRNAs.

Results
A murine heterogeneous inter-cross line was generated by inter-crossing four parental inbred mouse strains: MRL/MpJ, NZM2410/J and BXD2/TyJ as autoimmunityprone strains and Cast/EiJ for genetic heterogeneity. As a mouse model for an autoimmune blistering disease, 100 mice from generation 4 of the intercross-line were immunized with recombinant collagen type VII (COL7), an integral component of anchoring fibrils located at the dermal-epidermal junction. This induced a loss of tolerance and production of anti-COL7 autoantibodies in all immunized mice. The incidence of sub-epidermal blisters and the clinical phenotype of epidermolysis bullosa acquisita (EBA) were at 1/3 [21]. All mice were genotyped and miRNA expression profiling was acquired from skin tissue using the Affymetrix's GeneChip miRNA 2.0 array. To reveal genetic loci that regulate miRNA expression, we then performed an association study between the genotype and the respective expression levels of miRNAs. Further we performed co-expression analysis between expression levels of miRNAs and phenotype. The workflow is presented in a flowchart (Additional file 1: Figure S1).

miRNA expression is genetically controlled
We performed a genome-wide scan to detect genetic loci associated with miRNA expression. Expression levels of miRNA were treated as quantitative trait to yield expression quantitative trait loci (eQTL). Genome-wide significance was determined by a permutation test. Using an E-value cutoff of <0.05, 42 eQTL for 38 miRNAs were mapped to the genome, corresponding to 6.83 % of all murine miRNAs present on the above mentioned Affymetrix GeneChip (Table 1 and Fig. 1). Since the wild derived strain CAST/EiJ was incorporated into the advanced intercross line, we investigated the polymorphic sites located in the transcribed miRNAs that may have an effect of probe hybridization which is derived from C57B6/J [22]. We obtained the SNPs and indels in genome of CAST/EiJ from the database [23]. We found that 4/38 (10.53 %) miRNAs (miR-291a-3p, miR-341, miR-449b and miR-681) exhibits indels in CAST/EiJ strain on chromosome 7 (3.2 Mb), 12 (69 and 109 Mb) and 13 (113 Mb), which may suggest false positive associations for those loci ( Table 1). The highest -log P value of 6.57 was observed for miR-298 on chromosome 9 explaining 20.76 % of the variance. The peak SNP (rs3700596) was found within 1 kb of the Ube3D gene, an ubiquitinconjugating enzyme E2c binding protein. We found only trans -eQTL, except for miR-486 (−log P = 4.10, rs13479880), which was mapped to the same chromosome of its transcriptional site (Chr 8,~89 Mb). This indicates that other genes rather than their own transcript may contribute to the regulation of miRNA expression which could be tissue specific.
Prior knowledge suggests classes of gene families that play a major role in the regulation of miRNA expression i.e. Argonaute proteins and helicases which play a crucial role in its regulation [24]. To examine this assumption, we obtained the coordinates of all helicases and other genes involved in pathways of miRNA biogenesis from 'The RNA Helicase Database' and mapped them to the eQTL identified in our study [25]. Helicases like Ddx39, Ddx49, Cd97 and Upf1 were mapped within the Other genes which do not belong to these classes, such as lin28a and its homolog lin28b, have been shown to modulate let-7a [26]. These genes were mapped to the eQTL for miR-290 on chromosome 4 and miR-126 on chromosome 10, suggesting that Lin28 might modulate the expression of other miRNAs as well.
Some miRNA eQTL were constrained to a specific location in the genome, thereby suggesting potential eQTL hotspots for the regulation of multiple miRNAs. On chromosome 2, five miRNAs (miR-26a, mir-291a, miR-423, miR-671 and miR-23b) were mapped between 28-51 Mb (Fig. 2). Three nearby SNPs (rs6250599~48.5 Mb, rs6265423~47 Mb and rs13476472~45 Mb) showed significant association (−log P > 4) with all five miRNAs that were mapped to this region. A long non coding RNA 1700019E08Rik was located near SNP rs13476472 (~3 Kb upstream). As for the other two SNPs, the nearest gene to rs6265423 was mapped 6.8 kb apart, coding for the snRNA U7.39-201, while the nearest coding gene for SNP Similarly, on chromosome 8 we identified eQTL for three miRNAs (miR-486, miR487b and miR-501) in confidence interval 72-95 Mb. Three genes Gm1068, Siah1 and Dnaja2 were located near peak SNPs (rs13479880 and rs6257357).
To search the candidate interacting gene pairs we investigated the epistatic control of miR-501-3p. Using the Ingenuity pathway analysis (IPA), we searched for all the possible interacting genes present between loci present on chromosome 1 and chromosome 2 [27]. We found three putative interacting gene pairs; namely Commd3 with Cops5, Cacnb2 with Vopip1 and Commd3-Bmi1 with Rb1cc. Interestingly, Cops5 in the Nfkb1 pathway possibly be associated with miR-501 via Tp53 (Additional file 1: Figure S2).

Genetic overlap of QTL for clinical scores and miRNA expression
Any genetic variation with an effect on the clinical phenotype and miRNA levels should map to the same chromosomal region. Previously, genetic loci for the clinical phenotype of EBA, an autoimmune skin blistering disease,  [21]. In this study, we found 4 eQTL for miRNAs overlapping with the QTL for EBA (Fig. 4). The eQTL for miR130b (Chr 9: 36-44 Mb, −log P = 4.42), miR-542-3p (Chr 12: 79-103 Mb, −log P = 4.73) and Chr Chromosome, Pos (Mb) Position in million base pairs, Pval(log) Additive log 10 P value, I.Pval(log) Interacting log 10 P value Fig. 4 Overlapping QTL for EBA disease and miRNA eQTL. The circular plot shows all the eQTL for miRNAs (green) and QTL for EBA (red). It also presents the eQTL hot spots (dark green) and EBA QTL for onset (red) and severity (dark red). Each circular band represents a chromosome on which QTL and eQTL are mapped. The region within the chromosome which has either red or dark red and green or dark green bands is overlapping eQTL with EBA QTL miR-449b (Chr 14: 50-72 Mb, −log P = 4.49) were mapped on QTL for disease onset on chromosome 9, 12 and 14. Additionally, we mapped eQTL for miR-493 (50-60 Mb, −log P = 4.98) to the QTL for both disease severity and onset on chromosome 19. To further evaluate, if the variation in the genome is associated with the miRNA expression levels and also underlies the clinical phenotype, we investigated the genotype association of significant SNPs of miRNA eQTL with clinical phenotype (Fig. 5). We observed that the most significant SNP (rs6211533) mapped for the eQTL of miR-493 on chromosome 19 also shows variation for EBA severity score and onset. In this eQTL minor allele (BB) derived from MRL/MpJ was associated with higher expression of miR-493, late onset of the disease, and higher clinical score (Fig. 5b). Similar observation was also derived for eQTL mapped (peak SNP, rs6413270) for miR-130b where minor allele (BB) derived for MRL/MpJ and NZM/2410J was associated with lower expression of miR-130b, and late onset of disease phenotype (Fig. 5a).

Expression and co-expression of miRNAs
Immunization of mice with recombinant COL7 leads to development of subepidermal blisters and the clinical phenotype of EBA in 1/3 of the immunized mice, while 2/3 remain clinically healthy. To access differential expression of miRNAs, we divided the 4 th generation of our mouse cohort into two separate groups: affected and non-affected mice. Comparing the two groups, only two miRNAs were differentially expressed: miR-379 (adj. P value = 0.044) and miR-223 (adj. P value = 0.044) ( Fig. 6 and Additional file 1: Table S2). Both miRNAs were significantly over expressed in mice with disease. Additionally, an eQTL for miR-379 (−log P = 4.7, 98-116 Mb) was also mapped on chromosome 6, with its peak at~104 Mb (rs6208251, nearest gene: Cntn6). The comparison of diseased and non-diseased mice has its drawbacks; the disease severity observed for individual mice might differ depending on the differences in the genome. To investigate miRNA expression affecting disease severity, we performed co-expression analysis among the expression levels of miRNAs (co-expressed miRNAs are in same pathway) in correlation with quantitative scores for severity and onset week of disease. For this purpose, we employed the WGCNA R package, which clusters miRNAs into different modules and further associates them to a phenotypic score, such as EBA severity or disease onset. As a result, we identified 11 clusters (Fig. 7). Using this approach, only the 'black' module, consisting of 23 miRNAs, was significantly associated with EBA (ρ = 0.28, P = 0.005) (Additional file 1: Table S3). Additionally, it was significantly correlating to the onset of EBA (ρ = 0.29, P = 0.005). Due to the fact that the 'black' module stronger correlates with the onset of the disease than with the maximum score, one can speculate that miRNAs from this module is rather involved in the onset than in the severity of the disease. This would relate to our previous observation, in which miRNA eQTL were overlapping with EBA onset QTL. The pathways associated with the 'black' module were predominantly pathways which have been shown to play a crucial role in other autoimmune disorders, such as the MAPK signaling pathway, T-cell receptor and TGF-beta (Table 3).
To further validate this result we experimentally verified the expression of miR-223 in EBA skin by qRT-PCR. We found that miR-223 is upregulated in EBA skin in comparison to normal skin (Additional file 1: Figure S3). Another miRNA which was highly correlating with EBA was miR-21 (ρ = 0.36, P = 0.00036). Further, investigations of the co-expression module showed that miRNAs were potentially co-regulated i.e. possible overlapping genetic loci among co-expressed miRNAs. We observed that miRNAs in a given locus were either clustered within the same module or showed stronger inter-module membership for a specific module even if they were assigned to different modules. As an example, miR-322 and miR-431, that were mapped on chromosome 1 (51-59 Mb), are clustered in the 'red' module. In eQTL hot spots as on chromosome 2 (28-51 Mb), miR-423-3p and miR-23b are clustered in the 'yellow' module. Even though other miRNAs in this locus were assigned to a different module, they also show significance for the 'yellow' module. Examples are given by miR-671-5p (P = 1.224713e-12), miR-26a (P = 2.654271e-19) and miR-291a-3p (P = 3.530522e-02) (Additional file 1: Table S3). In line with this observation, for the eQTL mapped on chromosome 8, two miRNAs (miR-501-3p and miR-486) were clustered with the 'brown' module while miR-487 was assigns to the 'red' module but had significant module membership with the brown module as well (P = 0.004). Thus, these genomic loci may be considered as confirmed and strengthen the genetic contribution for the control of miRNAs expression levels.

Discussion
Small non-coding RNAs like miRNAs are known to contribute to the onset and severity of various diseases as well as the defense against them [28]. Thus, miRNAs function as tissue-specific key regulators, affecting some of the major pathways towards an aggravation of disease severity when aberrantly expressed [29]. Accordingly, it is not of surprise that miRNAs have been recently recognized as potential therapeutic targets [30,31].
However, the underlying mechanisms of such dysregulated miRNA expression patterns are not well characterized. Different studies have shown that gene expression alterations in different tissues are genetically derived [29]. Thus, it is plausible that not only the regulation of gene expression is genetically controlled, but also the expression of miRNAs. In this study we explore the diversity of miRNAs in inflamed skin tissue and genetic loci that control variations in miRNA expression levels across a mouse cohort. We provide evidence that miRNA levels in skin tissue are genetically controlled on transcriptional level by helicases and RNA polymerases that are important for the biogenesis of miRNA. Furthermore, we found that some of the miRNA eQTL are restricted to one particular locus in the genome (eQTL hot spots). Deeper investigation revealed that these miRNAs are under multi-locus and/or epistatic control.
Interestingly, eQTL hot spots were predominantly found in genomic regions coding for non-coding RNA. Hence, it is tempting to speculate miRNA expression might not necessarily be solely controlled by proteincoding RNA, but rather by non-coding RNA which would impose an additional level of post-transcriptional regulation. This in turn leads to the tempting hypothesis that non-coding RNAs do at least in part regulate miRNA expression. Such a scenario is supported by the fact that some non-coding RNAs have been shown to r Fig. 7 Module-trait relationship between clusters of miRNAs with EBA severity and onset. The graph is a representation of co-expression analyses of miRNAs. The clusters of miRNAs are called modules which are color coded on the y-axis. The disease phenotypes (traits) are EBA onset and EBA severity given on the x-axis. Blocks represent the correlation of module with phenotype using a Pearson correlation coefficient ranging from −1 to 1. The range is color coded with red representing positive correlation and blue representing negative correlation. P values are given in brackets below the correlation coefficients bind to miRNAs at functional level as demonstrated by an interaction of linc-MD1 with miR-133 and miR-135 [32]. Here, linc-MD1 works as a sponge and traps these miRNAs preventing the binding to the canonical targets. Moreover, a recent study even shows an interaction network between lncRNAs and miRNAs [32].
Based on our observed overlap between QTL controlling miRNA expression and EBA, a blistering phenotype in autoimmune skin disease, we conclude that there are interconnected pathways to simultaneously regulate both disease development and miRNA expression. This might explain the findings of earlier studies that show a clear correlation of aberrant miRNA expression and autoimmune diseases [33][34][35][36]. Accordingly, initiation and/or progression of the disease do not primarily appear to be caused by aberrant miRNA expression. Quite the contrary seems to be true; aberrant miRNA expression could be a consequence of the disease which in turn would lead to a downward spiral. Hence, miRNAs could provide a large and unexplored reservoir of potential biomarkers for EBA and related cutaneous autoimmune skin blistering diseases and an interesting target for therapeutic intervention.

Conclusion
Taken together, our study provides a complex framework of gene-gene and miRNA-gene-interactions, which eventually leads to disease development and progression. Our data provide evidence that miRNAs are important drivers of cutaneous autoimmune diseases by acting on various pathways. Moreover, the study strongly implies there is yet another, so far largely unexplored level of regulatory network, possibly comprised of by noncoding RNAs which on their part affect miRNA expression. In this sense, aberrant miRNA expression would indeed be one the responsible elements for disease progression, however, the driving force behind might be a different one.

Generation of a 4-way advanced intercross line
The out bred four-way autoimmune-prone advanced intercross line (AIL) was generated (in our group) from the parental mouse strains BXD2/TyJ, MRL/MpJ, NZM2410/J and CAST/EiJ [21,37,38]. All strains were purchased from The Jackson Laboratory (Bar Harbor, ME). The four inbred strains were intercrossed following

Induction of experimental EBA
Experimental EBA was induced by immunization with an immune-dominant peptide within the murine NC1 domain of Collagen type VII (GST-mCOL7C) as previously described [39]. In total, 600 mice of the AIL were immunized out of which 100 mice were randomly selected for miRNA expression profiling. Mice were evaluated every 4th week after immunization regarding their development and extent of skin disease, following an established scoring system for a total period of 12 weeks [40]. Murine skin tissues were obtained for analysis at the end of the experiment.

Expression profiling and statistical analysis
Total RNA was extracted from skin tissue and processed as previously described using the Flash Tag Biotin HSR RNA labeling kit and hybridized with the miRNA expression profiling GeneChip miRNA 2.0 Array according to the manufacturer's protocol [41]. Raw data was preprocessed using R packages and normalized using RMA (Robust Multi-array analysis) to generate expression levels across the samples following instruction of at bioconductors.org [42].

Genotyping and expression QTL analysis
Genomic DNA was isolated from tail clippings and incubated in 500 μl 50 mM NaOH at 95°C for 2 h. The reaction was neutralized by posterior addition of 50 μl 1 M Tris HCl (pH 8.0). DNA was further processed with DNeasy Blood & Tissue Kit according to manufacturer's instructions. The extracted DNA was quantified using Nanodrop and normalized to 50 ng/μl in TE buffer (10 mM Tris; 1 mM EDTA; pH = 8). Agarose gel electrophoresis was performed for quality control. A total of 1400 SNPs markers evenly spaced across 19 autosomes and the X chromosome were genotyped on 100 generation 4 th mice using Illumina mouse medium density array. 200 SNPs were discarded as they were non informative and 1200 markers were retained. The marker order and position in our map is provided in (Additional file 1). The SNP genotype information for each mouse from generation 4 and four founders is provided in (Additional file 1). These informative SNPs were used for performing eQTL analysis using HAPPY 2.3 on Debian Linux for miRNA expression levels [43,44]. The software infers the haplotype probabilities for each sample using the progenitor SNP information. The miRNA expression levels of mice (600 probes) were considered as quantitative traits and linkage analysis with haplotype probabilities was done using the additive gaussian model. We considered gender as additive covariate. 1000 permutations were performed across all miRNA probes. A significant threshold (α = 0.05) was defined genome-wide. Correction for family structure was done using a variance component model: A Kinship matrix was obtained by the EMMA R package to account for relatedness among the individual [45] and interpreted as a random effect with sex as fixed effect to perform single marker association with phenotype. We defined the confidence interval of the given eQTL based on a -log P value drop of 1.5. We considered local regulation if the eQTL for the probe was found on the same chromosome as its genomic location. A transregulation was presumed if the QTL controlling the trait was found in a different chromosome than its genomic location. The function "epistasis" from the HAPPY R package was used to find pair wise SNP interactions using F-score for the interaction model for SNPs which were significantly associated with miRNA expression levels in initial genomic scan [43]. All the P-values are corrected for multiple testing using the Bonferroni correction. Circos was used for visualization of eQTL and epistasis [46].

Co-expression analysis
The standard WGCNA procedure was used for module detection [47]. For miRNA we considered 97 samples. 3 samples were excluded as these could not phenotype for EBA score. A weighted adjacency matrix of pair-wise connection strengths (correlation coefficients of gene expression levels) was constructed using the soft-threshold approach with a scale independent topological power β = 7 (miRNA). For each probe, the connectivity was defined as the sum of all connection strengths with all others. Probes were aggregated into modules by hierarchical clustering and refined by the dynamic cut tree algorithm [48]. The Pearson correlation coefficient was determined for each phenotype-module pair. The representative module expression profile, or module eigengene value, is the first principal component of the gene expression profile within a module. The correlation between the module eigengene and the sample trait of interest yields the eigengene significance, as assessed by a correlation test. The modules were assigned by different colors where grey was assigned to traits that could not be clustered in any other module. To determine the enriched pathways involved in a cluster we used 'miRsystem' which uses the weighted pathway-ranking method for identifying enriched biological functions [49].