Construction of an lncRNA–mRNA Co-Expression Regulatory Network Mediating In ammation and Regeneration Following Acute Pancreatitis Injury


 Long noncoding RNAs (lncRNAs) play important roles in the occurrence and development of many diseases and can be used as targets for diagnosis and treatment. However, the expression and function of lncRNAs in the injury and repair of acute pancreatitis (AP) are unclear. To decipher lncRNAs’ regulatory roles in AP, we reanalyzed an RNA-seq dataset of 24 pancreatic tissues, including those of normal control mice (BL), those 7 days after mild AP (D7), and those 14 days after mild AP (D14). The results showed significant differences in lncRNA and mRNA expression of D7/D14 groups compared with the control group. Co-expression analysis showed that differentially expressed (DE) lncRNAs were closely related to immunity- and inflammation-related pathways by trans-regulating mRNA expression. The lncRNA–mRNA network showed that the lncRNAs Dancer, Gmm20488, Terc, Snhg3, and Snhg20 were significantly correlated with AP pathogenesis. WGCNA and cis regulation analysis also showed that AP repair-associated lncRNAs were correlated with extracellular and inflammation-related genes, which affect the repair and regeneration of pancreatic injury after AP. In conclusion, the systemic dysregulation of lncRNAs is strongly involved in remodeling AP’s gene expression regulatory network, and the lncRNA–mRNA expression network could identify targets for AP treatment and damage repair.


Introduction
Acute pancreatitis (AP) is an unpredictable and potentially fatal disease. Depending on its severity, it can be divided into mild acute pancreatitis (MAP), moderately severe acute pancreatitis (MSAP), and severe acute pancreatitis (SAP) according to the presence or absence of organ failure and its duration. Most AP cases clinically manifest as MAP and are self-limiting [1].
Pancreatic self-digestion causing unfolded protein reaction and the functions of autophagy, in ammatory mediator release, apoptosis, and other mechanisms have been shown to be involved in the occurrence and development of AP. Some potential molecular therapeutic targets of AP have been found, which are signi cant for the treatment of pancreatitis [2].
The endocrine and exocrine cells in the adult pancreas are not static, and can change their differentiation state and repair damaged pancreatic tissue during injury or stress response [3]. In human and animal models, pancreatic injury caused by pancreatitis can be repaired through the regeneration of acinar cells producing digestive enzymes, suggesting that revealing the molecular basis of pancreatic repair after AP may provide a new therapeutic target for the treatment and prevention of AP [4].
Long noncoding RNAs (lncRNAs) are a class of complex long-chain noncoding RNAs with little or no protein-coding ability, and are transcribed in most eukaryotic genomes [5]. Most lncRNAs have now been classi ed, but their functions require further study. Previous studies have shown that lncRNAs play an important role in the development of many human diseases [6][7][8].
Some studies have identi ed lncRNAs that are abnormally expressed in the process of AP, and found that some lncRNAs affect the occurrence and development of pancreatitis by targeting microRNAs or interacting with proteins, which are potential therapeutic targets [9][10][11][12][13]. However, the global expression pro le of lncRNAs in pancreatic tissue and their potential functions and mechanisms during repair after AP injury remain unclear. Therefore, there is a need for ongoing study of lncRNA-mediated pathogenesis and the mechanism of recovery from AP in order to formulate more appropriate treatment plans, prevent the aggravation of in ammation, and reduce the incidence of complications.
Our goal in this study is to analyze the DE lncRNAs and mRNAs expressed in a map based on APassociated RNA sequencing data, construct a novel co-expression network between DE lncRNAs and DE mRNAs, and predict the target genes and functions of trans-and cis-acting lncRNAs. The obtained ndings should contribute to a comprehensive understanding of the role of lncRNAs in AP and provide potential candidate lncRNAs for the diagnosis and treatment of AP.

