Integrated analysis of lncRNA and mRNA repertoires in lung from Pneumocystis infected mice

Background According to the widespread use of highly active antiretroviral therapy (HAART) to HIV patients with Pneumocystis pneumonia (PCP), the incidence of PCP associated with HIV has fallen. There is an enhancement of population of HIV-negative patients with PCP, who receiving immunosuppressive therapies or genetic primary immune deficiency disorders. Most of the previous studies have concentrated on the immune cell responses in PCP individuals, the present study aimed to integrate lncRNA and mRNA expression profiles in PCP patients without HIV. Methods The lung tissue of WT mice and BAFF-R –/– mice were harvested at 3 weeks after infected with pneumocystis. After total RNAs being extracted, transcriptome profiling was performed following the Illumina HiSeq 3000 protocol. The microarrays and quantitative RT-PCR analysis were conducted to evaluate the lncRNA and mRNA expression profiles in WT-PCP mice and BAFF-R –/– PCP mice. Results Compared with the control group, 166 mRNAs were observed to be aberrantly expressed (Fold change value≥2; P <0.05) in B cell-activating factor receptor deficient (BAFF-R –/– ) PCP group, including 99 up-regulated and 67 down-regulated, while there were 39 lncRNAs differently expressed in BAFF-R –/– PCP group, including 24 up-regulated and 15 down-regulated genes. In addition, gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed on all of the differently expressed genes and the data showed that BAFF-R deficiency played an important role in the primary and adaptive immune responses in PCP. Furthermore, the lncRNA and mRNA co-expression network was established. We noted that in the core network of lncRNA-TF (transcription factors) pairs could be classified into the categories including infection and immunity pathways. Conclusion In summary, in this study we furtherly explored the role of mature B cells in the pathogenesis and disease progression of PCP and the data demonstrated that BAFF-R deficiency play a significant role in immune regulation in PCP population. ubiquitin-like ubiquitin-protein


Introduction
Pneumocystis pneumonia (PCP) is a kind of devastating fungal disease in immunocompromised patients. Nowadays according to the widespread use of effective antiviral therapies have decreased the numbers of HIV positive patients with PCP, but there is a growing number of HIV negative PCP individuals who receiving immunosuppressive therapies (1) .
While non-HIV-infected PCP population is characterized by a more rapid disease progression and higher mortality (2,3) . Recently, the susceptibility of HIV-negative individuals for PCP has drawn the attention of researchers, but the pathogenesis and development of PCP remains to be clarified. Currently, however, the infected mice model remains the main source of laboratory studies and investigators have made improvement in understanding the host-pathogen relationship with Pneumocystis in immunology (4) . Macrophages play a significant role as the first line of host defense against Pneumocystis in the innate immune of pathogen invasion (5,6) .
In adaptive immunity, CD4 + T cells are the central factor for the clearance of Pneumocystis by the recruitment of the other effector cells (7) . Furthermore, there are some researchers verified that T helper cells including Th1, Th2, Th9 and Th17 cells could have the immune regulatory function in Pneumocystis infected mice model. Additionally, the role of B cells in PCP also has been reported in a few studies and our previous study has indicated that IL-1beta is a likely determinant of the IL-10-producing B cell-mediated suppression of Th1/Th17-cell immune responses in B cell-activating factor receptor deficient (BAFF-R -/-) PCP mice. It suggested that B cells play an important role in the regulation of Th cells in response to Pneumocystis infection (8) .
Noncoding RNAs (NcRNAs) were characterized as the deficiency of clear potential to encode proteins or peptides, while they could affect the expression of other genes via a variety of mechanisms (9). Researches over the last decades discovered many classes of ncRNA, including microRNA(miRNA), small nucleolar RNA (snoRNA), long ncRNAs (lncRNA), and circular RNA (circRNA). LncRNA are classified into ncRNA ＞ 200bp and have been demonstrated to play critical roles in the modification of genes and post-gene regulation mechanisms (10) . In recent decades, studies of lncRNA suggest that lncRNA acting as the crucial and diverse roles in the development and progression of various diseases (11) . While the function of lncRNA in the pathogenesis of PCP has not been clearly clarified yet. Our previous study has verified that in Pneumocystis infected mice, there were a great quantity of differentially expressed genes (DEGs) between WT mice and BAFF-R -/mice and IL-1bete related genes could play an immune regulatory role in Th1/Th17 cell immune responses (8) .
In the present study, RNA sequencing was performed to explore the lncRNA and mRNA expression profiles in lung tissue of WT PCP mice and BAFF-R -/-PCP mice. Additionally, to identify the underlying biological functions of differentially expressed RNAs, bioinformation analysis including GO and KEGG pathway analysis and lncRNA-mRNA co-expression network construction was conducted.

