An integrated analysis of lncRNA and mRNA expression profiles in the kidneys of mice with lupus nephritis

Animals

Mice of the strain MRL/MpJ-Fas-lpr (Mrl/lpr), an established model of LN, were used in our study. Control MRL/MpJ mice without the Fas-lpr mutation were used as healthy controls. Specific pathogen-free (SPF)-grade female MRL/lpr and MRL/MpJ mice (weighing 15–16 g and aged 4 weeks) were purchased from SLAC Laboratory Animal Co., Ltd. (Shanghai, China). The purchased animals were maintained in the SPF animal laboratory of Renmin Hospital of Wuhan University. The animals were housed using standard cages in a room with humidity of 50% ± 20%, temperature of 23 °C ± 3 °C, and a 12 h light/12 h dark cycle. The animals had free access to standard laboratory food and water. All animals were fed the same feed. The experiments were conducted after the mice reached the age 18 weeks. Animal studies were approved by the Laboratory Animal Welfare and Ethics Committee of Renmin Hospital of Wuhan University (IACUC Issue No. 20171205).

Measurements of 24-h UTP, SCR, and anti-dsDNA

According to previous studies, 24-hour urine protein, serum creatinine, and serum anti-dsDNA levels have been employed to estimate the severity of lupus nephritis and therapeutic efficacy of drugs in Mrl/lpr mouse model (Feng et al., 2019; Renner et al., 2015; Wu et al., 2015b).

The urine of mice was collected for 24 h. Urine protein concentration was measured by the Bradford method, and 24-hour urine protein (24-h UTP) was calculated for LN and control mice (N = 3 for each group).

Blood samples were collected from the abdominal aorta of LN and control mice (N = 3 for each group). Serum was obtained after centrifugation of blood samples for 10 min (3,000 r/min). Serum creatinine concentration (SCR) was measured by the peroxidase method, and serum anti–double-stranded DNA antibody titers (anti-dsDNA) were determined by ELISA kits (Bioyeargene Biotechnology Co., Ltd., Wuhan) according to the manufacturer’s instructions.

Collection of kidney tissue

Mice were anesthetized by intraperitoneal injection of chloral hydrate (400 mg/kg). The mice were then sacrificed by cervical dislocation and their kidneys were removed (N = 3, for each group). The process conformed to the criteria established by the Laboratory Animal Welfare and Ethics Committee of Renmin Hospital of Wuhan University. Fresh kidney tissue was immediately chilled in liquid nitrogen and kept at −80 °C until use.

RNA extraction, library construction and sequencing

Total RNA was extracted from each kidney tissue sample using TRIzol Reagent (Ambion) according to the manufacturer’s instructions. RQ1 DNase (Promega) was used to eliminate any potential DNA in the total RNA. The quality and concentration of the purified RNA were determined by measuring absorbance at 260 nm/280 nm (A260/A280) in SmartSpec Plus (BioRad). The integrity of the RNA was detected by agarose gel electrophoresis (1.5%).

Next, 10 µg of purified total RNA from each of the samples was purified and concentrated with oligo (dT)-conjugated magnetic beads (Invitrogen) to prepare a directional RNA-seq library. The polyadenylated RNAs were ion fragmented with end repair and 5′ adapter ligation. Primers with known 3′ adapter and random hexamers were used to reverse-transcribe the fragmented RNAs. The resulting cDNAs were subsequently purified and amplified. cDNAs amplification products (200–500 bp) were retained. After quantification, the cDNAs were stored at −80 °C until sequencing.

For high-throughput sequencing, libraries were prepared according to the manufacturer’s protocols. Specifically, 150-bp paired-end sequencing was performed using an Illumina HiSeq 2000 system (ABlife. Inc, Wuhan, China).

Raw data cleaning and mapping statistics

