Gene expression signatures of target tissues in type 1 diabetes, lupus erythematosus, multiple sclerosis, and rheumatoid arthritis

Abstract

Autoimmune diseases are typically studied with a focus on the immune system, and less attention is paid to responses of target tissues exposed to the immune assault. We presently evaluated, based on available RNA sequencing data, whether inflammation induces similar molecular signatures at the target tissues in type 1 diabetes, systemic lupus erythematosus, multiple sclerosis, and rheumatoid arthritis. We identified confluent signatures, many related to interferon signaling, indicating pathways that may be targeted for therapy, and observed a high (>80%) expression of candidate genes for the different diseases at the target tissue level. These observations suggest that future research on autoimmune diseases should focus on both the immune system and the target tissues, and on their dialog. Discovering similar disease-specific signatures may allow the identification of key pathways that could be targeted for therapy, including the repurposing of drugs already in clinical use for other diseases.

Fig. 1. Overview of the number of differentially expressed genes and the signaling pathways activated in the target tissues of four autoimmune diseases.

(A) Number of protein-coding genes differentially expressed in four autoimmune diseases. Each RNA-seq data set was quantified with Salmon using GENCODE 31 as the reference. Differential expression was assessed with DESeq2. The numbers within the bars represent the protein-coding genes with |fold change| >1.5 and an adjusted P value <0.05. RNA-seq sample numbers (n) are as follows: T1D (n = 4 for patients, n = 10 for controls), SLE (n = 20 for patients, n = 7 for controls), MS (n = 5 for patients, n = 5 for controls), and RA (n = 56 for patients, n = 28 for controls). Results for the RA samples were adjusted by age as an independent variable. (B to E) Gene set enrichment analysis (GSEA) of T1D (B), SLE (C), MS (D), and RA (E) target tissues. After quantification using Salmon and differential expression with DESeq2, genes were ranked according to their fold change. Then, the fGSEA algorithm (76) was used along with the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome databases to determine significantly modified pathways. Bars in red and blue represent, respectively, a positive and negative enrichment in the associated pathway. The x axis shows the normalized enrichment score (NES) of the fGSEA analysis, and the y axis the enriched pathways. The numbers at the end of the signaling pathway names represent, respectively, (i) the number of genes present in the leading edge of the GSEA analysis and (ii) the total number of genes present in the gene set.

Fig. 2. RRHO analysis comparing the gene expression signatures of target tissues among four autoimmune diseases demonstrates a high degree of similarity between up-regulated genes.

(A) Genes were ranked by their fold change from the most down- to the most up-regulated ones and then submitted to the RRHO algorithm. The level map colors display the adjusted log P values of the overlap (the P values were adjusted using the Benjamini and Yekutieli method) between genes up-regulated in both datasets (bottom left quadrant), down-regulated in both (top right quadrant), up-regulated in the left-hand pathology and down in the bottom part (top left quadrant), and down in the left-hand pathology and up-regulated in the bottom part pathology (bottom right quadrant). (B) The panel displays the number of genes significantly overlapping in each pairwise analysis (A). NS, not significant quadrant.

Fig. 3. Functional enrichment analysis of overlapping genes among four autoimmune diseases demonstrates signaling pathways concordance.

(A to F) Genes significantly overlapping between different pairs of autoimmune diseases in the RRHO analysis (Fig. 2B) were selected for enrichment analysis using the clusterProfiler tool with the Reactome database. The top 20 gene sets are represented according to their adjusted P values (Benjamini and Hochberg correction) and their gene ratio (no. of modified genes/total gene set size). Diseases were analyzed in pairs. Enrichment analysis of genes significantly up-regulated in the target tissues of both (A) T1D and SLE, (B) T1D and MS, (C) T1D and RA, (D) SLE and MS, (E) SLE and RA, and (F) MS and RA.

Fig. 4. Mining overlapping genes among target tissues in four autoimmune diseases allows the identification of potential common therapeutic targets.

(A to F) After determining statistically which genes were overlapped in pairs of autoimmune diseases from the RRHO analysis (Fig. 2), the top 150 overlapping genes were submitted to the Connectivity Map database to identify perturbagen classes driving an opposite signature (negative tau score) to the one present in the target tissues of the four autoimmune diseases. Only classes with a median tau score <−80 were considered. (A to F) Perturbagen classes driving down the genomic signatures of up-regulated genes. The same methodology and conditions have been applied for subsequent analysis: (A) T1D and SLE, (B) T1D and MS, (C) T1D and RA, (D) SLE and MS, (E) SLE and RA, and (F) MS and RA. EGFR, epidermal growth factor receptor. LOF, Loss of Function; GOS, Gain of Function; IAP, inhibitor of Apoptosis; FGFR, Fibroblast Growth Factor Receptor; MDM, Murine Double Minute; HIF, Hypoxia Inducible Factor; BCL, B-Cell Lymphoma.

Fig. 5. Venn diagrams of risk genes expressed in the target tissues of the four autoimmune diseases shows shared candidate genes among them.

Venn diagram representing risk genes identified in GWAS studies in target tissues for T1D, SLE, MS, and RA. For each disease, the risk genes were extracted from the GWAS Catalog (www.ebi.ac.uk/gwas/) and selected as described in Materials and Methods. In brief, each list was curated according to their relationship to the disease, and only genes with a P value <0.5 × 10⁻⁸ for their SNP-trait relationship were kept. Last, an intersection between the four lists was performed and represented as a Venn diagram. Numbers in the diagram represent the numbers of genes present in each subgroup, and genes overlapping between diseases are displayed by their HGNC symbols. A gene was considered as expressed if it presents a mean TPM > 0.5 in either the patient or control group. N/A, not applicable (no gene in common).