Mice
Female 6-to 8-wk-old C57BL/6 mice (Vital River Laboratory Animal Co., Ltd, Beijing, China) and BAFF-R -/mice (stock no. 007212, the Jackson Laboratory, Bar Harbor, ME) were used for this study. Mice were bred in our animal facility and kept in sterile isolated ventilated cages. All animal experiments were approved by the Capital Medical University Animal Care and Use Committee.

Infection with Pneumocystis
Pneumocystis murina was maintained in CB17 SCID mice as previously described with some modifications (8) . Briefly, mice were intratracheally inoculated with Pneumocystis cysts and control mice were administered with 100ul PBS. Mice survived at 3 weeks after infection and we obtained middle lobe from the right lung of each mouse.

Transcriptomics analysis of lung tissue
Transcriptome profiling was performed by RiboBio following the Illumina HiSeq 3000 protocol (Guangzhou, China). Firstly, total RNA was extracted from mixed samples of mice using TRIzol reagent (Invitrogen, Carlsbad, CA). Subsequently, ribosomes kit was used to remove rRNA and then mRNA was purified and fragmented into small pieces with elute, prime and fragment mix. After end-repair and adapter-ligation, the products were amplified through PCR and purified to construct the cDNA library.
The original sequencing data should be filtered, adaptor removed and low-quality reads processed to acquire clean data. For the detected mRNA, the gene expression was further calculated and the differentially expressed genes (DEGs) among the samples were analyzed.
After that, DEGs were performed for Gene Ontology (GO) analysis and Kyoto encyclopedia of genes and genomes (KEGG) pathway analysis. The expression of lncRNA was calculated and the differentially expressed lncRNA was analyzed among the samples, and specific annotations were given for the identification of the corresponding ORF, protein domain and coding potential, as well as the prediction of secondary structure and family analysis.

LncRNA target gene prediction
Target genes of lncRNA were predicted by establishing the relationship between lncRNA and mRNA and analyzing the mode of action of lncRNA. A large number of lncRNA play significant biological roles by regulating contiguous gene by integrating differentially expressed lncRNA with mRNA of their neighbors (10KB), thus obtained potential target genes of lncRNA. After the interaction between lncRNA and target genes was obtained, the network diagram of the interaction between lncRNA and target genes were drawn.

qPCR validation
RNAs were purified and quantified as previously described. RT-PCR was executed by SYBR PreMix Ex Taq TM II ROX (Takara) and the ABI Prism 7500 Sequence Detection System (Applied Biosystems, Foster City, CA). Gene expression was calculated using the 2 -delta delta CT method. GAPDH was used as a reference gene to analyze the gene quantitatively.

Statistical analysis
Data of this study was expressed as mean ± SE. Comparison of continuous data was performed by independent Student's t-test. The Prism 8 (GraphPad Software, San Diego, CA) was used for statistical analysis. P-value＜0.05 were considered statistically significant.

Identification of differentially expressed lncRNA and mRNA
A total of 6 samples from lung tissue of pneumocystis infected mice were used for this study. Removing reads with non-canonical letters or with low quality, and discarding the sequences shorter than 18 nt, the mean numbers of clean reads were 117,007,330. In total, 27343 transcripts were detected in the present study, in which 16344 transcripts were identified as mRNA and 5328 transcripts as lncRNA. The top 10 differentially expressed genes (DEGs) are shown in Table 1 (lncRNAs) and  Furthermore, we also found 3 novel lncRNAs increased in BAFF-R -/-PCP mice (Fig.3A).
On the analysis of the species on the basis of the known genetic model, using all the sample sequence data to forecast the new transcription area (FPKM > 1), and the expression level of new transcription area were analyzed (Fig.3B).

