Comprehensive transcriptome analysis of peripheral blood unravels key lncRNAs implicated in ABPA and asthma

Inclusion, sampling and deep RNA sequencing

This study was approved by the Medical Ethics Committee of First Affiliated Hospital of Guangzhou Medical University (ethics approval no. gyfyy-2016-73). All experiments were performed in accordance with relevant guidelines and regulations of the Ethics Committee of First Affiliated Hospital of Guangzhou Medical University. All participants provided written informed consent prior to the publication of clinical and sequencing data. Briefly, a total of 27 unique male individuals were enrolled for our study, including 7 asthma patients (diagnosed as allergic asthma), 12 ABPA patients and 8 healthy individuals considered as control group (Table 1). The diagnosis of asthma was according to the latest Global Initiative for Asthma (GINA) guidelines, and the diagnosis of ABPA was based on the criteria of The International Society of Human and Animal Mycology (ISHAM) working group, which contains two obligatory criteria and three additional criteria (Shah & Panjabi, 2016). Based on the clinical symptoms, 5 asthma patients were at exacerbation stage and two were at chronic stage (Table S1). And according to the proposed clinical staging of ABPA (Shah & Panjabi, 2016), among the included ABPA patients, seven were at acute stage, one was at exacerbation stage, two was at response stage, while the other two were not able to assess the clinical stage because of insufficient clinical data (Table S1).

