Comprehensive Gene Analysis of IgG4-Related Ophthalmic Disease Using RNA Sequencing

High-throughput RNA sequencing (RNA-seq) uses massive parallel sequencing technology, allowing the unbiased analysis of genome-wide transcription levels and tumor mutation status. Immunoglobulin G4-related ophthalmic disease (IgG4-ROD) is a fibroinflammatory disease characterized by the enlargement of the ocular adnexal tissues. We analyzed RNA expression levels via RNA-seq in the biopsy specimens of three patients diagnosed with IgG4-ROD. Mucosa-associated lymphoid tissue (MALT) lymphoma, reactive lymphoid hyperplasia (RLH), normal lacrimal gland tissue, and adjacent adipose tissue were used as the controls (n = 3 each). RNA-seq was performed using the NextSeq 500 system, and genes with |fold change| ≥ 2 and p < 0.05 relative to the controls were defined as differentially expressed genes (DEGs) in IgG4-ROD. To validate the results of RNA-seq, real-time polymerase chain reaction (PCR) was performed in 30 IgG4-ROD and 30 orbital MALT lymphoma tissue samples. RNA-seq identified 35 up-regulated genes, including matrix metallopeptidase 12 (MMP12) and secreted phosphoprotein 1 (SPP1), in IgG4-ROD tissues when compared to all the controls. Many pathways related to the immune system were included when compared to all the controls. Expressions of MMP12 and SPP1 in IgG4-ROD tissues were confirmed by real-time PCR and immunohistochemistry. In conclusion, we identified novel DEGs, including those associated with extracellular matrix degradation, fibrosis, and inflammation, in IgG4-ROD biopsy specimens. These data provide new insights into molecular pathogenetic mechanisms and may contribute to the development of new biomarkers for diagnosis and molecular targeted drugs.


