RNA sequencing reveals the emerging role of bronchoalveolar lavage fluid exosome lncRNAs in acute lung injury

Animals

The experiments were performed adhering to the institutional guidelines and approved protocols. The animal experiments were approved by the Institutional Animal Care and Use Committee of Nanjing Medical University (No. IACUC-2004021). All animal experiments were conducted at the Animal Core Facility of Nanjing Medical University. Twelve male Wistar rats (weighing 200 ± 20 g and aged 6–8 weeks) were purchased from Nanjing Qing Long Shan animal farm (Nanjing, China).

During the whole experiment, the rats in the control and the experimental groups had free access to food and water. The health and immune statuses of all rats used were normal, and they were not involved in any previous procedures. The rats were randomly grouped (random number table method) as ALI models and controls (n = 6 per group). ALI was modelled and the sample size was determined as previously described (Do-Umehara et al., 2013; Lu et al., 2012). Lipopolysaccharide (LPS) was dissolved in 0.5 mL of normal saline to obtain a solution at 10 mg/kg of body weight. After anesthesia with 3% sodium pentobarbital (50 mg/kg), the rats were placed in a supine position on the operating table and airway-injected with the LPS solution. An equal volume of normal saline was given to the rats in the control group. All Wistar rats were placed under the same conditions for 24 h and given the same anesthesia. All animals were given humane care.

For confirming the properties of the exosomes, the exosomes purified from the ALI group were resuspended in 200 µL of phosphate-buffered saline (PBS) and infused into the lungs of two rats. The rats in the control group was infused with PBS alone. The histological examination was performed 24 h later. Only investigator who performed modeling was aware of grouping but was not involved in the subsequent experiments or analyses.

BALF sampling and histopathological analysis

Twelve Wistar rats were divided into the experimental (n = 6) and control (n = 6) groups. Anesthesia was performed with 3% pentobarbital sodium (50 mg/kg). After successful confirmation of endotracheal intubation using alc-8 small-animal ventilator, 5 mL of normal saline (0.9%) was injected into the airway. Through airway intubation, the right lung was ligated and the left lung was irrigated with 4 °C pre-cooled saline. This was repeated four times, and the BALF was collected in centrifuge tubes. Once the BALF was obtained, the rats were sacrificed by cervical dislocation and the lungs were harvested. The left lung was weighed (wet weight), placed in an oven at 65 °C for 7 d, and then weighed again to determine the dry weight. The dry-to-wet ratio was calculated. The right lung was formalin-fixed and paraffin-embedded. The sections (4 µm) were cut and stained with hematoxylin and eosin.

Enzyme-linked immunosorbent assay

Commercial enzyme-linked immunosorbent assay were used to measure the levels of interleukin (IL)-6 and tumor necrosis factor (TNF)-α from 12 rats (n = 6/group) (R&D Systems, Minneapolis, MN, USA) following the manufacturer’s protocol.

Extraction of exosomes from BALF

The BALF exosomes were purified following the ISEV guidelines (Deady et al., 2014; Théry et al., 2018). This includes determining the speed of ultracentrifugation based on rotor type, tube/adapter, and centrifuge speed. Second, the pore size of the matrix should be considered. For example, a group of vesicles may be excluded if the pore size does not include EVs >70 nm in diameter. As well as EV characterization based on protein content, at least one of 1a (CD63, CD81, CD82, etc.) or 1b (ERBB2, EPCAM, CD90, etc.), 2a (TSG101, HSPA8), 3a (APOA1/2, APOB; APOB100, etc.) or 3b (Tamm-Horsfall protein, UMOD) class proteins must be analysed to demonstrate the properties of EVs and the purity of EV preparations. For this, 15 mL of BALF samples were centrifuged for 10 min at 2,000g. The supernatant was centrifuged for 20 min at 12,000g (Optima L-100XP Ultracentrifuge, Beckman Coulter, Brea, CA, USA). The supernatant was centrifuged again for 70 min. After centrifugation, the supernatant was discarded, and the precipitate (exosomes) was resuspended in 200 µL of PBS in a 1.5-ml eppendorf tube and stored at –80 °C.

Exosome properties

A Tecnai G2 Spirit BioTwin Nano Transmission Electron Microscope (FEI, Hillsboro, OR, USA) detector was used to examine the exosome morphology. A nanoparticle size detector was used to detect exosome particle size.

Western blotting