GO and pathway analysis of differentially expressed lncRNA-targeted mRNA
Since most of the functional lncRNA are mainly executed in protein-coding target genes, we applied the co-expression patterns of mRNA and lncRNA to further study the biological function of these differentially expressed lncRNA. Using GO and KEGG pathway analysis, we aimed to identify the speculated functions about the differentially expressed lncRNAs' target genes between WT-PCP group and BAFF-R -/-PCP group.
GO analysis results demonstrated that changes in biological processes of DEGs were significantly enriched in telomerase RNA binding, actin binding, cytokine receptor activity, ubiquitin-like protein transferase activity, ubiquitin-protein transferase activity, peptide antigen binding, endodeoxyribonuclease activity, drug binding, deoxyribonuclease activity and K63linked polyubiquitin biding. Changes in cellular component of DEGs were mainly enriched in nuclear DNA-directed RNA polymerase complex, catalytic complex, spindle, RNA polymerase complex, nuclear lumen, nuclear part, histone pre-mRNA 3-end processing complex, mast cell granule, transferase complex and telomerase holoenzyme complex. Changes in molecular function of DEGs were significantly enriched in regulation of organelle organization, interstrand cross-link repair, DNA repair, regulation of vascular smooth muscle cell proliferation, vascular smooth muscle cell proliferation, organelle organization, L-amino acid import, regulation of telomere maintenance, telomere organization and telomere maintenance ( Fig.6). KEGG pathway analysis demonstrated that the DEGs were mainly enriched in Chemokine signaling pathway, HTLV-I infection and endocytosis (Fig.7).

Co-expression analysis of lncRNA and mRNA
We performed a co-expression network of lncRNA and mRNA to infer the underlying regulatory function and potential target genes of lncRNA. LncRNA which locate within 300K windows downstream or upstream of the given mRNAs were subsequently further studied.
Correlated protein coding genes and lncRNA were chosen to visualize the correlation of clusters (Fig.8a). Furthermore, we analyzed the functional regulatory relationship between the genes. Functional enrichment in each cluster indicated that lncRNA were functionally associated with mRNAs in particular biological processes, such as B cell receptor signaling pathway, Herpes simplex infection, Staphylococcus aureus infection, Choline metabolism in cancer, Measles, Ubiquitin mediated proteolysis, Fanconi anemia pathway and Endocytosis (Fig.8b).

qRT-PCR validation
To validate the credibility of the results of microarray analysis, we performed real-time

