A network-based approach reveals long non-coding RNAs associated with disease activity in lupus nephritis: key pathways for flare and potential biomarkers to be used as liquid biopsies

Abstract

Objective: A blood-based biomarker is needed to assess lupus nephritis (LN) disease activity, minimizing the need for invasive kidney biopsies. Long non-coding RNAs (lncRNAs) are known to regulate gene expression, appear to be stable in human plasma, and can serve as non-invasive biomarkers.

Methods: Transcriptomic data of whole blood samples from 74 LN patients and 20 healthy subjects (HC) were analyzed to identify differentially expressed (DE) lncRNAs associated with quiescent disease and flares. Weighted gene co-expression network analysis (WGCNA) was performed to uncover lncRNAs with a central role (hub lncRNAs) in regulating key biological processes that drive LN disease activity. The association of hub lncRNAs with disease activity was validated using RT-qPCR on an independent cohort of 15 LN patients and 9 HC. cis- and trans-targets of validated lncRNAs were explored in silico to examine potential mechanisms of their action.

Results: There were 444 DE lncRNAs associated with quiescent disease and 6 DE lncRNAs associated with flares (FDR <0.05). WGCNA highlighted IFN signaling and B-cell activity/adaptive immunity as the most significant processes contributing to nephritis activity. Four disease-activity-associated lncRNAs, namely, NRIR, KLHDC7B-DT, MIR600HG, and FAM30A, were detected as hub genes and validated in an independent cohort. NRIR and KLHDC7B-DT emerged as potential key regulators of IFN-mediated processes. Network analysis suggests that FAM30A and MIR600HG are likely to play a central role in the regulation of B-cells in LN through cis-regulation effects and a competing endogenous RNA mechanism affecting immunoglobulin gene expression and the IFN-λ pathway.

Conclusions: The expression of lncRNAs NRIR, KLHDC7B-DT, FAM30A, and MIR600HG were associated with disease activity and could be further explored as blood-based biomarkers and potential liquid biopsy on LN.

Keywords: RNA-sequencing; SLEDAI-2K; WGCNA; blood-based biomarker; ceRNA; disease activity; long non-coding RNAs; lupus nephritis.

Figure 1

(A) Graphical overview of our research steps. (B) Volcano plot of the three differential expression analyses (DEA) performed comparing (left to right) Lupus Nephritis patients vs. Healthy controls, High Disease Activity vs. Inactive Disease patients, and Inactive Disease patients vs. Healthy controls. Upregulated genes are colored violet, and downregulated genes are colored blue. Genes not reaching our significance thresholds (|log2FC| >0.58 and FDR <0.05) are shown in gray. (C) Bubble plot showing inflammation-related Gene Ontology terms found as significantly enriched in each of the three DEA when performing gene set enrichment analysis (GSEA). Color represents FDR values whereas bubble size represents the Normalized Enrichment Score of each term. (D) Venn diagram comparing the DE genes of each DEA. Color gradient corresponds to the gene count in each compartment. (E) Venn diagram comparing the enriched terms of each GSEA. Color gradient corresponds to the term count in each compartment.

Figure 2

(A) Volcano plot of the long non-coding RNAs (lncRNAs) in each of the three differential expression analyses (DEA) performed comparing (left to right) LN patients vs. HC, HDA vs. InaD patients, and InaD patients vs. HC. Upregulated lncRNAs are colored violet, and downregulated lncRNAs are colored blue. LncRNAs not reaching our significance thresholds (|log2FC| >0.58 and FDR <0.05) are shown in gray. (B) Heatmap showing the expression profile of the top 250 lncRNAs with the highest absolute log 2 fold change value found as DE between LN and HC. Expression values were z-score normalized. Top annotation row shows the condition of each sample, colored green for LN patients and light blue for HC. (C) Heatmap showing the expression profile of the top 250 lncRNAs with the highest absolute log 2 fold change value found as DE between InaD patients and HC. Expression values were z-score normalized. Color scale of top annotation is the same as (B). (D) Heatmap showing the expression profile of the six lncRNAs found as DE between HDA and InaD patients. Expression values were z-score normalized. Top annotation row shows the disease activity group of each sample with black representing HDA patients and gray representing InaD patients. (E) Venn diagram comparing the DE lncRNAs of each DEA. Color gradient corresponds to the term count in each compartment.

Figure 3

(A) Heatmap showing the correlation of the eigengene of each module found to be significantly correlated with SLEDAI-2K. Color corresponds to correlation level with purple for positive and blue for negative correlation. (B) Bubble plot of Gene Ontology terms found as significantly enriched in the top three correlated modules (Lightgreen, Salmon, Lightyellow). Color represents adjusted p-values, and size represents the number of genes related to a term found in each module. (C) Scatterplot of SLEDAI-2K Gene Significance against Module Membership for each gene in the (left to right) Salmon, Lightgreen, and Lightyellow modules. Genes with MM >0.8 (Hub genes) are shown in color, with violet for lncRNAs and blue for other RNA types. Genes with MM ≤0.8 (non-hub genes) are shown in gray. (D) Heatmap showing the correlation of the RNA-Seq-based expression values of the nine hub lncRNAs that were significant when tested using the Spearman correlation coefficient. Color corresponds to correlation level with purple for positive and blue for negative correlation. (E) Scatterplots of expression levels of FAM30A (top left), KLHDC7B-DT (top right), MIR600HG (lower left), and NRIR (lower right) normalized using z-score scaling per experiment type (qPCR, RNA-Seq) against SLEDAI-2K values. Boxes on top of each plot show the Spearman correlation coefficient and the associated p-value. Colors correspond to experiment type with blue for qPCR, gold for RNA-Seq.

Figure 4

(A) Plot of the genomic region surrounding FAM30A. The genomic region depicted corresponds to 12 kbp upstream of the FAM30A start position and 12 kbp downstream its end position. Identified transcripts of genes found in the area are shown in gold (exons) connected by gray lines with arrows (introns). The exact position of the locus in the human genome is marked by the red line on the right side of the Chromosome 14 ideogram on the top of the figure. (B) Bar plot showing the percentage of immunoglobulin (IG) genes found in each WGCNA module. Each bar is colored according to the module name. (C) Network of ceRNA interactions of FAM30A and MIR600HG. Node fill color corresponds to RNA type with gold for LncRNA, red for miRNA, and light blue for mRNA. The node outline is colored depending on whether the node is connected to both FAM30A and MIR600HG (common—purple) or just one of the two lncRNAs (unique—black). Node size is a function of each degree with highly connected nodes shown as bigger points. The network layout was created using the Davidson and Harels simulated annealing algorithm of the igraph package.