Exosome surface proteins (CD63 and TSG101) were examined by Western blot (Keerthikumar et al., 2016; Lee, El Andaloussi & Wood, 2012; Logozzi et al., 2009). Equal amounts of proteins from the samples were resolved by sodium dodecyl sulfate–polyacrylamide gel electrophoresis and transferred on to a poly-polarized PVDF membrane. The blot was incubated with the primary antibodies overnight: mouse anti-GAPDH (1:5,000; ab8245; Abcam, Cambridge, United Kingdom), mouse anti-CD63 (1:1,000; sc-5275; Santa Cruz Biotechnology, Santa Cruz, CA, USA), and rabbit anti-Tsg101 (1:5,000; ab125011; Abcam, Cambridge, United Kingdom). The secondary antibody was HRP-conjugated goat anti-rabbit IgG or goat anti-mouse IgG (both 1:50,000; Wuhan Boster Biological Technology, Ltd., Wuhan, China). The bands were revealed using an ECL reagent (Pierce Chemical, Dallas, TX, USA). The film was scanned and analysed using BandScan.

RNA-seq

Four samples were randomly selected from the two sets of samples for high-throughput transcriptome sequencing. We carried out quality inspection on the sample RNA, and explained the detection index RIN (RNA Integrity Number) of RNA integrity. RIN ranges from 0 to 10. The higher the score, the better the integrity of the RNA. The RNA quality inspection result of our sample is >=7.0, which is a qualified sequencing sample, and the base distribution was balanced. For raw reads that might contain unqualified reads with low overall quality, sequencing primers, low end quality, and so forth, we applied Seqtk (https://github.com/lh3/seqtk) to filter them to obtain clean reads that could be used for data analysis. The RNeasy mini kit (Qiagen, Venlo, The Netherlands) was used to isolate total RNA. The TruSeq RNA Sample Preparation Kit (Illumina, Inc., San Diego, CA, USA) was used to synthesize paired-end libraries. The poly-A-containing mRNA molecules were purified using poly-T oligo-attached magnetic beads. A Qubit 2.0 Fluorometer (Life Technologies Co., Grand Island, NY, USA) was used to quantify the purified libraries, which were validated using an Agilent 2100 bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The cluster was generated using cBot with the library diluted to 10 pM. The cluster was sequenced on an Illumina HiSeq X-ten (Illumina, Inc., San Diego, CA, USA). Shanghai Biotechnology Corporation (Shanghai, China) performed library construction and sequencing.

Unqualified reads were filtered to obtain clean reads for data analysis using Seqtk (https://github.com/lh3/seqtk) for filtering (version 2.2.8). The reads were preprocessed by filtering out rRNA reads, sequencing adapters, short-fragment reads, and other low-quality reads using Hisat2 (version 2.0.4) (Kim, Langmead & Salzberg, 2015) to map the cleaned reads to the human GRCh38 reference genome with two mismatches. The novel lncRNA and NONCODE database (version: NONCODE 2016; http://www.noncode.org/) were predicted using Stringtie (version:1.3.0) (Pertea et al., 2015, 2016), and the known data in the Ensembl database lncRNA were used for expression quantification. The ID starting with MSTRG was novel lncRNA, the ID starting with NON was the known lncRNA in the database, and the ID starting with ENS was the known lncRNA in the Ensembl database.

Stringtie (version 1.3.0) was run with a reference annotation to generate fragments per kilobase of exon model per million mapped reads (FPKM) values for known gene models. Differentially expressed genes were identified using edgeR (Robinson, McCarthy & Smyth, 2010). The P value was set using the false discovery rate (FDR) (Benjamini et al., 2001; Benjamini & Hochberg, 1995; Benjamini & Yekutieli, 2001). The fold-changes were also estimated according to the FPKM in each sample. The differentially expressed genes were selected using the following filtering criteria: FDR ≤0.05 and fold-change (FC) 195 ≥ 2.

StringTie (Pertea et al., 2015, 2016) (version: 1.3.0) was applied to quantify the expression of novel lncRNAs and NONCODE databases predicted using 2.2.11 (version: NONCODE 2016; http://www.noncode.org/), as well as examine known lncRNAs in Ensemble database.

GO and KEGG analysis of differentially expressed lncRNAs

The reads were converted into FPKM for standardized gene expression levels (Mortazavi et al., 2008; Robinson, McCarthy & Smyth, 2010) for comparisons between groups. The differentially expressed lncRNAs were used for Gene Ontology (GO) enrichment analysis (http://geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment (http://www.genome.jp/kegg). All experiments were performed three times independently. Subsequently, five differentially expressed lncRNAs were randomly selected for validation.

Real-time quantitative reverse transcription–polymerase chain reaction

The miRNeasy Micro Kit kit (Qiagen, Venlo, The Netherlands) was used to extract exosomal total RNA. Table 1 presents the primers for quantitative polymerase chain reaction (qPCR). The amplification parameters were 95 °C for 10 s and 60 °C for 34 s, for a total of 40 cycles. The relative expression levels of lncRNAs were calculated using the 2^–ΔΔCt method (Livak & Schmittgen, 2001).

Prediction of the lncRNA–miRNA interaction networks

Five lncRNAs were selected to construct lncRNA–miRNA networks using the miRDB Database to investigate the regulation network between lncRNAs and their target miRNAs. The Cytoscape software (version 3.7.1, https://cytoscape.org/) was used for network visual representation.

Statistical analysis

Data were tested for normal distribution using the Shapiro–Wilk test. They were presented as means ± standard deviations and analyzed using Student t test. All analyses were performed using SPSS 17.0 (IBM, New York, NY, USA). P values < 0.05 indicated statistically significant differences. The rats that did not meet the ALI standard were excluded.

ALI modeling and the pro-inflammatory effects of exosomes

Compared with the control group (n = 6), the lung tissues in the ALI group (n = 6) showed a significantly smaller alveolar cavity, more extensive alveolar space, and infiltration of many neutrophils in the alveolar wall (Fig. 1A). Compared with the control group, the ALI group showed a higher wet/dry ratio, TNF-α levels, and IL-6 levels (all P < 0.001). These results indicated that exosomes were involved in mediating inflammatory responses in ALI (Fig. 1B). In addition, our study found that the airway injection of exosomes in rats led to significant infiltration by neutrophils, with smaller alveolar cavities and full alveolar septum (Fig. 2A). Also, the wet and dry lung weight of rats was significantly higher in the exosome group than in the PBS group (Fig. 2B).

Exosome confirmation

Nano transmission electron microscopy showed that the diameter of the exosomes, shown as clear vesicle-like structures, was mainly between 40 and 200 nm, primarily around 100 nm; also, they were larger in the ALI group (n = 2) (Figs. 3A and 3B). Both exosome surface proteins (CD63 and Tsg101) were shown as positive by Western blot (Fig. 3B). All these findings confirmed the successful extraction of exosomes from BALF.

High-throughput sequencing results and analysis

The RNA of the exosomes extracted from the BALF in the ALI and control groups was sequenced using high-throughput sequencing (uploaded to NCBI, #SUB7338616). A total of 2,958 differentially expressed lncRNAs were identified, including 2,524 upregulated and 434 downregulated ones. The results were summarized as scatter diagram (Fig. 4A), volcano plot (Fig. 4B), and heatmap (Fig. 4C).

GO and KEGG database analyses

We conducted an in-depth analysis of the sequencing results. Our analysis showed that there were 5,500 differentially expressed mRNAs between the two groups, of which 2,717 were differentially up-regulated and 2,783 were down-regulated. The mRNAs directly bound to lncRNA and the differentially expressed mRNAs downstream of the differentially expressed lncRNAs (including some novel lncRNAs) were analyzed, and GO and KEGG pathway enrichment analyses were performed on the results. GO enrichment analysis was performed on 2,958 differentially expressed lncRNAs identified. The gene number distribution of top 30 genes in GO analysis is shown in Fig. 5A. As can be seen from the scatter diagram, the three functions with the most significant number of genes included phagocytic vesicle membrane, regulation of receptor biosynthesis process, and I-SMAD binding. Using the same screening criteria as GO analysis, differentially expressed genes for signaling pathways were analyzed using the KEGG database analysis. Salmonella infection, Toll-like receptor signaling pathway, and osteoclast differentiation were the most enriched pathways (Fig. 5B). In addition, we list the top 30 differentially up-regulated and differentially down-regulated lncRNAs (Table 2).

qRT-PCR verification

Of the first 30 differentially expressed lncRNAs, four lncRNAs (NONRATT002967.2, NONRATT003362.2, NONRATT004060.2, NONRATT025040.2, and NONRATT025699.2) starting with NONRATT (i.e., novel lncRNAs) were randomly selected for RT-PCR validation. The PCR validation and sequencing results were consistent (n = 6) (Fig. 6).

We analyzed the target genes of NONRATT002967.2, NONRATT004060.2, and NONRATT025699.2 as Tpbg1, Tceb2, and Igf1, respectively. In addition, we also used RT-PCR to verify NONRATT025040.2, whose mechanism will be explored in the future. Its target gene was Foxa1, which was mainly responsible for regulating the differentiation of lung epithelial cells. The results showed that NONRATT025040.2 decreased in the ALI group compared with the control group.

LncRNA–miRNA interaction networks

Given that lncRNAs can bind to miRNAs and work as a miRNA “sponge”, the relation of the five lncRNAs and possible binding miRNAs was investigated. The lncRNA–miRNA interaction network of the five lncRNAs was predicted using miRDB (http://mirdb.org/) and visualized using Cytoscape (Fig. 7). NONRATT004060.2, NONRATT002967.2, NONRATT025699.2, NONRATT025040.2, NONRATT003362.2 interact with 10 miRNAs respectively. NONRATT025699.2 is closely related to NONRATT025040.2, and there are eight miRNAs that work together, including hsa-miR-1237-5p, hsa-miR-4488, hsa-miR-365a-3p, hsa-miR-25-5p, hsa-miR-365b-3p, hsa-miR-6746-5p, hsa-miR-6763-3p, hsa-miR-4697-5p. However, NONRATT003362.2 has no interaction relationship with the other four lncRNAs.