Results
1. Analysis of lncRNA and mRNA expression pro le of pancreatic tissue from mouse with AP We extracted RNA sequencing data of pancreatic tissue samples [14]. The lncRNA prediction and functional analysis pipeline was performed using previously reported methods [19]. Figure 1A shows the pipeline of predicted and annotated lncRNAs based on the analytical results. A total of 199 new lncRNA genes were found, of which 156 could be compared to the noncoding RNA database, and 4,898 annotated lncRNAs were found. We then analyzed the overlapping and speci cally detected lncRNAs among the three groups, showing that a large number of lncRNAs were detected in all three groups ( Figure 1B). Based on the analysis of the lncRNAs differentially expressed among the three groups, there were 147 DE lncRNAs between the D7 group and BL group (D7 vs. BL), of which 86 were upregulated and 61 were downregulated. There were only 57 DE lncRNAs between the D14 group and BL group (D14 vs. BL), of which 36 were upregulated and 21 were downregulated ( Figure 1C).
To further explore the expression features of the DE lncRNAs from these three groups, we performed principal component analysis (PCA) for lncRNAs ( Figure 1D, E). According to the expression of all identi ed DE lncRNA genes, sample cluster analysis was carried out. We found that the samples of BL, D7, and D14 could be completely separated. Based on the expressed lncRNAs, we found D14, D7, and BL speci cations could be separated, so they may be involved in the process of injury recovery after AP ( Figure 1F). The heat map shows the expression pro le of the top 20 DE lncRNAs with the highest average expression ( Figure 1G).
The mRNA expression pro le of pancreatic tissue from mouse with AP is presented in Supplementary  Figure 1. Gene differential expression analysis showed that there were 1,425 differentially expressed genes between the D7 group and BL group (D7 vs. BL), including 1085 upregulated genes and 340 downregulated genes. There were only 504 DE genes between the D14 group and BL group (D14 vs. BL), including 406 upregulated genes and 98 downregulated genes (FC ≥ 1.5 or ≤ 0.67, P value < 0.05) ( Figure S1A). There were also found to be 252 common differentially upregulated and 40 common differentially downregulated genes between the D7 vs. BL and D14 vs. BL comparisons. The results showed that, with the recovery from AP, the gene expression in pancreatic tissue gradually recovered, but there was still differential expression ( Figure S1D, E). To further explore the sample clustering among the three groups, we performed PCA on the D7 group, D14 group, and BL group according to the expression levels of all detected genes. The results showed that the con dence ellipsis of the three groups did not overlap, and the samples could be signi cantly separated ( Figure S1B). The expression heat map of DE mRNA showed high consistency in these groups ( Figure S1C). These results showed that mRNA expression was active and differed among the D7, D14, and BL control groups, indicating that mRNA may play an important role in the recovery from AP. GO functional pathway analysis was carried out for the genes abnormally expressed in the D7 and D14 groups. The signal pathways signi cantly associated with the genes expressed in the D7 group include immune system process, innate immune response, in ammatory response, positive regulation of tumor necrosis factor production, tissue of extracellular matrix, bacteria, collagen bril tissue, proteolysis, and cell adhesion. The genes upregulated in the D14 group were particularly associated with acute phase reaction, peptide hormone reaction, sodium ion transport, collagen bril tissue, positive regulation of phosphatidylinositol 3-kinase signal, ossi cation, wound healing, steroid metabolism, ion transport, and cell adhesion, among others. The signal pathways signi cantly associated with the genes downregulated in the D7 group include ion transport, transmembrane transport, regulation of ion transmembrane transport, potassium ion transport, potassium ion transmembrane transport, exercise behavior, response to hypoxia, lipid decomposition process, protein complex assembly, and lipid metabolism process. In the D14 group, the positive regulation of downregulated gene enrichment gene expression, transmembrane transport cell differentiation. Further GO function analysis was conducted on the genes that were upregulated and downregulated. The genes that were upregulated and downregulated were particularly associated with collagen bril tissue, proteolysis, ossi cation, wound healing, cell adhesion, mechanical stimulation response, positive regulation of cell population proliferation, extracellular matrix tissue, lung development, cell oxidative detoxi cation, and downregulated genes. Owing to the small number of genes, the relevant functional pathway results could not be output in the GO analysis (Figure S1F-I).
2. Analysis of correlation of lncRNA and mRNA expression in pancreatic tissue from mouse with AP By analyzing the correlation of expression of DE lncRNAs and mRNAs in different samples, it was found that 19 DE lncRNAs in the D7 vs. BL comparison were correlated with 243 DE mRNAs, while only 5 DE lncRNAs in the D14 vs. BL comparison were correlated with 55 DE mRNAs (Figure 2A, B). Upon further overlapping analysis on the known differences of upregulation and downregulation in D7 vs. BL and D14 vs. BL (removing the newly predicted RNA), it was found that there were 16 upregulated lncRNA genes and 5 downregulated lncRNA genes ( Figure 2C, D). Bar plots show the top 10 enriched KEGG pathways of the co-expressed mRNAs of up and down overlapped DE lncRNAs (Supplementary Figure 2A, B). Further functional analysis was conducted on 21 DE lncRNA genes from the D7 vs. BL and D14 vs. BL comparisons and the mRNA genes related to the expression of lncRNA genes speci cally differentially expressed in the D7 vs. BL comparison. It was found that the mRNA genes related to the 21 overlapping DE lncRNA genes were particularly associated with extracellular matrix recombination, while the genes co-expressed with lncRNAs speci cally differentially expressed in the D7 vs. BL comparison were particularly associated with immunity-and in ammation-related pathways, suggesting that there may be a certain in ammatory reaction in 7 days ( Figure 2E, F). By constructing a regulatory network diagram, it was found that Dancr and Gm20488 may affect the tissue recovery and regeneration after AP by affecting the expression of extracellular matrix-related genes, while TERC, Snhg3, and Snhg20 may participate in the in ammatory response during the recovery period by affecting the expression of in ammation-related genes ( Figure 2G, H).
3. Weighted gene co-expression network analysis (WGCNA) of lncRNAs and mRNAs in pancreatic tissue from mouse with AP Through the co-expression analysis of all DE lncRNA and mRNA genes by WGCNA, ve co-expression modules were found. Among them, lncRNA and mRNA genes in the blue, turquoise, yellow, and gray modules were signi cantly correlated with at least one treatment group factor, and the number of differentially expressed lncRNAs and mRNAs varied greatly among the different modules ( Figure 3A-C). Through functional enrichment analysis of genes in the blue and turquoise modules with a large number of differentially expressed mRNAs, it was found that the genes in blue were particularly associated with immunity-and in ammation-related pathways, while the turquoise module was also associated with extracellular matrix recombination, wound recovery and other pathways, in addition to immunity-and in ammation-related pathways, indicating that there are certain differences in the functions of lncRNAs and mRNAs in the different modules ( Figure 3D

Cis-target analysis of DE lncRNAs in pancreatic tissue from mouse with AP
For the DE lncRNAs, we extracted the information on genes within 10 kb upstream and downstream, calculated the Pearson's correlation coe cient between DE lncRNAs and mRNAs by co-expression analysis, and screened the lncRNA-target relationship pairs that met the criteria of correlation coe cient > 0.6 and P value < 0.05. Then, the two datasets were combined to obtain the cis targets of the DE lncRNAs ( Figure 4A). To further reveal the function of DE lncRNAs in AP, the cis regulatory targets of lncRNAs were predicted and analyzed. It was found that 18 DE lncRNAs may affect the expression of 25 mRNAs in a cis regulatory manner. There was a certain correlation between the expression of these lncRNAs and cis regulatory mRNA genes in the D7 vs. BL and D14 vs BL comparisons. The lncRNA Gm28177 may regulate the expression of STAT4 and STAT1; the lncRNA 9630013d21rik may regulate the expression of Cdkn2c, and the lncRNA Mir217hg may regulate the expression of Mir217 and Mir216a ( Figure 4B-G).
KEGG function analysis of 25 mRNAs cis-regulated by lncRNAs showed that they were particularly associated with in ammatory bowel disease (IBD), Th1 and Th2 cell differentiation, hepatitis B, JAK-STAT signal pathway, and necrotic apoptosis, among which STAT4 and STAT1 involved in the in ammation-related pathway and Cdkn2c in the cell cycle pathway warrant particular attention ( Figure  4H-J).

Discussion
Understanding the molecular mechanism of self-repair and regeneration after AP is helpful for its diagnosis and treatment. To date, few studies have focused on the relationship between lncRNAs and the repair/regeneration of AP, and the role of lncRNAs in AP has not been fully clari ed.
In the current study, RNA-seq data analysis showed that a large number of genes were differentially expressed in the pancreas of D7 and D14 groups compared with the expression levels in the BL group. GO functional enrichment analysis showed that they have a variety of different functions. At the stage of pancreatic recovery from 7 to 14 days, a large number of repair-related genes such as those involved in wound repair, cell proliferation, and matrix recombination were continuously expressed to promote the repair and regeneration of the pancreas.
At the same time, we found a total of 5,097 lncRNAs from the RNA-seq data. One hundred and forty-seven lncRNAs were differentially expressed between the D7 group and BL group, while only 57 lncRNAs were differentially expressed between the D14 group and BL group. The results showed that, with the recovery of AP, the expression of lncRNA genes in pancreatic tissue gradually recovered like that of mRNA genes, but differential expression still occurred. Based on the expression of these DE lncRNAs in various samples, we focused on the top different lncRNAs, investigated their functions, analyzed their potential involvement in pancreatitis, and found that differentially expressed lncRNAs such as Snhg3, Dancr, TERC, and Snhg20 warrant particular attention. Snhg3 can participate in the activation of microglia and then affect cerebral ischemia-reperfusion injury [20]. Meanwhile, Dancr promotes the proliferation and metastasis of pancreatic cancer by regulating miRNA-33b [21]. TERC alleviates the progression of osteoporosis through the uptake of miRNA-217 by upregulating RUNX2 [22]. Finally, Snhg20 participates in angiotensin II-induced cardiac brosis and hypertrophy [23].
WGCNA has been used to reveal the expression relationships between protein-coding genes and lncRNAs [24]. Therefore, the differentially expressed lncRNAs identi ed in this study may play important roles in the process of recovery after pancreatitis and warrant further study. By constructing a regulatory network diagram, it was found that Dancr and Gm20488 may affect the expression of extracellular matrix-related genes, while TERC, Snhg3, and Snhg20 may affect the expression of in ammation-related genes. Therefore, these lncRNAs may affect the repair and regeneration of pancreatic injury after pancreatitis by regulating the expression of extracellular matrix and in ammatory response genes. lncRNAs can affect the expression of distant genes in a cis manner. Several AP-associated lncRNAs have been demonstrated to directly bind target proteins or mRNAs and result in post-transcriptional modi cation [25]. Therefore, we predicted the function of lncRNAs by analyzing cis-acting and coexpression networks. We found that 18 DE lncRNAs may affect the expression of 25 mRNA genes through a cis regulatory mechanism. There was a certain correlation between the expression of these lncRNAs and cis regulatory mRNA genes in the D7 vs. BL and D14 vs BL comparisons. We investigated the literature related to the lncRNAs mRNA gene encoding above. The genomic loci of miR-216a and miR-217 are located in Mir217HG, which is involved in regulating the formation of pancreatic ducts, but does not affect the development of pancreatic cancer [26]. Meanwhile, the Stat4 and Stat1 signal pathway is involved in regulating pancreatitis and pancreatitis lung injury [27][28][29]. Moreover, Cdkn2c may be involved in regulating the formation of pancreatic endocrine tumors [29]. However, at present, there is no literature on lncRNAs such as Gm28177 and 9630013d21rik. It is speculated that they may affect the development of pancreatitis and the recovery from injury due to it via cis regulation of the expression of the above mRNA genes.
In conclusion, this study identi ed DE lncRNA and mRNA genes in pancreatic tissues of mice in the D7, D14, and BL groups, established the expression correlation and co-expression relationship between DE lncRNAs and mRNA genes, and analyzed the cis regulatory targets of lncRNAs. The results provide clues for further verifying and revealing the function and mechanism of action of lncRNAs in the repair and regeneration of pancreatic injury in AP, and are of great signi cance for identifying key targets in the treatment and repair of AP.

Material And Methods
Retrieval and processing of publicly available data The publicly available sequence data les of GSE99774 deposited by Boggs et al. [14] were downloaded from the Sequence Read Archive (SRA). The GSE99774 dataset contains RNA sequencing data of 24 pancreatic tissues, 8 in each group. A mouse model of mild AP was induced by cerulein. SRA Run les were converted to fastq format with NCBI SRA Tool fastq-dump. The raw reads were trimmed of lowquality bases using FASTX-Toolkit (v.0.0.13; http://hannonlab.cshl.edu/fastx_toolkit/). Then, the clean reads were evaluated using FastQC (http://www.bioinformatics. babraham. ac.uk/projects/fastqc).

Read alignment and differentially expressed gene (DEG) analysis
Clean reads were aligned to the human GRCh38 genome by HISAT [15] allowing no more than four mismatches. Uniquely mapped reads were ultimately used to calculate read number and reads per kilobase of exon per million fragments mapped (RPKM) for each gene. The expression levels of genes were evaluated using RPKM. For performing gene differential expression analysis, we chose the software DEseq2. DEseq2 models the original reads and uses a scale factor to explain the difference of library depth. Then, DEseq2 estimates the gene dispersion and reduces these estimates to generate more accurate dispersion estimates to model the read count. Finally, the model of a negative binomial distribution is tted by DEseq2, and the hypothesis is tested by Wald test or likelihood ratio test. DEseq2 can be used to analyze the differential expression between two or more samples, and the analytical results can be used to determine whether a gene is differentially expressed based on fold change (FC; the absolute ratio of expression change) and false discovery rate (FDR). In this study, the criteria for a signi cant difference in expression were as follows: FC ≥ 1.5 or ≤ 1.5, and P < 0.05.

lncRNA prediction and direction identi cation
To systematically analyze the lncRNA expression patterns, we used a pipeline for lncRNA identi cation similar to that reported previously [16], which was constructed based on Cu inks software [17]. All steps of the pipeline are shown in Figure 1A.

Functional enrichment analysis
To clarify the functional categories of the DEGs, Gene Ontology (GO) terms and KEGG pathways were identi ed using KOBAS 2.0 server [18]. Hypergeometric test and Benjamini-Hochberg FDR controlling procedure were used to de ne the enrichment of each term. Reactome (http://reactome.org) pathway pro ling was also used for functional enrichment analysis of the sets of selected genes.

Other statistical analyses
Principal component analysis (PCA) was performed using the R package factoextra (https://cloud.rproject.org/package=factoextra) to show the clustering of samples with the rst two components. After normalizing the reads by tags per million (TPM) of each gene in the samples, an in-house script (sogen) was used to visualize the next-generation sequence data and genomic annotations. The pheatmap package (https://cran.r-project.org/web/packages/pheatmap/index.html) in R was used to perform clustering based on Euclidean distance. Student's t-test was used for comparisons between two groups.  Expression correlation analysis of lncRNAs and mRNAs in pancreas tissue from mouse with acute pancreatitis.
A. Scatter plot showing DElncRNAs enrichment by D7 compared with BL samples and its number of coexpressed DEmRNAs. Red points denote up-regulted lncRNAs involved in co-expression pairs and blue points denote down-regulated lncRNAs. Cutoffs of pvalue < 0.01 and Pearson coe cient > 0.9 were applied to identify the co-expression pairs. B. Scatter plot showing DElncRNAs enrichment by D14 compared with BL samples and its number of coexpressed DEmRNAs. Red points denote up-regulted lncRNAs involved in co-expression pairs and blue points denote down-regulated lncRNAs. Cutoffs of pvalue < 0.01 and Pearson coe cient > 0.9 were applied to identify the co-expression pairs.   Cis-targets analysis of De-lncRNAs in pancreas tissue from mouse with acute pancreatitis.
A. For the signi cantly differentially expressed lncRNA(DElncRNA), extract the gene information within 10KB upstream and downstream, calculate the Pearson correlation coe cient between DE lncRNA and mRNA for co expression analysis, and screen the lncRNA-target relationship pairs that meet the absolute value of correlation coe cient greater than 0.6 and P < 0.05, Then the two data sets were combined to obtain the cis target of lncRNA. DElncRNA, the number of signi cantly differentially expressed lncrnas (Union set) detected in the comparison of multiple groups of the project; Expressed gene, the number of expressed genes detected in all samples of the project; Candidate lncRNA, the number of differentially expressed lncRNAs within 10KB upstream and downstream; Up / down stream mRNA, the number of genes within 10KB upstream and downstream of lncrna; Remaining lncRNA, the number of lncRNAs that meet the above location requirements and co expression relationship; cis target, the number of lncRNA cis regulatory targets that meet the above location requirements and co expression relationship, and enter the subsequent GO and KEGG analysis.