For each enrolled subject, the peripheral whole blood was extracted and the peripheral blood mononuclear cell (PBMC) was separated immediately by ficoll-paque. The total RNA was extracted using trizol (invitrogen) method and the library for RNA sequencing was prepared based on a standard protocol established by RiboBio Company in Guangzhou and sequenced on Illumina. All the raw RNA-seq data used in the present study are deposited at Short Read Archive (SRA) database of NCBI (https://trace.ncbi.nlm.nih.gov/Traces/sra/) and are assigned the accession number PRJNA582337.

Raw data processing

The raw data was trimmed by Trimmomatic v0.36 (ILLUMINACLIP: TruSeq3-PE.fa:2:30:10:8:true SLIDINGWINDOW:4:15 LEADING:3 TRAILING:3 MINLEN:50) (Bolger, Lohse & Usadel, 2014) and the cleaned reads were aligned against Ensembl hg38 human genome via STAR (v020201) (Dobin et al., 2013). The transcriptome was re-constructed using StringTie (v1.3.3b) (Pertea et al., 2015). A stringent stepwise pipeline (Fig. 1) which had been applied in our previous studies (Huang et al., 2019a; Leng et al., 2019; Zheng et al., 2020) was utilized for identification of high-confidence dataset of lncRNAs.

Identification of lncRNAs

Firstly, the known lncRNAs were picked out from complete set of assembled transcripts according to the “biotype_transcript” of reference gtf file of Homo sapiens from Ensembl database (Zerbino et al., 2018). The transcript was detected as lncRNA when its “biotype_transcript” tagged as “Long_non-coding_RNA”, “Non_coding”, “3prime_overlapping_ncRNA”, “Antisense”, “lincRNA”, “Retained_intron”, “Sense_intronic”, “Sense_overlapping”, “macro_lncRNA” and “bidirectional_promoter_lncRNA”.

Then the human known mRNAs derived from Ensembl database were excluded from the remaining transcripts. Afterwards, the alignment was performed for the remaining transcripts against known protein sequences from NCBI nr database (Pruitt, Tatusova & Maglott, 2007) and Uniprot database (UniProt Consortium, 2018). Successfully aligned transcripts were excluded from potentially non-protein-coding sequences. For the unmapped transcripts, filtering was conducted to get rid of the sequences with length less than 200 nt and the longest ORF longer than 100 residues. Finally, the qualified sequences were further filtered to remove the protein-coding sequences according to the results from Pfamscan (El-Gebali et al., 2019) and CPC (Kong et al., 2007). Eventually, the remaining transcripts were reserved as final dataset of lncRNAs.

Transcriptome variations analysis

For single nucleotide variants, i.e., SNPs and small insertion and deletion (INDEL) detection, we applied Genome Analysis Toolkit (GATK) analysis pipeline (Van der Auwera et al. 2013) for variant calling. All the variants were annotated by ANNOVAR (Wang, Li & Hakonarson, 2010). Simple sequence repeats (SSRs) detection was performed using the Perl script of MISA–MicroSAtellite identification tool (http://pgrc.ipk-gatersleben.de/misa/) with default parameters. Statistical hypotheses were used to investigate whether there were significant differences in the observed proportion of SNPs, INDELs and SSRs among populations of ABPA, asthma and healthy control. The detailed process sees Supplementary Methods.

Transcriptome expression analysis

Differential expression analysis was performed using R package DESeq2 (adjusted p-value ≤ 0.05) (Love, Huber & Anders, 2014). Functional enrichment analysis was performed via ClueGO V2.5.2 (Bindea et al., 2009) plugin in Cytoscape V3.6.1 (Shannon et al., 2003), which could generate a dynamical network structure composed of functionally grouped terms based on the gene list of interest, i.e., differentially expressed genes. The differentially expressed transcripts were subjected to R programming language for hierarchical cluster analysis to detect distinct expression patterns.

Differential interaction network analysis was conducted based on three filtering criteria, which generated three sub-networks for each group, including loss-of-function network, gain-of-function network, and anti-function network. A total of 12 algorithms (MCC, DMNC, MNC, Degree, EPC, BottleNeck, EcCentricity, Closeness, Radiality, Betweenness, Stress, ClusteringCoefficient) of cytoHubba (Chin et al., 2014) plugin in Cytoscape were used to identify hub node for the established network. Detailed procedure see Supplementary Methods. For the hub node validation, the expression profile of four data sets (GSE35571, GSE473, GSE31773, GSE2125) from NCBI GEO database were downloaded and analyzed (Table S2).

Statistics of clinical data

For the basic statistics of clinical data (File S3), all the data were shown in median (interquartile range). Non-parametric Kruskal–Wallis test was performed to compare the differences between three groups while Mann–Whitney test was performed to compare the differences between two groups. Regarding the ratio data, Chi-square test was applied for the comparison. Statistical analysis was performed by SPSS 22.0. A value of P < 0.05 was considered as statistically significant.

Characteristics of study subjects

As shown in Table 1, the asthma group was elder compared with either ABPA group or Healthy control group (P < 0.05). There was no significant difference in eosinophil or neutrophil counts of peripheral blood in the three groups. Regarding the percentages of inflammatory cells in induced sputum, no statistical difference was found in any type of inflammatory cells between asthma and ABPA. There was also no statistical difference in the results of lung function between the two patient groups. In the comparison of IgE levels, patients with ABPA had higher total IgE levels (P < 0.05), higher positive rate of specific IgE (P < 0.05), in particular A. fumigatus specific IgE (P < 0.01) than patients with asthma. In addition, use of glucocorticoids (either oral or inhaled) did not differ significantly between asthma and ABPA.

LncRNAs identification

RNA-seq via Illumina sequencing platform yielded on average 130 million reads per sample (Table S3). After quality trimming, the clean reads were mapped to human genome and subsequently assembled into 283,883 transcripts. A total of 94,155 lncRNAs across 27 samples, including 52,791 known lncRNAs and 41,364 novel lncRNAs were detected using a stringent stepwise filter pipeline (Fig. 1).

No significant difference in variants across ABPA, asthma and health controls

A total of 4,521,799 SNPs and 671,669 INDELs were identified across 27 samples via GATK analysis pipeline (Fig. 1, Table S4). Majority of variants locates in intronic (~61%) and intergenic regions (~26.2%) (Table 2, Fig. S1). A total of 62,733 variants are in genic regions, and 42.9% of them are synonymous SNP (Table S5). 44,090 variants significantly existed differentially across groups by Fisher Exact probability test (File S1). In addition, many genes (including few lncRNAs) exhibited distinct variants cross groups (Fig. 2A).

For SSRs analysis, a total of 4,740,915 different types of SSRs were identified (175,589 per sample on average) (Table S6). No significant difference in SSRs across these samples was detected (Figs. 2D–2I), despites two statistical approaches (details process see Supplementary Methods) revealed around 150 SSRs were differentially enriched in total (Figs. 2B, 2C, File S2).

Differential expression analysis reveals immune-related processes implicated in both diseases

A total of 2,481 differentially expressed transcripts (DETs, including 554 lncRNAs) were found in the comparison of ABPA vs healthy controls (Figs. 3A–3D), and 2,706 transcripts (including 529 lncRNAs) were differentially expressed in the comparison of asthma vs healthy controls, respectively. Functional analysis via ClueGO for both DET datasets demonstrated that many immune-related biological processes were significantly enriched, including leukocyte activation, activation of innate immune response and positive regulation of immune effector process, etc. (Figs. 3B, 3C). Furthermore by the comparison of two DET datasets, distinct DETs involved in diverse biological processes were detected among three different groups (Figs. 3D, 3E).

Additionally, hierarchical cluster analysis via R language based on DETs revealed three main expression patterns of DETs cross these samples (Fig. 4). Concretely, many transcripts (including mRNAs and lncRNAs) are found to be down-regulated in healthy controls but were up-regulated in ABPA and asthma patients (Fig. 4A). Functional analysis showed that they might be involved in myeloid leukocyte mediated immunity, regulation of defense response to virus, etc. (Fig. 4B). On the contrary, several transcripts are down-regulated in healthy controls but are up-regulated in ABPA and asthma patients (Fig. 4E), which might be associated with columnar/cuboidal epithelial cell differentiation, regulation of apoptotic signaling pathway, etc. (Fig. 4F). A few transcripts are specifically up-regulated in ABPA patients (Fig. 4C), and their functions might be correlated to supramolecular fiber organization and positive regulation of organelle organization (Fig. 4D).

Differential network analysis reveals crucial lncRNAs correlated to both diseases

The aforementioned analysis has revealed several lncRNAs and mRNAs differentially expressed in ABPA and asthma patients. Next to investigate the possible roles of these lncRNAs implicated in the pathogenesis of ABPA and asthma, we proposed a differential network analysis. Initially, we built an RNA–RNA interaction network based on the expression profile data of eight healthy controls using the Pearson Correlation Coefficient (PCC) method. Then the interactions of these RNAs in the networks were re-assessed based on the expression profiles of 12 ABPA patients and 7 asthmatic patients. To assess the changes of network from healthy status to disease status, three different filtering criteria (details see Supplementary Methods) was proposed to construct three sub-networks for each disease, namely loss-of-function network, gain-of-function network, and anti-function network.

In loss-of-function network, for instance, 293,926 RNA–RNA interactions with strong correlation (|PCC value| ≥ 0.70 & P-value < 0.05) are found in health controls, whereas these interactions are shown having weak correlation (|PCC value| ≤ 0.30) in ABPA patients. It is suggested these RNA–RNA interactions might be disrupted by the occurrence of ABPA. 243,399 RNA–RNA interactions exhibit similar situation in asthma. For each disease, the immune-related RNA–RNA interactions in two loss-of-function networks were extracted to establish the loss-of-immune-related-function network, respectively (Figs. 5A, 5B). Key elements (hub nodes of the network) in both networks were detected by topological network analysis via Cytohubba. Our analysis revealed several hub lncRNAs which might play key roles in the immune dysfunction of ABPA and asthma. In this way, we built the gain-of-function and anti-function networks for both diseases (Figs. 5C–5F), and relevant key lncRNAs and mRNAs were disclosed (Table S7, S8).

To further investigate convincing RNA–RNA interactions implicated in ABPA and asthma, we integrated the above-mention networks (Fig. S2). Notably, the RNAs in the network can be categorized into four main immune-related groups, including leukocyte involved immune process, cytokine related immune process, autophagy apoptosis regulation and response to exogenous invasion. Topological analysis identified many key mRNAs/lncRNAs implicated in the regulation of diverse immune-related processes in response to ABPA and asthma (the frequency of counting as top 20 hub nodes among 12 algorithms ≥ 3, Table 3, details see Supplementary Methods), including five lncRNAs, namely AL139423.1-201, AC106028.4-201, HNRNPUL1-210, PUF60-218 and SREBF1-208, and five known mRNAs (PRR14L-201, CCRL2-202, GDAP2-202, SPAG1-206, SLC35B2-207). Notably, many mRNAs have already been demonstrated to be associated to asthma or ABPA by previous studies. For instance, Inhibition of PARP14 (Top one hub node in Table 3, it was found to significantly up-regulated in both diseases) was reported to reduce allergic airway diseases, it thereby was proposed as potential therapeutic target for asthma (Mehrotra et al., 2013).

Moreover, all the RNA–RNA interactions involved the five key lncRNAs were observed in the loss-of-function networks in both diseases and exhibited down-regulated expression in the peripheral blood of patients in both diseases (Table 3), suggesting the regulations of hub lncRNAs in these interacted target RNAs were disrupted after the occurrence of asthma or ABPA. Their targets were found to be involved in many immune-related processes and pathways (Figs. 6A–6E). Permutation test shows majority of hub nodes are high-confidence (for details see Supplementary Methods). These findings indicate that the expression of these five lncRNAs would be repressed during the occurrences of ABPA or asthma, many immune-related processes mediated by them might be disrupted, simultaneously.

Validation of hub nodes using GEO datasets

Four hub nodes including two mRNAs (ENST00000327423, ENST00000369443) and two lncRNAs (ENST00000587128, ENST00000595806) from the two networks were selected for validation. The expression profiles of four GEO datasets (GSE35571, GSE473, GSE31773 and GSE2125) were used. The results indicated that the expression level of the selected hub nodes are consistent with the result of GEO datasets (Figs. 6F–6I). Also, many interactions involved in these hub nodes could be found in these GEO datasets (Table S9), implying our predictions are reliable and convincing.