Introduction
Immunoglobulin (Ig) G4-related disease (IgG4-RD) is a relatively new disease concept that was originally proposed in Japan. This disease causes fibrosis in the organs throughout the body including the eyes, and the organs become swollen [1]. Involvement of follicular T cells (Tfh) and regulatory T cells (Tregs) has been reported [2,3], but the pathogenetic mechanisms remain to be elucidated. Various etiologies, such as infectious diseases [4], have been proposed.
IgG4-related ophthalmic disease (IgG4-ROD) is an orbital lymphoproliferative disorder that show enlargement of the ocular appendages, including the lacrimal gland. IgG4-ROD is classified as a benign disease, though malignant lymphoproliferative diseases, such as mucosa-associated lymphoid tissue (MALT) lymphoma, can also occur in the orbit. Clinical symptoms, diagnostic imaging, histology, molecular analysis, and flow cytometry using biopsy specimens are the most commonly used diagnostic tools for orbital lymphoproliferative disorders [5,6]. Differentiation of these orbital lymphoproliferative disorders is often difficult, as these lesions share clinical, imaging, and histologic features, as well as molecular markers [5,[7][8][9][10]. In addition, several reports have suggested that 12% of orbital MALT lymphoma result from a background of IgG4-ROD [11,12].
Differential diagnosis between IgG4-ROD and orbital MALT lymphoma, both of B-cell lineage, remains complicated and challenging in routine clinical practice because of the lack of specific diagnostic biomarkers. Elevated serum IgG4 is not sufficiently sensitive or specific enough for this purpose [13,14]. Therefore, approximately one third of IgG4-ROD patients do not completely meet the clinical criteria and are diagnosed with 'possible' or 'probable' IgG4-ROD, which contributes to clinical confusion and delayed diagnosis [6,15]. Reactive lymphoid hyperplasia (RLH) is another benign lymphoproliferative disorder. To the best of our knowledge, there is also no report on the differentiation between IgG4-ROD and RLH based on gene profiling, making differential diagnosis difficult.
RNA sequencing (RNA-seq) is a type of next-generation sequencing (NGS) technology that allows the comprehensive evaluation and quantification of all subtypes of RNA expressed in tissues [16]. Compared with the microarray analysis conventionally used in transcriptome analysis, the correlation coefficient of the fold change (FC) of the gene expression level between RNA-seq and the microarray was approximately 0.7 [17]. Moreover, RNA-seq detected expression differences in some genes that were not significantly different in the microarray, and the differences in the expression of these genes were confirmed by a quantitative polymerase chain reaction (PCR) analysis [17]. Other advantages of RNA-seq include a greater dynamic range and the ability to identify abnormally altered genes or molecular pathways that may lead to the discovery of novel diagnostic biomarkers. Therefore, RNA-seq is an ideal tool for defining the gene expression profiles or transcriptomes of orbital lymphoproliferative disorders.
Several previous studies used NGS to analyze IgG4-RD for the purposes of elucidating the effects of steroid treatment and discovering for biomarker [18,19], but these studies did not compare IgG4-ROD and other lymphoproliferative disorders. In this study, we investigated the expression of RNA in lesions of IgG4-ROD in biopsy specimens compared with various controls comprising adipose tissues adjacent to the IgG4-ROD lesion, orbital MALT lymphoma, RLH, and normal lacrimal gland tissues to search for new biomarkers and elucidate the pathophysiology of IgG4-ROD.

Patients
This study included 71 patients with lymphoproliferative disorders who visited Tokyo Medical University Hospital of Ophthalmology between 2016 and 2020.
RNA-seq was conducted using biopsy specimens containing typical lesions obtained from IgG4-ROD patients (n = 3). As controls, biopsy specimens of adjacent adipose tissue from the same IgG4-ROD patients, orbital MALT lymphoma possibly derived from the lacrimal gland, RLH, and normal lacrimal gland tissues were used (n = 3 each). The normal lacrimal gland tissues were collected to examine them for infiltration at the time of removing other orbital tumor (schwannoma, 1 case) or for the purpose of diagnosing lacrimal gland enlargement (lacrimal gland dislocations, 2 cases). All biopsy specimens were stored frozen within 30 minutes after collection and each sample was stored at −80 • C until analysis. The demographic and laboratory features of the subjects and controls studied by RNA-seq are summarized in Table 1. The mean (±SD) serum IgG4 level was higher in the IgG4-ROD patients than in the patients with MALT lymphoma and RLH (1253.3 ± 770.3, 22.7 ± 10.1 and 53.1 ± 53.1 mg/dL, respectively; not measured in lacrimal gland cases). Real-time PCR was performed using biopsy specimens obtained from 30 patients with IgG4-ROD and 30 patients with orbital MALT lymphoma. The demographic and laboratory features of the study population in real-time PCR study are summarized in Table 2. The 30 samples of MALT lymphoma included the one sample used for RNA-seq, but the 30 samples of IgG4-ROD did not include those used in RNA-seq. The diagnoses of IgG4-ROD, orbital MALT lymphoma, and RLH were based on clinical, radiographic, histologic and flow cytometric studies, and molecular genetic analyses, such as gene rearrangement, in the biopsies. In particular, the diagnosis of IgG4-ROD was made in accordance with the published criteria [6]. Briefly, the criteria were based on three characteristic findings: (1) enlargement of the orbital tissues, marked diffuse lymphoplasmacytic infiltrate with either fibrosis or sclerosis, (2) tissue IgG4/IgG ratio > 40% and tissue IgG4 + plasma cells >50 cells/high power field, and (3) serum IgG4 level above 135 mg/dL. Definite IgG4-ROD was diagnosed when (1), (2) and (3) were present. Probable IgG4-ROD was diagnosed when (1) and (2) were present. Twenty-eight patients had definitive and five had probable IgG4-ROD.
Written informed consent was obtained from all the participants of the study. The study was also approved by the Ethics Committee of the Tokyo Medical University Hospital, Tokyo, Japan (number: SH3281). All investigations were conducted according to the principles of the Helsinki declaration.

RNA-Seq Analysis for Differentially Expressed Genes and Pathway Analysis
Total RNA was extracted from each tissue using TRIZOL (Thermo Fisher Scientific, Waltham, MA, United States of America (USA)), and RNA integrity was measured by a Bioanalyzer RNA 6000 Pico Kit (Agilent Technologies, Santa Clara, CA, USA). The depletion of ribosomal RNA (rRNA) and RNA-seq library preparation was performed with a NEBNext rRNA Depletion Kit and a NEBNext Ultra Directional RNA Library Prep Kit, respectively (New England Biolabs, Ipswich, MA, USA). Sequencing for 2 × 36-bp paired-end reads was performed with NextSeq500 (Illumina, San Diego, CA, USA). FASTQ files were imported to the CLC Genomics Workbench (Qiagen, v10.1.1, Hilden, Germany). Reads were mapped to the hg19 human reference genome and quantified against 57,773 genes. Differentially expressed genes (DEGs) in IgG4-ROD were defined as p < 0.05 and absolute FC (|FC|) ≥ 2 and were calculated using Empirical Analysis with the DEG tool.
Pathways associated with up-regulated DEGs were analyzed using Reactome [20]. Up to 3000 genes can be input into Reactome. When there were 3001 or more DEGs, the top 3000 genes according to their p values were analyzed. The Reactome provides many cellular processes as a network of molecular transformations under a single consistent data model [21].

Real-Time Quantitative PCR
The total RNA was extracted using a miRNAeasy Mini Kit (Qiagen GmbH, Hilden, Germany) and reverse transcribed using a High Capacity RNA-to-cDNA kit (Thermo Fisher Scientific) according to the manufacturer's instructions. Matrix metallopeptidase 12 (MMP12), secreted phosphoprotein 1 (SPP1), and β-actin were detected using TaqMan probes (Assay ID: Hs00159178_m1, Hs00959010_m1, and Hs99999903_m1, respectively; Thermo Fisher Scientific) with a 7900 HT Fast Real-Time PCR system (Thermo Fisher Scientific). The PCR cycling conditions were 95 • C for 10 min, followed by 40 cycles of 15 s at 95 • C and 60 s at 60 • C. Data were normalized to the expression of the β-actin. Relative quantification was performed using a calibration curve with placental cDNA.

Statistical Analysis
Real-time PCR data were analyzed with SPSS statistical software version 26. The statistical significance of difference was determined using a Mann-Whitney U test. A p value less than 0.05 was considered statistically significant. Principal component analysis (PCA) was performed using R (3.6.2.) [22]. PCA was used to discriminate the different biological samples based on the distances of a reduced set of new variables (principal components). PCA was performed using the normalized expression levels of 8856 RNAs obtained from RNA-seq for each sample. Analysis of variance (ANOVA) was also performed, and RNAs with p < 0.05 were extracted. The results of the first and second principal components were used for plotting the results in two dimensions. Unsupervised hierarchical clustering analysis using DEGs was performed using an algorithm based on Pearson correlation and the average-linkage method.

DEGs of IgG4-ROD and Other Lymphoproliferative Disorders Using RNA-Seq
RNA-seq detected 35 up-regulated DEGs and nine down-regulated DEGs in IgG4-ROD when compared with all the controls ( Figure 1). An unsupervised hierarchical cluster analysis was then performed on these 44 genes to investigate variations in the DEGs between IgG4-ROD and the other controls ( Figure 1, Table 3). These genes were segregated clearly into low and high expression groups for IgG4-ROD specimens but showed widely variable distributions with no consistent patterns in all the control tissues, indicating that the cluster analysis using DEGs separated IgG4-ROD from other diseases and normal tissues. These results indicate a possibility that these genes may contain diagnostic and/or therapeutic biomarkers or important factors contributing to the pathological condition of IgG4-ROD.  Table 3. List of DEGs when the IgG4-ROD tissue was compared to all control tissues. Table 3. Cont.
When compared to adjacent adipose tissue, 2698 DEGs were up-regulated and 3057 were down-regulated in IgG4-ROD. Pathway analysis using the up-regulated genes revealed that 59 pathways (including duplications) were related to IgG4-ROD ( Figure 2, Table 4). The predominant pathways were related to the immune system, such as the "endosomal/vacuolar pathway", "antigen presentation: folding, assembly, and peptide loading of class I MHC", "immunoregulatory interactions between a lymphoid and a non-lymphoid cell", "interferon alpha/beta signaling", and "ER-phagosome pathway".  Compared with orbital MALT lymphoma, 1718 DEGs were up-regulated and 705 were down-regulated in IgG4-ROD. Pathway analysis using the up-regulated genes revealed that 64 pathways were related to IgG4-ROD ( Figure 3, Table 5). The predominant pathways were related to the immune system, such as "classical antibody-mediated complement activation", "FCGR activation", "regulation of complement cascade", "CD22 mediated BCR regulation", and "complement cascade".   IgG4-ROD: Immunoglobulin G4 related ophthalmic disease, MALT: mucosa-associated lymph tissue.
When compared with RLH, 439 DEGs were up-regulated and 342 were down-regulated in IgG4-ROD. Pathway analysis using the up-regulated genes revealed that 43 pathways (including duplications) were related to IgG4-ROD ( Figure 4, Table 6). The predominant pathways were related to the immune system, such as "classical antibody-mediated complement activation", "FCGR activation", "CD22 mediated BCR regulation", "regulation of actin dynamics for phagocytic cup formation", and "role of phospholipids in phagocytosis". IgG4-ROD and RLH were distributed close to each other in the PCA plot ( Figure 5), suggesting a possibility that the two diseases share similarities.    In total, 2897 up-regulated and 3554 down-regulated DEGs were detected in RLH compared with normal lacrimal gland tissue; among them, 2262 and 2875 DEGs, respectively, were also detected in IgG4-ROD. Table S1 lists the down-regulated and up-regulated DEGs showing low p values in RLH compared with the lacrimal gland; all these genes showed similar dysregulated patterns in IgG4-ROD. However, RNA-seq detected 781 DEGs in IgG4-ROD compared with RLH, which could be used for further analysis. The normalized mRNA expression levels in IgG4-ROD and the other lymphoproliferative disorders obtained from RNA-seq were subject to PCA, and the plot of the first principal components versus that of the second principal components is shown in Figure 5. IgG4-ROD, MALT lymphoma, and RLH were distributed relatively close to each other, while adipose tissue and lacrimal gland tissue were distributed in distinctly different quadrants. However, it is possible to explicitly divide these groups into five groups (IgG4-ROD, adipose tissue in IgG4-ROD, MALT lymphoma, RLH, and normal lacrimal gland tissue) by the first component and second component. When compared with normal lacrimal gland tissue, 3026 DEGs were up-regulated, and 3513 were down-regulated, in IgG4-ROD. Pathway analysis using the up-regulated genes revealed that 188 pathways (including duplications) are likely related to IgG4-ROD ( Figure 6, Table 7). The predominant pathways were related to the cell cycle, such as "cell cycle", "cell cycle, mitotic", and "cell cycle checkpoints", and to the immune system, such as "cytokine signaling in immune system" and "generation of second messenger molecules".  An unsupervised hierarchical cluster analysis was then performed to investigate the variations in DEGs between IgG4-ROD and each control (Figures 2, 3, 4 and 6b). In all the cluster analyses, IgG4-ROD was clearly separated from each control.
Pathway analysis identified many pathways related to the immune system ( Table 8), suggesting that immune reaction is strongly involved in the pathology of IgG4-ROD compared with other diseases and non-disease controls. The complete results of all pathway analyses are shown in Tables S2-S5.

Immunostaining
Representative micrographs of immunostaining for MMP12 and SPP1 in IgG4-ROD, MALT lymphoma, RLH, and lacrimal gland tissues are shown in Figure 8. In the IgG4-ROD tissue, many MMP12-positive cells were present not only in lymphoid follicles but also in fibrotic areas, and SPP1 was intensely immunostained in the fibrotic part. On the other hand, MMP12 and SPPI reactivities were weak or absent in the other three tissue samples. These findings suggest that MMP12 is involved in fibrosis and follicle formation, while SPP1 is involved in fibrosis. Immunostaining indicated expression of MMP12 and SPP1 at the protein level in IgG4-ROD biopsy tissues, suggesting that these molecules play a role in pathogenesis and could thus serve as new biomarkers for IgG4-ROD.

Discussion
In this study, we performed a comprehensive RNA analysis of lacrimal gland-derived IgG4-ROD using adipose tissue adjacent to the IgG4-ROD lesion, orbital MALT lymphoma, RLH, and normal lacrimal gland tissue as controls. Orbital MALT lymphoma and RLH were disease controls, while adipose tissue from the same patients with IgG4-ROD controlled for individual differences in gene expression, and normal lacrimal gland tissue controlled for tissue-specific gene expression. The results of the RNA-seq and pathway analysis suggest that IgG4-ROD is strongly related to the immune system and the up-regulation of MMP12 and SPP1.
The DEGs and pathways detected by RNA-seq for IgG4-ROD versus the lacrimal gland may provide insight into the pathology of IgG4-ROD. Fibrosis and class switches to IgG4 are characteristics of the pathology of IgG4-RD. Fibrosis in IgG4-RD is considered to be caused by transforming growth factor-β (TGF-β) secreted by Tregs and CD4 + cytotoxic T cells (CD4 + CTLs) [3]. CD68 + CD163 + alternatively activated (M2) macrophages are abundantly present in IgG4-RD tissues and are involved in the orchestrated immune reaction by regulating cytokine production and fibrosis [23][24][25].
Moreover, binding of the dendritic cells and T cells to B cells via BAFF and CD40 together with the activation of AID by various cytokines has been reported to induce class switching [26]. Furthermore, class switching to IgG4 may be induced by the interleukin (IL)-21 secreted by Tfh in IgG4-RD [27]. Compared to lacrimal tissue, not only TGFB1 but also genes related to the surface markers of cells secreting TGF-β (CD4 + CTL marker, CD4; Treg markers, CD25 and FOXP3), as well as multiple genes related to class switching (BAFF, BAFFR, CD40, CD40L, AID, IL21) and Tfh (ICOS, CXCR5), were up-regulated in IgG4-ROD. Thus, the DEGs detected by RNA-seq demonstrate the existence of multiple genes that are implicated in the pathology of IgG4-RD in previous reports, as well as novel DEGs that should be further investigated.
In the pathway analysis using DEGs detected in IgG4-ROD versus orbital MALT lymphoma, the pathways related to the complement system dominated. Some patients with IgG4-RD have hypocomplementemia [28], suggesting the occurrence of antigen-antibody reactions in the body. Since Chlamydia psittaci DNA has been found in MALT lymphoma samples [29], antigen-antibody reactions to infection may also take place in MALT lymphoma. However, infection only serves as a trigger, and the leading cause of the pathology seems to be the self-proliferation of tumor cells due to MALT1 translocation under transcriptional control in the IgH enhancer region [30,31], activated by infection. In IgG4-ROD, however, the complement pathway is activated, and the complement is depleted. Thus, the leading cause of the pathology is likely antigen-antibody reactions to infections [4]. Since MALT lymphoma is speculated to arise from a background of IgG4-RD [11,12], both diseases may have infection as a common background.
Notably, IgG4 is regarded as an anti-inflammatory antibody because of its poor capacity to bind to Fc receptors [32], Therefore, one possible mechanism for activation of the complete pathway is that other IgG subclasses, such as IgG1 and IgG2a, bind to Fc receptors. A similar mechanism was proposed to explain why hypocomplementemia occurs frequently in IgG4-ROD patients despite the relative inability of IgG4 to fix the complement [33]. In this study, FCGR1A was up-regulated in IgG4-ROD. This gene encodes a receptor with high affinity for the Fc region of IgG, suggesting that the complement's pathway may be activated by binding to a subclass other than IgG4.
RLH is part of the spectrum of orbital lymphoproliferative disorders [34] and is considered to be a consequence of the chronic inflammatory response of lymphoid cells to irritating or antigenic stimuli [35,36]. Differentiation between IgG4-ROD and RLH is often difficult in the results of histology, gene rearrangement, and flow cytometry [5]. Indeed, fewer DEGs were detected in IgG4-ROD compared to RLH than with the other control tissues. Furthermore, PCA using the normalized expression levels of the 8856 RNAs obtained from RNA-seq showed a close distribution of IgG4-ROD and RLH in terms of their genetic features, which may reflect the similar clinical phenotypes of the two conditions. This finding suggests that the altered genes are almost the same in both diseases. Indeed, some genes related to autoimmune diseases, such as BACH2, C17orf99, CCL19, and IL6R, showed high expression in IgG4-ROD when compared to tissues other than RLH but no difference when compared to RLH, suggesting that the two diseases may involve an autoimmune etiology [35]. If this hypothesis holds true, it is not surprising that there is considerable overlap in the gene alterations between the two conditions. An interesting question is how IgG4-ROD differs from RLH at the molecular level. To the best of our knowledge, there is no report comparing RLH with IgG4-ROD. In this study, we found for the first time that RNA-seq was able to identify 781 DEGs for IgG4-ROD versus RLH. Until now, IgG4 was the only factor distinguishing between the two, and it was difficult to search for other biomarkers. RNA-seq detected 781 DEGs, indicating the possibility of searching for new biomarkers. When we performed a preliminary pathway analysis using genes that were up-regulated in RLH relative to IgG4-ROD, pathways including "collagen degradation" were identified (Table S6). In addition, the increased expression of BMP2 with anti-fibrotic action in RLH has been reported [37]. These findings suggest that anti-fibrosis may be promoted in RLH (or attenuated in IgG-ROD) with differential biomarkers between IgG4-ROD and RLH.
When common DEGs were extracted while excluding adjacent adipose tissue, 96 genes remained as DEGs. The genes obtained from this analysis are putatively involved in diseases involving IgG4 and IgE. However, there was no difference in the expression of these genes in comparison with the adjacent adipose tissue. This suggests that the surrounding tissues may also be affected by IgG4-ROD. However, the adipose tissue was a control tissue obtained from the same individuals, so the possibility of individual differences may be eliminated. Since there are no reports featuring analyses using different tissue samples collected from the same patient, the present findings could offer valuable data. In addition, although tissues from the lacrimal gland were also used in this study as a control, IgG4-ROD also occurs in tissues other than the lacrimal gland [6]. The presence of lymphocyte clusters in adipose tissue has been reported [38] and may be responsible for lymphoproliferative disorders. Based on the above, adipose tissue from the same patient with IgG4-ROD remains an important control, but the possibility that this tissue is affected by the lesion should be fully addressed.
This study identified two novel genes involved in IgG4-ROD: MMP12 and SPP1. Up-regulation of these genes in IgG4-ROD relative to orbital MALT lymphoma was validated by real-time PCR, and increased expression at the protein level in the IgG4-ROD lesion compared to the control tissues was verified by immunostaining. SPP1 encodes osteopontin, which is involved in the attachment of osteoclasts to the mineralized bone matrix. MMP12 encodes an enzyme that decomposes the extracellular matrix. SPP1 was reported to be involved in fibrosis in various diseases [39][40][41][42]. The histopathology of IgG4-RD is characterized by fibrosis and lymphoid follicle formation. Until now, the TGF-β secreted by CTLs and Tregs was implicated in fibrosis [3], but it is possible that the presence of SPP1 strongly activates TGF-β-induced fibrosis. This is inferred from our finding that SPP1 is strongly expressed in the fibrotic part of IgG4-ROD tissue.
In this study, commonly up-regulated DEGs among IgG4-ROD and control tissues, MMP12 showed low p values when IgG4-ROD was compared to each control tissue. It was also reported that M2 macrophages secrete MMP12, which is an elastase, and are involved in the pathological conditions of contact dermatitis [43]. Generally, M2 macrophages have an anti-inflammatory effect [44][45][46] and are related to IgG4-RD [47]. The exhibition of an anti-inflammatory effect suggests the presence of inflammation. Some IgG4-ROD patients resist systemic steroid treatment. If IgG4-ROD is indeed an immune reaction to a pathogen, the disease would recur as the effect of steroids gradually decreases, resulting in failure to suppress the immune reaction. In the pathway analysis, pathways related to the immune system were extracted when IgG4-ROD was compared with any control tissue, strongly suggesting the possibility that immunity is involved in the pathological condition. We hypothesize that in IgG4-ROD, an infection triggers an inflammatory reaction in the lacrimal gland, and M2 macrophages secrete elastase with an anti-inflammatory effect. MMP12 is included in the pathways involved in the extracellular matrix and collagen degradation, although there was no significant difference in these pathways observed in this study. MMP12 promotes fibrosis in Schistosoma mansoni infection [48] and may also be involved in fibrosis in IgG4-ROD.
In our experiment, in which we identified the biopsy gene profile of RLH, we found that the gene profile of RLH was remarkably similar to that of IgG4-ROD and was characterized by decreased expression of MMP12 and SPP1 compared to IgG4-ROD (Table S1). Notably, we were able to differentiate these two very similar diseases at the RNA expression level. Pathologically, IgG4-ROD and RLH are distinctly different in terms of the presence or absence of fibrosis; at the gene level, IgG4-ROD differs from RLH in not only showing the up-regulation of genes related to fibrosis, such as SPP1 and MMP12, but also the attenuation of anti-fibrotic elements, such as collagen degradation, that may lead to fibrosis.
Limitations of this study include its retrospective design and relatively small number of cases collected from a single institution, which might have produced a selection bias and confounding bias. Due to the lack of comparable studies on IgG4-ROD including IgG4-RD with other diseases, it is not possible to compare the present findings with those of other studies. A prospective study with a large number of clinical samples from multiple centers is required to validate the present findings and further investigate other clinical manifestations of IgG4-ROD and other orbital tumors to diagnose clinical subtypes and elucidate their pathogeneses. In this study, we used real-time PCR to verify only two genes with enhanced expression in RNA-seq. However, it is necessary to carry out additional studies by quantifying other genes and analyzing their protein expression levels. Furthermore, MMP12 and SPP1 should be analyzed in larger samples of IgG4-ROD and other IgG4-RDs, and it is also necessary to examine the dynamics of these genes and proteins in IgG4-positive MALT lymphoma.

Conclusions
To the best of our knowledge, this is the first report of a comprehensive mRNA analysis of IgG4-ROD using RNA-seq. We succeeded in detecting not only previously reported genes but also novel genes that that may serve as biomarkers or may be involved in the pathogenesis of IgG4-ROD. In addition, these results may provide novel insights for finding new therapeutic target molecules.
Supplementary Materials: The following are available online at http://www.mdpi.com/2077-0383/9/11/3458/s1, Table S1: DEGs of RLH versus normal lacrimal glands: top 50 in descending order of p value. Table S2: Results of the pathway analysis using DEGs for IgG4-ROD tissue compared with adipose tissue. Table S3: Results of the pathway analysis using DEGs for IgG4-ROD tissue compared with MALT lymphoma tissue. Table S4: Results of the pathway analysis using DEGs for IgG4-ROD tissue compared with RLH tissue. Table S5: Results of the pathway analysis using DEGs for IgG4-ROD tissue compared with lacrimal tissue. Table S6: Results of the pathway analysis using DEGs for RLH tissue compared with IgG4-ROD.