Discussion
Pneumocystis pneumonia (PCP) is an opportunistic, continuously evolving disease challenging clinicians who have a life-threatening illness that inhibits their immune systems (12) .
PCP is initially observed in acquired immunodeficiency syndrome (AIDS) patients, however, because of the new antiretroviral medications recently approved in HIV positive patients, pneumocystis has been prevalence declines in these individuals (13) . Accumulating evidence revealed that the increasing use of immunosuppressive agents in haematological malignancy, autoimmune conditions and transplant recipients has resulted in an increase in PCP in HIVnegative populations (12) . According to the previous study, HIV-negative patients with PCP are related to more severe outcomes and rapid disease deteriorate (3) . Multiple immune cells such as macrophages, CD4 + T cells, CD8 + T cells and B cells participate in the development of PCP.
BAFF-R is early expressed in circulating B cells and play an important role in the maturation of B cells (14) . Using BAFF-R -/mice, our previous study focused on the function of B cells and the results showed that mature B cells play a vital role in the regulation of Th cells in PCP. We also performed mRNA expression profile analysis of WT-PCP mice and BAFF-R -/mice at 2wk post-infection and found that there were multiple differently expressed genes between the two groups. Also, GO pathway analysis indicated that a number of genes associated with immune modulatory pathways (8) .
To the best of our knowledge, there are few studies using microarray technology to reveal the immune regulatory function of cells in PCP. Advances in genetic characterization of pneumocystis indicates that specific polymorphisms are related to the epidemiology of this pathogen, containing the geographic distribution, medicine resistance, toxicity of several genotypes and genetics of this population (4) . To clarified the characteristics of the genetic diversity of PC, Magdalena et al used the sputum specimens of patients with PCP, and the results of the study showed that specific genes such as mtLSUrRNA, CYB, DHPS and SOD were correlated with certain patient characteristics (15) . Using RNA sequencing, some scholars found that Arp9, Sp and Gsc1 were differentially regulated in the predominant life forms in PCP population, providing a preliminary annotation of the Pneumocystis murina genome (16) .
Also, a study focused on the antifungal agents of PCP individuals and noted that there were 80 genes up-or downregulated among uninfected, untreated mice and mice treated with anidulafungin to study the metabolism of persisting forms (17) .
Noncoding RNAs are defined as short noncoding RNAs or lncRNA based on the sequencelength cutoff of 200 nucleotides (9) . The biological roles of lncRNA are recently beginning to be investigated, and critical functions for lncRNA have been confirmed in almost every biological system studied. Up to now, the important function of lncRNA are mainly studied in malignant diseases and recent studies showed that lncRNA play important roles in different phases of cancer immunity via molecular mechanisms (18) . However, the significance of The above genes are mostly related to immune regulatory systems: Nppa plays an important role in the control of extracellular fluid volume and electrolyte homeostasis and participates in the pathogenesis and disease progression of cardiac diseases (19) . Shisa8 is a kind of the cystineknot AMPA receptor-modulating proteins, which could mediate most of the fast-excitatory transmission in the central nervous system of vertebrates (20) . Trim10 is play a critical role in the globin gene transcription and its upregulation is a key factor required for terminal erythroid cell differentiation (21) . Ocstamp plays a key role in bone resorption and its partners CD9 could promote periodontal bone destruction by up-regulation of OC-genesis (22) . Pou2af1, which regulates B cell homeostasis and controls humoral immune, is proved to enhance IL-17A expression in Th17 cells via RORgamat (23) . Pax5 target genes are mostly expressed in leukemic cells and have been demonstrated as a tumor suppressor in lymphoblastic leukemia. (24) . This suggests that BAFF-R deficiency has an effect on the expression of a number of protein coding genes. We also performed q-PCR experiment of the DEGs to further confirmed the results of gene sequencing and the data could support the microarray data. The different express of these genes could influence the modulation of immune cells and signaling pathways. showed that differentially expressed mRNAs were mainly involved in the molecular function terms of regulation of organelle organization, telomere organization and telomere maintenance and biological process terms of cytokine receptor activity and peptide antigen binding. Also KEGG pathway analysis results indicated that pathways were principally enriched in chemokine signaling pathway, AMPK signaling pathway, etc. According to these data, we found that the deficient of BAFF-R could have a significant influence on the immune microenvironment, thus might lead to the immunosuppressive patients' susceptible to PCP.
Also, non-HIV PCP patients demonstrated the worsen disease progression, however, the treatment of non-HIV PCP has almost only concentrated on TMP-SMZ and caspofungin (25) .
The other more efficient therapy is essentially required to be explored. The differentially

Ethics approval and consent to participate
This study was approved by the ethics committee of Beijing Chao-Yang Hospital, Capital Medical University. All methods were carried out in accordance with relevant guidelines and regulations.

Consent for publication
The Author confirms: • that the work described has not been published before (except in the form of an abstract or as part of a published lecture, review, or thesis); • that it is not under consideration for publication elsewhere; • that its publication has been approved by the responsible authorities at the institution where the work is carried out.

Availability of data and materials
The datasets used or analyzed during the current study are available from the corresponding author on reasonable request.

Competing interests
The authors declare no competing non-financial/financial interests.