First, we discarded the raw reads with more than 2-N bases. Adapter were clipped, and low-quality bases were removed. We further discarded short reads (<16 nt) using FASTX-Toolkit (Version 0.0.13, http://hannonlab.cshl.edu/fastx_toolkit/).

Next, TopHat2 was used to map the resulting clean reads to the mouse genome by allowing 2 mismatches (Kim et al., 2013). We discarded the reads that mapped to multiple genomic locations. The unique mapped reads were employed to calculate the FPKM value (FPKM represents reads per kilobase and per million) of each gene.

Prediction of novel LncRNA in mouse

LncRNA prediction was conducted according to a previous study (Cabili et al., 2011). In brief, Cufflinks (V2.2) was used to assemble the mapped reads with default parameters (Trapnell et al., 2012). Novel transcripts with intergenic and antisense regions were regarded as candidate lncRNAs, which should have less than 1000-bp overlap with known coding genes. The coding potential score (CPS) of candidate lncRNAs should be less than zero, which was calculated by coding potential calculator (CPC) software (Kong et al., 2007). The resulting transcripts that had more than one exon and were longer than 200 bases were defined as lncRNAs. After the antisense reads were discarded, the FPKM value of each lncRNA gene were recalculated.

Differentially expressed mRNA and lncRNA Genes (DEGs and DElncRNAs)

In this analysis, edgeR was used to identify DEGs and DElcnRNAs according to the FPKM value of each gene in each sample (Robinson, McCarthy & Smyth, 2010). DEGs and DElncRNAs were defined according to the criteria of 2-fold changes and p-values <0.05.

Coexpression network analysis

According to a previous study (Liu et al., 2017), we constructed a coexpression network between lncRNAs and mRNAs. In brief, the correlation coefficients and p-values were obtained for each mRNA-lncRNA pair based on the expression of each mRNA and lncRNA. A given threshold (absolute correlation coefficients no less than 0.7 and p-values less than 0.01) was utilized to filter the results. The filtered gene pairs form the expression network.

GO and KEGG analysis

Gene Ontology (GO) and KEGG pathway analyses were used to predict the function of DEGs, which were identified using the KOBAS 2.0 server (Xie et al., 2011). The enrichment of genes in each term or pathway was determined by the hypergeometric test and Benjamini–Hochberg FDR at corrected p-value <0.05.

RT-qPCR validation of DElncRNAs and DEGs

Quantitative real-time PCR (RT-qPCR) was performed to validate the RNA-seq data in this study. We selected some differentially expressed genes and lncRNAs with relatively high and reproducible expression for three biological replicates of each group for qPCR verification. The expression of some selected DEGs and DElncRNAs was quantified by RT-qPCR by using GAPDH as a reference. RT-qPCR was conducted on the RNA samples that were the same as those employed for RNA-seq. The RT-qPCR conditions were as follows: denaturing at 95 °C for 10 min, 40 cycles of denaturing at 95 °C for 15 s, and annealing and extension at 60 °C for 1 min. RT-qPCR of each sample employed three technical replications. The primers for each gene for RT-qPCR are listed in Table S1.

Statistical analysis

An unpaired two-tailed t-test was performed to compare the pathological evaluation and RT-qPCR data from two different groups. SPSS software for Version 19 (IBM, USA) was utilized to perform the statistical test. A significant difference was set at probability (P) values less than 0.05 (p < 0.05). Each group had at least three biological replicates (n ≥ 3). The data are presented as the mean ± standard deviation (SD).

Detection of key physiological indices demonstrated that the lupus nephritis model was established

To ensure the establishment of the LN mouse model, several biochemical indicators of LN mice were detected and compared to those of control mice. Proteinuria usually occurs after 8 weeks of age in MRL/lpr mice. We first measured the 24-hour urinary total protein (24-h UTP) of both LN and control mice. The results showed that the 24-h UTP was significantly higher in LN mice than in control mice (Fig. 1A). Serum creatinine concentration (SCR) and serum anti-double-stranded DNA antibody titers (anti-dsDNA) were also indications of LN. Compared with control mice, SCR was significantly elevated in LN mice (Fig. 1B). In addition, LN mice had a higher concentration of anti-dsDNA in serum than control mice (Fig. 1C). Taken together, these results indicted a lupus nephritis mouse model was established.

Both the mRNA and lncRNA expression profiles of kidney tissue differed between LN and control mice

To elucidate the molecular mechanisms governing lncRNAs that may account for lupus nephritis, we profiled the transcriptomes of kidney tissue from LN and control mice (N = 3 in each group). In total, six transcriptome data sets were obtained, with each composing an average of approximately 80 million PE-end reads for analysis (Table S2). Next, we aligned the clean reads to the reference sequence of mice with TopHat2 with two mismatches to detect and characterize the expression patterns of annotated genes. In addition, Cufflinks was used to perform ab initio transcript assembly and reconstructed transcripts, and 19,327 novel lncRNA genes were identified after a series of filtering steps. In total, there were 28,173 mRNA genes, and 22,504 lncRNA genes were expressed in at least one sample (with FPKM >  0) (Fig. 2A). These expressed mRNA and lncRNA genes were used to characterize the gene expression profile of LN and control mice.

Principal component analysis (PCA) was performed to explore the temporal expression patterns associated with all mRNAs and lncRNAs in our datasets. Both mRNA and lncRNA expression levels could be distinguished between the LN and control samples (Fig. 2B and 2C). Pearson correlation tests for all RNA-seq sample pairs were performed and indicated similar outcomes for mRNA expression and lncRNA expression (Fig. 2E and 2D). These results indicated that the expression profiles, including mRNAs and lncRNAs, of the kidney tissues of LN mice were significantly different from those of the control mice.

There were more upregulated lncRNAs and mRNAs than downregulated lncRNAs and mRNAs in the kidney tissue of LN mice

To explore the difference in mRNA and lncRNA expression profiles, edgeR was performed to identify the differential expression of mRNAs (DEGs) and lncRNAs (DElncRNAs) (≥2-fold up or down, FDR < 0.05) between LNs and controls. In total, there were 11,258 DElncRNAs and 5,266 DEGs (Fig. 3A). In particular, compared to controls, we identified 11,084 upregulated lncRNAs but 174 downregulated lncRNAs in LN (Fig. 3A). Notably, we also identified more unregulated mRNAs than downregulated mRNAs in mice (Fig. 3A).

Next, heatmaps were plotted based on the normalized FPKM values of DElncRNAs. Hierarchical cluster analysis showed a clear separation of LN and control samples with a high consistency for the three replicate data sets (Fig. 3B). These results were in keeping with the resutls obtained for the DEGs (Fig. 3C). Finally, the expression of five lncRNA and ten mRNA was validated by qPCR (Fig. 4). All five lncRNAs and nine out of ten mRNAs showed consistency for qPCR and RNA-seq. These results indicated that the expression at the transcript level of a larger number of mRNAs and lncRNAs was affected in LN kidney tissues in mice.

Downregulated DEGs but not upregulated DEGs were significantly enriched in immune-and SLE-associated pathways in the kidneys of LN mice

GO analysis was performed to demonstrate the potential function of all DEGs. There were 625 upregulated and 1,045 downregulated genes that were annotated with GO biological process terms. The upregulated DEGs were enriched in 180 GO terms, and the downregulated DEGs were enriched in 358 GO terms (Tables S3 and S4). Notably, the upregulated DEGs were enriched in terms related to transport and transmission (Fig. 5A). However, the downregulated DEGs were enriched in such terms as immune response and inflammatory response (Fig. 5B). Although the upregulated DEGs were also enriched in immune-and inflammatory-associated pathways, the gene number was significantly lower than that of downregulated genes (Fig. 5C). Next, we merged genes from the four immune-and inflammatory-associated pathways, including “immune system process”, “inflammatory response”, “immune response”, and “innate immune response”, which resulted in 129 downregulated genes and 17 upregulated genes (Table S5). We plotted a heatmap of these genes, which showed that the expression of these genes differed between the LN and control samples (Fig. 5D).

To further demonstrate the functional roles of these DEGs, KEGG pathway analysis was performed for the DEGs. The results showed that upregulated genes were enriched in neuroactive ligand–receptor interactions and glutamatergic synapses (Fig. 5E). The downregulated genes were enriched in such pathways as the response to Staphylococcus aureus infection and antigen processing and presentation (Fig. 5F). Notably, both downregulated and upregulated genes were enriched in systemic lupus erythematosus (Fig. 5E and 5F). However, there were 28 downregulated genes associated with systemic lupus erythematosus, but only 2 upregulated genes. A heatmap of these genes could separate the LN and control samples (Fig. 5G). In addition, approximately half of these SLE-associated genes were inflammatory factors (Fig. 5H, Table S6). These results showed that the expression of many immune and inflammatory factors was repressed in kidney tissue with LN.

LncRNAs may regulate the expression of immune and SLE-associated genes in the kidneys of LN mice

To determine the function of these DElncRNAs, we first identified the DElncRNAs with expressed mRNAs located within 10 kb upstream or downstream. There were 10,462 of 11,257 DElncRNAs with 2994 expressed mRNA genes. Next, the Pearson correlation coefficient for the expression of the DElncRNAs and mRNAs was calculated. DElncRNAs and mRNAs were considered to be coexpressed when the absolute value of the correlation coefficient was more than 0.6 with P -values <0.05. We identified 7,225,223 coexpressed pairs by 104,62 DElncRNA and 1,455 mRNA. Finally, we identified 778 mRNAs that were located within 10kb upstream or downstream of 2,181 coexpressed DElncRNAs. These mRNAs were considered to be the potential targets of the DElncRNAs (Fig. 6A). The expression of two positive and two negative pairs of DElncRNAs and mRNAs was validated by qPCR (Fig. 6B–E). This result was observed by both qPCR and RNA-seq.

Next, GO analysis was performed to annotate the function of DElncRNA-target mRNAs. The results showed that these genes were enriched in such pathways as protein ubiquitination, protein phosphorylation and the transforming growth factor beta receptor signaling pathway (Fig. 6F). In particular, these genes were also enriched in the immune system process pathway, which includes six genes (Fig. 6F). We then annotated these DElncRNA-target mRNA genes by KEGG pathway analysis. The results showed that these genes were enriched in pathways including ubiquitin-mediated proteolysis and phosphatidylinositol signaling system (Fig. 6G). Notably, one gene (H2-Aa) was annotated with systemic lupus erythematosus according to KEGG analysis (Fig. 6G). Next, we plotted the coexpression network of these SLE and immune genes and DElncRNAs. The results showed that 11 DElncRNAs were coexpressed with these genes and may have targeted these genes (Fig. 6H). qPCR confirmed that lncRNA Gm20513 positively regulated the expression of the SLE-associated gene H2-Aa in kidney tissue (Fig. 6I–L). The qPCR results were consistent with those obtained by RNA-seq. All these results indicated that some lncRNAs may regulate the expression of immune and SLE-associated genes in the kidneys of LN mice.