Integrative genetic and genomic networks identify microRNA associated with COPD and ILD

Abstract

Chronic obstructive pulmonary disease (COPD) and interstitial lung disease (ILD) are clinically and molecularly heterogeneous diseases. We utilized clustering and integrative network analyses to elucidate roles for microRNAs (miRNAs) and miRNA isoforms (isomiRs) in COPD and ILD pathogenesis. Short RNA sequencing was performed on 351 lung tissue samples of COPD (n = 145), ILD (n = 144) and controls (n = 64). Five distinct subclusters of samples were identified including 1 COPD-predominant cluster and 2 ILD-predominant clusters which associated with different clinical measurements of disease severity. Utilizing 262 samples with gene expression and SNP microarrays, we built disease-specific genetic and expression networks to predict key miRNA regulators of gene expression. Members of miR-449/34 family, known to promote airway differentiation by repressing the Notch pathway, were among the top connected miRNAs in both COPD and ILD networks. Genes associated with miR-449/34 members in the disease networks were enriched among genes that increase in expression with airway differentiation at an air-liquid interface. A highly expressed isomiR containing a novel seed sequence was identified at the miR-34c-5p locus. 47% of the anticorrelated predicted targets for this isomiR were distinct from the canonical seed sequence for miR-34c-5p. Overexpression of the canonical miR-34c-5p and the miR-34c-5p isomiR with an alternative seed sequence down-regulated NOTCH1 and NOTCH4. However, only overexpression of the isomiR down-regulated genes involved in Ras signaling such as CRKL and GRB2. Overall, these findings elucidate molecular heterogeneity inherent across COPD and ILD patients and further suggest roles for miR-34c in regulating disease-associated gene-expression.

Figure 1

Heterogeneity of miRNA expression profiles associated with COPD or ILD. (A) The expression profiles of 255 miRNAs were significantly associated with the presence of disease (ANOVA FDR q-value < 0.10 and fold change > 1.25 in either disease group compared to controls). Consensus clustering was used to identify 5 distinct samples clusters and 4 distinct miRNA clusters. (() Stacked barplots display the proportion of disease and control samples within each sample cluster. The majority of control samples (53%) fell into cluster S1. Clusters S2, S4, and S5 were enriched with ILD patients compared to cluster S1. (C) Meta-scores were derived for each miRNA module with PCA and used to demonstrate the relationship between miRNA module expression and sample clusters. Differences in the module meta-scores between sample clusters were determined using a linear model controlling for reads aligned, protocol, smoking status, gender, age, and overall disease status. Sample cluster S1 was treated as the baseline group since it had the highest proportion of Control samples.

Figure 2

Association of sample clusters with clinical features. Differences in clinical variables between sample clusters were determined using a linear model controlling for reads aligned, protocol, smoking status, gender, and age. Sample cluster S1 was treated as the baseline group since it had the highest proportion of Control samples. (A) When including all samples, sample clusters S2-S5 has significantly lower DLCO compared to S1. (B–D) Sample cluster S3 had significantly lower FEV1 Percent Predicted as well as significantly higher Percent Emphysema and BODE compared to S1. (E–J) When examining differences within a disease, ILD samples in clusters S4 and S5 has significantly lower DLCO compared to ILD samples in cluster S1. COPD samples in cluster S3 had significantly lower DLCO, FEV1 percent predicted, FEV1/FVC ratios and significantly higher percent emphysema and BODE scores compared to COPD samples in cluster S1. Asterisks indicate significance: *p < 0.05; **p < 0.01; ***p < 0.001.

Figure 3

Examining top miRNA in disease-specific integrative networks. (A) We select those SNP-miRNA-mRNA triplets where the SNP-mRNA relationship is defined by a miRNA mediator; we filter out independent relationships and those triplets where the SNP is not associated with the miRNA. (B) The CIT networks follow a power law. The negative correlation between the frequency of node degree and the node degree indicates that the networks are scale-free. (C) Number of genes regulated by each miRNA. miR-449/34 family members were found to be among the top 20 differentially connected in COPD and ILD compared to control group. The red dots indicate the significantly differentially connected miRNAs by a Fisher’s exact test (FDR < 0.2).

Figure 4

Enrichment of miR-449/34 modules. (A) Clustering of miRNA modules based on the Jaccard index revealed a group of strongly overlapping miRNA in the miR-449/34 family. (B) GSVA was used to predict the activity of the miR-449/34 family in a gene-expression dataset of airway epithelial differentiation. The set of genes that positively correlated with miR-449/34 family (406 genes) were enriched among genes that increase in expression with the airway epithelial cells differentiation (p < 0.05; Linear mixed-effects model). (C) Similarly, enrichment of miR-449/34 gene set family with the airway cells differentiation is shown by GSEA (FDR q-value < 0.001).

Figure 5

Functional roles of an isomiR for hsa-miR-34c. (A) The top 25 highest expressed sequences are shown for the hsa-miR-34c locus. Three of these sequences represented an isomiR on the 5′ end which contained a non-canonical miRNA seed (purple). (B) Predicted mRNA targets for the miR-34c-5p canonical seed (red), the 5′ isomiR (blue), or both seeds (green) were more negatively correlated than genes that were not predicted targets of either seed (black). P-values were determined using Wilcoxon rank-sum tests that compared the correlation coefficients for each group of predicted miR-34c-5p targets (blue, red, or green) to the correlation coefficients for genes without any target seed (black). (C) The overlap of negatively correlated (FDR < 0.25) and 5′ UTR-site predicted targets of miR-34c and 5′ isomiR. miR-34c targets are significantly enriched for Notch signaling pathway by Enrichr (p < 0.02); 5′ isomiR targets are significantly enriched for Ras signaling pathway by Enrichr (p < 0.002). (D) IMR90 fibroblast cells show significant repression of NOTCH1 with all experimental transfections (p < 0.005, p < 0.005, p < 0.005). (E) NOTCH4 expression is only downregulated by miR-34c mimic transfection (p < 0.005) and not by transfection of the isomiR. (F) GRB2 is significantly downregulated with the miR-34c 5′ isomiR mimic transfection (p < 0.005, p < 0.005) but not the miR-34c canonical transfection. (G) CRKL is also significantly downregulated with the miR-34c 5′ isomiR mimic transfection (p < 0.05, p < 0.05) but not the miR-34c canonical transfection.