Molecular Evidence for Precursors of Sjögren’s Foci in Histologically Normal Lacrimal Glands

Understanding the formation of Sjogren’s lymphocytic infiltrates could permit earlier diagnosis and better outcomes. We submitted gene transcript abundances in histologically normal rabbit lacrimal glands to principal component analysis. The analysis identified a cluster of transcripts associated with Sjögren’s foci, including messenger RNAs (mRNAs) for C–X–C motif chemokine ligand 13 (CXCL13) and B-cell activating factor (BAFF), which dominated the major principal component. We interpreted the transcript cluster as the signature of a cluster of integrally functioning cells. Pregnancy and dryness increased the likelihood that the cluster would develop to high levels, but responses were subject to high levels of stochasticity. Analyzing microdissected samples from high- and low-cluster-level glands, we found that certain transcripts, including mRNAs for C–C motif chemokine ligand 21 (CCL21), CXCL13, cluster of differentiation 4 (CD4), CD28, CD25, BAFF, and interleukin 18 (IL-18) were significantly more abundant in immune cell clusters (ICs) from the high-cluster-level gland; mRNAs for CCL2, CD25, and IL-1RA were significantly more abundant in acinus-duct axis samples; mRNAs for CCL4, BAFF, IL-6, and IL-10 were more abundant in some acinus-duct samples; cells with high prolactin immunoreactivity were more frequent in interacinar spaces. In conclusion, integrated functional networks comprising Sjögren’s infiltrates, such as ICs, acinar cells, ductal cells, and interacinar cells, can form in histologically normal glands, and it is feasible to detect their molecular signatures.


Introduction
Sjögren's syndrome is an autoimmune epithelitis characterized by production of antibodies and focal lymphocytic infiltrates that cause severe dysfunction of the lacrimal glands and salivary glands. The secretory dysfunction leads to ocular surface and oral pathology and symptoms. When manifestations of the autoimmune processes arise in extraglandular sites, such as the kidneys, vasa vasorum, and peripheral vessels, they further impair quality of life and significantly increase risks of life-threatening vasculitides and non-Hodgkin's mucosa-associated lymphoid tissue (MALT) lymphomas.
An understanding of how the disease processes develop and evolve could improve diagnosis and treatment. However, the etiopathogenesis of Sjögren's syndrome remains a formidable challenge.
Sjögren's syndrome comprises two major subtypes, i.e., primary Sjögren's syndrome and secondary Sjögren's syndrome [1], and both subtypes present in multiple clinical phenotypes. The cellular and molecular processes that underlie the clinical subtypes are diverse. Primary Sjögren's infiltrates may be either T-cell predominant or B-cell predominant [2], and may have either positive-or negative type I interferon (IFN) signatures [3]. B-cell predominant infiltrates may have or lack germinal centers [4]. These features imply that Sjögren's syndrome may have multiple etiologies and may develop along multiple pathogenic pathways.
Identified risk factors offer clues to Sjögren's syndrome etiology. It is known for some time that certain human leukocyte antigen (HLA) class II molecule variants favor B-cell responses to the anti-Sjögren's-syndrome-related A (Ro/SSA) and anti-Sjögren's-syndrome-related B (La/SSB) autoantigens, formation of lymphocytic infiltrates with germinal centers, and risks for MALT lymphomas and vasculitides, evidently because the high-risk HLA class II variants have high affinities for Ro/SSA and La/SSB epitopes [5][6][7]. Additional genetic polymorphisms reported to be associated with risk for Sjögren's syndrome include variants of interleukin 10 (IL-10) [8][9][10], tumor necrosis factor (TNF) [11], antigen peptide transporter 2 (TAP2) [12], and 2 -5 -oligoadenylate synthetase 1 (OAS1) [13]. Female sex is a well-known risk factor [14], and parity was also reported to be a risk factor [15,16]. The concordance rate for monozygotic twins was estimated to be low [17]. Therefore, interactions between genetic predisposition, sex, and parity are not sufficient to determine Sjögren's etiology. Other phenomena, assumed to be environmental in nature, must be involved.
Many of the immune response-related molecules expressed by epithelial cells of salivary glands with Sjögren's lesions also are expressed by epithelial cells in salivary glands from healthy control subjects, but generally at much lower levels than in affected glands. Therefore, the data suggest that exposures to certain risk factors increase the probability that exocrine gland epithelial cells will upregulate their expression of genes that cause them to either (1) function as surrogate antigen presenting cells, or (2) create local microenvironments where interactions between professional antigen-presenting cells, autoreactive T cells, and autoreactive B cells are likely lead to immune cell activation and self-organization of ectopic lymphoid structures.
Working with histologically normal lacrimal glands from healthy rabbits, we confirmed Frey and coworkers' finding that lacrimal gland epithelial cells express exhibit immunohistochemical positivity for prolactin (PRL) [40], which has pleiotropic actions as a cytokine [41][42][43][44]. We found that corneal injury, corneal adenovirus infection, and pregnancy altered levels of PRL immunoreactivity [45,46]. We also found that lacrimal glands express messenger RNA (mRNA) for PRL and also mRNAs for class II major histocompatibility complex (MHC II), CD80, and CD86, numerous cytokines, and numerous chemokines. The abundances of many transcripts exhibited considerable gland-to-gland variability. Some of the abundance variations appeared to be related to the dryness, while others were related to the temperatures and to the environment the animals experienced prior to study [47,48]; additional variability was induced by the hormonal environment of pregnancy [49]. However, much of the variability appeared related to stochastic phenomena localized within glands. Pearson's analysis of the variability showed that variations of clusters of transcripts were significantly correlated. Principal component analyses confirmed the major correlation cluster and indicated that certain transcripts might be expressed by two or more correlation clusters in the same gland. We interpreted the transcript correlation clusters as signatures of clusters of cells that were interacting with each other coordinately.
The major correlation cluster was of interest, as it included mRNA for BAFF, a B-cell mitogen and activating factor, and mRNAs for the B-cell chemokine, CXCL13, and the T-cell cytokine, CCL21, all associated with Sjögren's lymphocytic foci. A preliminary laser capture microdissection of glands from a nulliparous animal confirmed that epithelial cells expressed a number of immune response-related gene transcripts, that immune cells in small clusters expressed certain other transcripts, and that both epithelial cells and immune cells expressed certain transcripts. Notably, the immune cell clusters, or "accumulations", were too small to be identified as Sjögren's foci. These characteristics suggest that epithelial cells and immune cells can engage in multiple signaling interactions with each other, forming localized networks that might expand over time to eventually manifest as histopathologically identifiable Sjögren's lesions.
We collated and analyzed data from our previously published studies to discern how age, pregnancy, and exposures varying ambient conditions influence the formation of correlation clusters. We also analyzed the abundances of selected transcripts in samples of acinar cells, intralobular duct cells, interlobular duct cells, intralobar duct cells, and clustered immune cells obtained by laser capture microdissection of three glands from term-pregnant animals. Two glands were representative of the subgroup of glands in which the Sjögren's-resembling transcript correlation cluster was present at low levels, and one gland was representative of the subgroup in which the Sjögren's-resembling transcript correlation cluster was present at high levels.
We found that epithelial cells rarely express mRNA for MHC II, but frequently express mRNA for CD1d, which presents both bacterial and autologous glycolipids invariant α-chain natural killer (NK) T cells. They even more frequently express mRNA for MHC I, which presents epitopes of intracellular proteins, both viral and autologous, to CD8 + T cells. We also found that samples of the epithelial segments heterogeneously express costimulatory molecules, chemokines, and cytokines. These findings support the hypothesis that epithelial cells in histologically normal glands can interact with clusters of immune cells to comprise networks with stochastically varying cellular compositions and transcript expression profiles. Environmental and hormonal exposures increase the likelihoods that small networks will expand and, perhaps, evolve as ectopic lymphoid structures of different phenotypes. The ability to detect molecular signatures of the small networks in histologically normal glands may have important implications for the future diagnosis and proactive treatment of incipient Sjögren's syndrome.

Systemic and Stochastic Variations in Immune Response-Related Gene Transcript Expression
To compare responses to exposures to varying combinations of ambient dryness and temperature, the hormonal environment of pregnancy, and the transition from young adulthood to mature adulthood, we collated abundance data for transcripts that were measured in all three previous studies and submitted the data to principal component analysis (PCA). Groups V.G2, V.G3, V.G5, P.G5, and V.G6 were 20 weeks old and were raised and housed in a barrier-free facility, where they experienced different environmental exposures. Group V.G7 was 52 weeks old and was raised and housed in barrier facilities. Group P.G5 was term-pregnant; the other groups were nulliparous. Figure 1 presents the average daily high temperatures and average maximum degree of dryness that the six groups experienced. from group P.G5 (i.e., the term-pregnant animals) exhibited the most dispersed distributions with respect to PC1; one clearly-defined subgroup of glands exhibited positive PC1 projections within the empirical normal range; the other subgroup exhibited negative PC1 projections that extended well beyond the ranges of the glands from any of the other groups. We designated the respective subgroups P.G5.B and P.G5.A, and we selected one gland from each subgroup for further analysis.  The individual glands' projections with respect to the significant principal components are presented in Figure 2. The glands from the two groups that experienced the most benign ambient conditions (i.e., V.G7 and V.G3) clustered closely together within narrow ranges of PC1, PC2, and PC3 projections, indicated by brackets in Figure 2A,B. We interpreted these as empirical normal ranges. The glands from groups V.G2, V.G5, and V.G6 (i.e., nulliparous groups that experienced hotter, drier, or hotter and dryer conditions) were distributed more dispersedly with respect to PC1. Projections of some V.G2 and V.G5 glands fell within the empirical normal range; projections of other V.G2 and V.G5 glands, and of all V.G6 glands, fell well below the empirical normal range. The glands from group P.G5 (i.e., the term-pregnant animals) exhibited the most dispersed distributions with respect to PC1; one clearly-defined subgroup of glands exhibited positive PC1 projections within the empirical normal range; the other subgroup exhibited negative PC1 projections that extended well beyond the ranges of the glands from any of the other groups. We designated the respective subgroups P.G5.B and P.G5.A, and we selected one gland from each subgroup for further analysis.
The V.G3, V.G5, V.G6, and P.G5.A glands clustered together with respect to their PC4 and PC5 projections, again defining empirical normal ranges ( Figure 2C). The V.G7 glands fell outside the empirical normal ranges, possibly reflecting an influence of the animals' transition from young adulthood to maturity. Four of the P.G5 glands fell outside or well outside the empirical normal ranges, one in each quadrant, possibly reflecting stochastic responses to the hormonal environment of pregnancy.
Median abundances of numerous transcripts in the glands from groups V.G2, V.G3, V.G5, and V.G6 could be described by heuristics related to the average daily high temperature or the average daily maximum degree of dryness the animals experienced [47,48]. Therefore, we analyzed the relationships between the environmental variables and the principal component projections of the glands from all five groups of nulliparous animals described in this study. As shown in Figure 3A, the median PC2 projections could be modeled as decreasing exponentially with exposure to increasing daily high temperatures (R 2 = 0.981). As shown in Figure 3B, the median PC1 projections could be modeled as decreasing exponentially with exposure to increasing degrees of dryness (R 2 = 0.960). The median PC1 projection of group V.G6 glands was notably displaced from the exponential growth model prediction. The V.G6 glands showed a significant exponential relationship between decreasing PC1 projections and increasing PC3 projections ( Figure 2B, R 2 = 0.891), suggesting that a phenomenon related to the PC3 projections displaced their PC1 projections above the value predicted in Figure 3B. A total of 72 transcripts were assayed in whole-gland samples from groups V.G5 and P.G5 [49]. Discrepancies between principal component values from that dataset and the collated dataset resulted primarily from the data from groups V.G2 and V.G6. Arrows indicate PC1 projections of glands that were selected for microdissection.
The V.G3, V.G5, V.G6, and P.G5.A glands clustered together with respect to their PC4 and PC5 projections, again defining empirical normal ranges ( Figure 2C). The V.G7 glands fell outside the empirical normal ranges, possibly reflecting an influence of the animals' transition from young A total of 72 transcripts were assayed in whole-gland samples from groups V.G5 and P.G5 [49]. Discrepancies between principal component values from that dataset and the collated dataset resulted primarily from the data from groups V.G2 and V.G6. Arrows indicate PC1 projections of glands that were selected for microdissection.
could be modeled as decreasing exponentially with exposure to increasing degrees of dryness (R = 0.960). The median PC1 projection of group V.G6 glands was notably displaced from the exponential growth model prediction. The V.G6 glands showed a significant exponential relationship between decreasing PC1 projections and increasing PC3 projections ( Figure 2B, R 2 = 0.891), suggesting that a phenomenon related to the PC3 projections displaced their PC1 projections above the value predicted in Figure 3B.  Figure 3A).
The abundances of numerous transcripts exhibited low concordances between right eye (oculus dextrus (OD))-associated and left eye (oculus sinister (OS))-associated lacrimal glands from the group V.G5 animals [47]. To assess the relative contributions that systemic factors and strictly local stochastic factors made to PC1 projections of the P.G5.B and P.G5.A glands, we plotted PC1 projections of the companion right eye OD-associated and left eye OS-associated glands from each group P.G5 animal ( Figure 4). Glands P.G5.05.OD and P.G5.05.OS were the only companions that exhibited similar PC1 projections. The poor concordances between companion glands indicate that strictly local stochastic factors contribute substantially to PC1 projection variations. Table 1 presents the transcript loadings identified by principal analysis of the collated data. Messenger RNAs for IL-1⍺, IL-1β, IL-6, IL-10, CCL2, CCL4, CCR5, CXCL13, CD4, CD8, CD28, cytotoxic T-lymphocyte-associated antigen 4 (CTLA-4), BAFF, MHC II, and PRL contributed negative loadings to PC1 in this analysis. Many of these transcripts correlated strongly with each other when we submitted the complete datasets for the V.G2, V.G3, V.G5, and V.G6 glands to separate Pearson's tests [47]. They also tended to contribute strong negative loadings to the first principal component when we submitted the complete datasets for the V.G7 glands [50] and the V.G5 and P.G5 glands to separate to principal component analyses [49]. Exceptions to this generalization resulted from having collated data from groups V.G2 and V.G6 together with data from the other groups.  Figure 3A).
The abundances of numerous transcripts exhibited low concordances between right eye (oculus dextrus (OD))-associated and left eye (oculus sinister (OS))-associated lacrimal glands from the group V.G5 animals [47]. To assess the relative contributions that systemic factors and strictly local stochastic factors made to PC1 projections of the P.G5.B and P.G5.A glands, we plotted PC1 projections of the companion right eye OD-associated and left eye OS-associated glands from each group P.G5 animal ( Figure 4). Glands P.G5.05.OD and P.G5.05.OS were the only companions that exhibited similar PC1 projections. The poor concordances between companion glands indicate that strictly local stochastic factors contribute substantially to PC1 projection variations.   LG.OD LG.OS  Table 1 presents the transcript loadings identified by principal analysis of the collated data. Messenger RNAs for IL-1α, IL-1β, IL-6, IL-10, CCL2, CCL4, CCR5, CXCL13, CD4, CD8, CD28, cytotoxic T-lymphocyte-associated antigen 4 (CTLA-4), BAFF, MHC II, and PRL contributed negative loadings to PC1 in this analysis. Many of these transcripts correlated strongly with each other when we submitted the complete datasets for the V.G2, V.G3, V.G5, and V.G6 glands to separate Pearson's tests [47]. They also tended to contribute strong negative loadings to the first principal component when we submitted the complete datasets for the V.G7 glands [50] and the V.G5 and P.G5 glands to separate to principal component analyses [49]. Exceptions to this generalization resulted from having collated data from groups V.G2 and V.G6 together with data from the other groups.

Distribution of Transcript Expression among Epithelial Segments and Immune Cells
We used laser capture methodology to microdissect gland P.G5.06.OS, a representative subgroup P.G5A gland, and glands P.G5.01.OD and P.G5.03.OS, representative of subgroup P.G5.B. The selected glands are indicated by arrows in Figure 2A. We assayed abundances of selected transcripts in all samples from gland P.G5.06.OS and in immune cell clusters from gland P.G5.01.OD; we also assayed abundances of a subset of the selected transcripts in samples of acinar cells, intralobular duct cells, interlobular duct cells, and intralobar duct cells from gland P.G5.03.OS. The transcripts we assayed in the microdissected samples were a subset of the transcripts measured in whole-gland samples of the V.G5 glands, P.G5.A glands, and P.G5.B glands [49]. As this subset was larger than the subset analyzed in Figures 2 and 3, we interpreted the findings with reference to the large sets of principal component projections and loadings identified by principal component analysis of the complete set of transcripts in samples of the V.G5 glands, P.G5.A glands, and P.G5.B glands [49], rather than the projections in Figures 2 and 3 and the loadings in Table 1.
Many transcripts were present at detectable levels in acinar cell and duct cell samples, as well as in immune cell cluster samples. To estimate approximate relative amounts in acini, the duct segments, and immune cells, we used a model that weighed abundances in acinar cell, intralobular duct cell, interlobular duct cell, intralobar duct cell, and immune cell cluster samples by factors of 0.80, 0.04, 0.04, 0.04, and 0.08. The modeled amounts of mRNAs for CCL2 and CD25 were higher for gland P.G5.06.OS than the modeled amounts of P.G5.B glands ( Figure 5). These relationships recapitulated the relationships between the transcripts' abundances in the whole-gland samples. Therefore, the differences between abundances in the P.G5.A and P.G5.B whole-gland samples can be attributed to significantly higher abundance of mRNA for CCL2 in acinar cells from gland P.G5.06.OS, and to significantly higher abundances of mRNA for CD25 in both immune cells and acinar cells from gland P.G5.06.OS. These findings indicate that acinar cells contributed mRNA for CCL2, and both acinar cells and immune cells contributed mRNA for CD25 to the negative PC1-loading transcript cluster. Additional transcripts were also present at detectable levels in epithelial segment and immune cell samples (Figure 6). Messenger RNAs for a proliferation-inducing ligand (APRIL), polymeric immunoglobulin receptor (pIgR), epidermal growth factor (EGF), and lipophilin CL were predominantly expressed in acinar cells; none of these transcripts contributed to the PC1-negative loading transcript cluster [49]. Messenger RNAs for CD3ε, CCL21, CD4, CD3ζ, CXCL13, and MHC II were predominantly expressed in immune cell clusters, as discussed further below.  Figure 2) and a representative subgroup P.G5.A gland (negative PC1 projection in Figure 2). Error bars indicate pooled standard errors for IC and the four epithelial elements. The transcripts selected showed statistically significant differences between epithelial segments from the P.G5.A gland and the P.G5.B glands (significant p-values and p-values indicating trends (t) are as indicated). (B) Measured abundances of the transcripts in whole-gland samples of the P.G5.B and P.G5.A subgroup glands. Glands that were microdissected are indicated by large-font symbols.
The modeled amount of mRNA for TGF-β 2 was lower for gland P.G5.06.OS than for the P.G5.B glands, again recapitulating the relationship between the transcript's abundances in the whole-gland samples. The difference can be attributed to significantly lower abundances in intralobular duct cells, interlobular duct cells, intralobar ducts cells, and immune cells from gland P.G5.06.OS.
In contrast to the abundances of mRNAs for CCL2, CD25, and TGF-β 2 , the modeled abundances of mRNA for TGF-β 1 in gland P.G5.06.OS were significantly higher than the modeled abundances in the P.G5.B glands, contrary to the relationships between the abundances measured in the whole-gland samples ( Figure 5). It is possible that a cell type that we did not sample by laser capture microdissection of the P.G5.B gland expressed high levels of the mRNA for TGF-β 1 .
Additional transcripts were also present at detectable levels in epithelial segment and immune cell samples ( Figure 6). Messenger RNAs for a proliferation-inducing ligand (APRIL), polymeric immunoglobulin receptor (pIgR), epidermal growth factor (EGF), and lipophilin CL were predominantly expressed in acinar cells; none of these transcripts contributed to the PC1-negative loading transcript cluster [49]. Messenger RNAs for CD3ε, CCL21, CD4, CD3ζ, CXCL13, and MHC II were predominantly expressed in immune cell clusters, as discussed further below.  Figure 2); values are normalized to total mean modeled abundances. Error bars indicate pooled normalized standard errors for IC and the four epithelial elements. Transcripts depicted in Figure 5 are not included.
Other transcripts were expressed at appreciable levels in both epithelial segments and in immune cell clusters. In general, the modeled abundances of the transcripts that were expressed by both epithelial cells and immune cells exhibited much larger variances than the transcripts that were predominantly expressed by either acinar cells and or immune cells. Consequently, apparent differences between modeled abundances in the epithelial segment samples were often not statistically significant. Therefore, the findings indicate that acinar cells, intralobular duct cells, interlobular duct cells, and intralobar duct cells expressed detectible levels of the immune responserelated gene transcripts, but they tended to be extremely heterogeneous with respect to the levels at which they expressed the transcripts. Messenger RNA for CCL4 appeared to be more abundant in acinar cells from gland P.G5.06.OS than in acinar cells from the P.G5.B gland. Similarly, mRNAs for BAFF, IL-6, and IL-10 appeared to be more abundant in both acinar cells and duct segment cells from gland P.G5.06.OS, but they exhibited large variances, and differences from the corresponding epithelial segments from the P.G5.B gland were not statistically significant. Therefore, it appears that some acini and duct segments contributed mRNAs for CCL4, BAFF, IL-6, and IL-10 to the PC1negative loading transcript cluster. Acini and duct segments also expressed mRNAs for CCL28 and may have contributed it to the PC1-negative loading transcript cluster.

Transcript Expression in Immune Cell Clusters
Immune cell cluster samples from gland P.G5.06.OS differed significantly from immune cell cluster samples from the P.G5.B gland with respect to the levels at which they expressed 20 of the assayed transcripts (Figure 7). Messenger RNAs for CCL21, CXCL13, CD3ε, CD3ζ, CD25, CD80, IL-1RA, IL-18, pIgR, matrix metalloproteinase 9 (MMP-9), CD4, CD28, CTLA-4, BAFF, and TGF-β1 were significantly more abundant in immune cell clusters from gland P.G5.06.OS. Messenger RNAs for TGF-β2, MHC I, APRIL, and EGF were significantly less abundant in immune cell clusters from gland P.G5.06.OS, and there was a trend (p = 0.06) toward a lower abundance of mRNA for CCL4 in immune cells from gland P.G5.06.OS. Notably, the transcripts that were less abundant in immune cell clusters from gland P.G5.06.OS were predominantly expressed by epithelial cells, rather than by immune cell clusters ( Figure 5,6).
Of the transcripts that were more highly expressed in immune cells from gland P.G5.06.OS, mRNAs for CCL21, CD4, CD25, CD25, CTLA-4, BAFF, and MMP-9 contributed strong negative loadings to PC1 [49]. These findings indicate that immune cell clusters contributed mRNAs for CCL21, CD4, CD25, CD25, CTLA-4, BAFF, and MMP-9 to the PC1-negative loading transcript cluster and, therefore, participated in a functional network with acinar cells contributing mRNA for CCL2, additional mRNA for CD25 and, presumably, mRNAs for CCL4, BAFF, IL-6, and IL-10.  Figure 2); values are normalized to total mean modeled abundances. Error bars indicate pooled normalized standard errors for IC and the four epithelial elements. Transcripts depicted in Figure 5 are not included.
Other transcripts were expressed at appreciable levels in both epithelial segments and in immune cell clusters. In general, the modeled abundances of the transcripts that were expressed by both epithelial cells and immune cells exhibited much larger variances than the transcripts that were predominantly expressed by either acinar cells and or immune cells. Consequently, apparent differences between modeled abundances in the epithelial segment samples were often not statistically significant. Therefore, the findings indicate that acinar cells, intralobular duct cells, interlobular duct cells, and intralobar duct cells expressed detectible levels of the immune response-related gene transcripts, but they tended to be extremely heterogeneous with respect to the levels at which they expressed the transcripts. Messenger RNA for CCL4 appeared to be more abundant in acinar cells from gland P.G5.06.OS than in acinar cells from the P.G5.B gland. Similarly, mRNAs for BAFF, IL-6, and IL-10 appeared to be more abundant in both acinar cells and duct segment cells from gland P.G5.06.OS, but they exhibited large variances, and differences from the corresponding epithelial segments from the P.G5.B gland were not statistically significant. Therefore, it appears that some acini and duct segments contributed mRNAs for CCL4, BAFF, IL-6, and IL-10 to the PC1-negative loading transcript cluster. Acini and duct segments also expressed mRNAs for CCL28 and may have contributed it to the PC1-negative loading transcript cluster.

Transcript Expression in Immune Cell Clusters
Immune cell cluster samples from gland P.G5.06.OS differed significantly from immune cell cluster samples from the P.G5.B gland with respect to the levels at which they expressed 20 of the assayed transcripts ( Figure 7). Messenger RNAs for CCL21, CXCL13, CD3ε, CD3ζ, CD25, CD80, IL-1RA, IL-18, pIgR, matrix metalloproteinase 9 (MMP-9), CD4, CD28, CTLA-4, BAFF, and TGF-β 1 were significantly more abundant in immune cell clusters from gland P.G5.06.OS. Messenger RNAs for TGF-β 2 , MHC I, APRIL, and EGF were significantly less abundant in immune cell clusters from gland P.G5.06.OS, and there was a trend (p = 0.06) toward a lower abundance of mRNA for CCL4 in immune cells from gland P.G5.06.OS. Notably, the transcripts that were less abundant in immune cell clusters from gland P.G5.06.OS were predominantly expressed by epithelial cells, rather than by immune cell clusters (Figures 5 and 6).
Of the transcripts that were more highly expressed in immune cells from gland P.G5.06.OS, mRNAs for CCL21, CD4, CD25, CD25, CTLA-4, BAFF, and MMP-9 contributed strong negative loadings to PC1 [49]. These findings indicate that immune cell clusters contributed mRNAs for CCL21, CD4, CD25, CD25, CTLA-4, BAFF, and MMP-9 to the PC1-negative loading transcript cluster and, therefore, participated in a functional network with acinar cells contributing mRNA for CCL2, additional mRNA for CD25 and, presumably, mRNAs for CCL4, BAFF, IL-6, and IL-10. The findings in Figure 5,6,7 indicate that immune cell clusters in gland P.G5.06.OS contained an admixture of cells from other transcript clusters in addition to the PC1-negative loading cluster. Transcripts that contributed relatively weak negative loadings to PC1, i.e., mRNAs for CD3ζ, CD80, and CXCL13 [49], were also significantly more abundant in immune cell clusters in gland P.G5.06.OS. Notably, mRNA for CD3ζ contributed a strong positive loading, and mRNAs for CD80 and CXCL13 contributed positive, albeit weaker, loadings to PC3 [49]. In addition to its large negative projection with respect to PC1, gland P.G5.06.OS also had a large positive projection with respect to PC3, indicating that cells of the PC3-positive loading cluster were present in that gland. In contrast, gland P.G5.01.OS had a negative projection with respect to PC3. Moreover, principal component analysis of transcript abundances in the microdissected immune cell cluster samples only ( Figure 8 and Table  2) indicated that mRNAs for CD3ζ, CD80, and CXCL13 contributed strong loadings to the major PC, of the same sign as the loadings contributed by mRNAs for CCL21, CD4, CD25, CD25, CTLA-4, BAFF, and MMP-9. The findings in Figures 5-7 indicate that immune cell clusters in gland P.G5.06.OS contained an admixture of cells from other transcript clusters in addition to the PC1-negative loading cluster. Transcripts that contributed relatively weak negative loadings to PC1, i.e., mRNAs for CD3ζ, CD80, and CXCL13 [49], were also significantly more abundant in immune cell clusters in gland P.G5.06.OS. Notably, mRNA for CD3ζ contributed a strong positive loading, and mRNAs for CD80 and CXCL13 contributed positive, albeit weaker, loadings to PC3 [49]. In addition to its large negative projection with respect to PC1, gland P.G5.06.OS also had a large positive projection with respect to PC3, indicating that cells of the PC3-positive loading cluster were present in that gland. In contrast, gland P.G5.01.OS had a negative projection with respect to PC3. Moreover, principal component analysis of transcript abundances in the microdissected immune cell cluster samples only ( Figure 8 and Table 2) indicated that mRNAs for CD3ζ, CD80, and CXCL13 contributed strong loadings to the major PC, of the same sign as the loadings contributed by mRNAs for CCL21, CD4, CD25, CD25, CTLA-4, BAFF, and MMP-9.  Figure 2) and (B) the subgroup P.G5.A gland (negative PC1 projection in Figure 2).

Cells Outside the Acinus Duct/Immune Cell Cluster Axis
The modeled abundances of mRNA for PRL were remarkably disparate from the measured abundances of mRNA for PRL in the whole-gland samples (Figure 9). The modeled abundance of mRNA PRL for gland P.G5.06.OS samples was much lower than the modeled abundances for  Figure 2) and (B) the subgroup P.G5.A gland (negative PC1 projection in Figure 2).

Cells Outside the Acinus Duct/Immune Cell Cluster Axis
The modeled abundances of mRNA for PRL were remarkably disparate from the measured abundances of mRNA for PRL in the whole-gland samples (Figure 9). The modeled abundance of mRNA PRL for gland P.G5.06.OS samples was much lower than the modeled abundances for samples that were microdissected from a V.G5 gland [47]. This finding is consistent with the hypothesis that epithelial cells in gland P.G5.06.OS downregulated their expression of PRL in response to the high levels of PRL that the anterior pituitary contributes to the circulation during pregnancy. The large variances of the modeled abundances in microdissected samples from the P.G5.B gland precludes a comparison. However, the measured abundance of mRNA for PRL in gland P.G5.06.OS sample was much higher than in the P.G5.B whole-gland samples and the V.G5 whole-gland samples. This finding suggests that high levels of PRL were expressed by cells that we did not sample by microdissection of gland P.G5.06.OS. To test this hypothesis, we stained frozen sections of three P.G5.B glands and three P.G5.A glands for PRL immunopositivity. Cells, possibly macrophages [51], that were intensely positive for PRL immunoreactivity were scattered in interacinar spaces, where they would not have been sampled by our protocol, and they were 8.15-fold more frequent (p < 0.001) in gland P.G5.06.OSs than in the P.G5.B glands (Figure 9). The finding of a high frequency of PRL High cells in P.G5.A glands is consistent with the hypothesis that cells outside the acinus duct axis participated in the Sjögren's-like epithelium immune cell network.

Discussion
Our findings indicate that clusters of immune cells expressing high levels of many of the same transcripts that are highly expressed in Sjögren's infiltrates are present in histologically normal lacrimal glands. As the clusters can be detected before classical Sjögren's infiltrates are identified, we suggest that they may be early precursors of Sjögren's infiltrates. This conjecture is consistent with the finding that levels of the typical Sjögren's autoantibodies can be elevated before Sjögren's

Discussion
Our findings indicate that clusters of immune cells expressing high levels of many of the same transcripts that are highly expressed in Sjögren's infiltrates are present in histologically normal lacrimal glands. As the clusters can be detected before classical Sjögren's infiltrates are identified, we suggest that they may be early precursors of Sjögren's infiltrates. This conjecture is consistent with the finding that levels of the typical Sjögren's autoantibodies can be elevated before Sjögren's syndrome symptoms become clinically significant [52].
Our findings are consistent with the hypothesis that signaling interactions with epithelial cells in all levels of the acinus duct axis, as well as with cells outside the acinus duct axis, influence formation of the immune cell clusters that evade peripheral tolerance mechanisms. Acinar cells maintain high baseline levels of APRIL. Acinar and ductal cells appear to contribute CCL2, CCL4, CCL28, BAFF, IL-6, and IL-10. Immune cell clusters appear to contribute CXCL13, CCL21, MHC II, and MMP-9 [53], as well as additional BAFF. The levels at which mRNAs for these proteins are expressed vary coordinately, and through a wide range, as indicated by the range of gland projections with respect to PC1 (Figure 2).
Because mRNA for MHC II was predominantly expressed in immune cell clusters, it appears that induction of epithelial HLA class II molecule expression may be a downstream event in Sjögren's pathogenesis, rather than a triggering event in Sjögren's etiology. Notably, however, acinar epithelial cells expressed mRNAs for MHC I, CD1d, CD80, and CD86. Therefore, lacrimal gland epithelial cells might contribute to Sjögren's etiology not only by expressing chemokines that recruit T cells, B cells, and professional antigen-presenting cells, but also by expressing mitogenic factors that support immune cell survival and CD4 + cell proliferation in response to MHC II-restricted epitopes presented by the professional antigen-presenting cells. They might also by displaying costimulatory signals, presenting MHC I-bound autoantigen epitopes to CD8 + T cells, and presenting CD1d-bound glycolipids to invariant α-chain NK T cells. Presumably, the antigen receptors of CD4 + or CD8 + T cells that participate in such interactions would have autoantigen epitope affinities strong enough to support positive selection during thymic maturation, but too weak to determine activation-induced deletion.
A corollary to this hypothesis is that, in addition to the coreceptors, CD80 and CD86, epithelial cell also express additional coreceptors and receptors, typically expressed by T cells, which mediate immune cell-to-epithelial cell paracrine signaling, as well as epithelial cell autocrine signaling and epithelial cell-to-epithelial cell paracrine signaling. As indicated in Figure 5, epithelial cells expressed mRNA for CD25, the α-subunit of the IL-2 receptor; mRNA for the CCL4 receptor, CCR5; mRNAs for CD8 and CD28 ( Figure 6); and mRNA for IL-2 [47]. These transcripts all contributed to the PC1-negative loading transcript cluster. These findings are not unprecedented, as CD8 is known to be expressed by NK cells, dendritic cells, and cortical thymocytes, as well as by T cells; CD28 is known to be expressed by plasmacytoid dendritic cells, plasmacytes, NK cells, eosinophils, and neutrophils, as well as by T cells [54]; and IL-2 receptors are known to be expressed by renal tubular epithelial cells [55,56]. Their expression by epithelial cells in the lacrimal gland might contribute to the positive feedback loop which coordinates epithelial cell expression and immune cell expression of the transcripts of the PC1-negative loading cluster.
The resemblance between the transcript expression profiles of immune cell clusters and Sjögren's foci is likely even greater than revealed by the present microdissection data, as PC1 also receives strong negative loadings from additional transcripts that also are expressed in Sjögren's foci. The transcripts that contribute strong negative loadings to PC1 include CD40L and CD40 [49], a cognate coreceptor pair expressed respectively by T cells and by B cells and professional antigen-presenting cells. They also include mRNAs for CD19 and CD72 [49], expressed by resting B cells; mRNAs for IL-2, IL-4, IL-5, IL-7, IL-13, IL-17A [49], which contribute to B-cell activation; and mRNA for IL-21 [49], which supports T-cell survival in T-cell zones and dark zone bases.
Several transcripts that are particularly associated with active lymphoid follicles notably do not contribute strong negative loadings to PC1. These include mRNA for lymphotoxin beta (LT-β), expressed by lymphoid tissue-organizing cells; mRNA for CXCL12, which recruits B cells to germinal center dark zones; mRNA for amyloid precursor protein intracellular domain (AICD), which mediates immunoglobulin gene somatic hyper mutation; mRNA for CD22, expressed by mature B cells; and mRNA for CD138, expressed by terminally differentiated plasmacytes [49]. Each of these transcripts contributes a strong or moderately strong negative loading to PC2. Gland P.G5.06.OS had a positive projection with respect to PC2, but gland P.G5.03.OD, another subgroup P.G5.A gland that had an even larger negative projection with respect to PC1, had a large negative projection with respect to PC2. Our finding (Figure 7) that immune cell cluster samples from gland P.G5.06.OS had significantly higher abundances of transcripts that contributed strong loadings to other PCs is informative in this context, as it indicates that transcripts which contribute to other PCs can be expressed in the same immune cell clusters as the PC1-negative loading transcripts. In other words, immune cell clusters in gland P.G5.03.OD may have been relatively further along on a trajectory toward manifestation as recognizable ectopic lymphoid structures with germinal centers. Moreover, because autoreactive B cells escape self-recognition checkpoints when they are activated in ectopic lymphoid structures [57][58][59], one might ask whether B cells in the immune cell clusters we described contribute autoantibodies to the circulation even before germinal centers are formed.
The remarkably broad distribution of lacrimal glands' projections with respect to PC1 as depicted in Figure 2 may offer a hint for answering the question of why the prevalence of primary Sjögren's syndrome is relatively low (0.02-0.1%) [60] That is, immune cell clusters are only likely to expand into Sjögren's foci in the glands in which they are already well developed, i.e., in glands with the largest negative projections with respect to PC1. In view of the broad distributions of lacrimal glands' projections with respect to PC2 and PC3, and the broad distributions of subgroup P.G5.A glands' projections with respect to PC4 and PC5 as depicted in Figure 2, the conclusion that immune cell clusters express transcripts from multiple transcript clusters may have general implications for understanding the heterogeneities of both primary and secondary Sjögren's syndrome phenotypes.
In rabbits, the likelihood that a gland will achieve a large negative PC1 projection is increased by exposure to high degrees of environmental dryness, as well as by the hormonal environment of pregnancy ( Figure 3). Environmental dryness is a risk factor for dry-eye disease in humans [61]; to our knowledge, it is yet to be reported as a risk factor for Sjögren's syndrome. However, the findings in Figures 2-4 demonstrate that responses to both environmental dryness and the hormonal environment of pregnancy are subject to very high levels of stochasticity.
In addition to the genetic risk factors we discussed in Section 1 of this paper, mechanisms related to viral infections were proposed [62][63][64][65][66]; these include mimicry of autoantigens by viral proteins, apoptosis-associated release of autoantigens, apoptosis-associated release of ligands for molecular pattern recognition receptors, and expression of proliferation factors and chemokines in latently infected T cells or B cells [67,68]. Occupational exposure to organic solvents [69] and exposure to psychosocial stress [70] were also proposed to be risk factors. Our findings indicate that epithelial cells in histologically normal lacrimal glands are able to recruit immune cells, provide mitogenic support for them, and shape their activities. These abilities may contribute to the responses to such risk factors, and the inherent stochasticity of the development of immune cell/epithelial cell networks would contribute to the stochasticity of responses to the risk factors.

Animals
We analyzed lacrimal glands from six groups of female rabbits, all described previously [47][48][49][50]. Groups V.G2, V.G3, V.G5, and V.G6 were 20-week-old nulliparous adults. Incomplete data were obtained from a sixth group of 20-week-old rabbits (V.G4) and were not included in the analysis. Group P.G5 was term-pregnant 20-week-old rabbits. Group V.G7 was 52-week-old nulliparous adults [50]. Groups V.G2, V.G3, V.G5, P.G5, and V.G6 were obtained from a barrier-free facility (Irish Farms, Norco, CA, USA). They were raised at different times and, thus, were exposed to different environmental temperatures and different degrees of environmental dryness prior to arriving at the University of Southern California Health Sciences Campus Vivaria (USC HSC, Los Angeles, CA, USA) for a four-day acclimation period. Group V.G7 was raised and maintained in barrier facilities (Covance Research Products, Denver, PA, and the USC HSC Vivaria), where the animals were exposed to consistently mild temperatures and mild degrees of dryness. All procedures with animals conformed with the Association for Research in Vision and Ophthalmology (ARVO) statement on the Use of Animals in Vision Research. All procedures with animals were approved by the University of Southern California Institutional Animal Care and Use Committee, protocol #20057, 24 August 2016.

Tissue Collection
We euthanized rabbits with Euthasol ® after sedating with ketamine/xylazine. We removed lacrimal glands at necropsy in RNAse-free conditions and divided each gland into three parts. We placed one part in RNALater ® for RNA extraction, one part in formalin for paraffin embedding and hematoxylin and eosin (H&E) and trichrome staining studies (not reported), and one part in OCT ® for rapid freezing in liquid nitrogen, sectioning, immunohistochemical staining, and laser capture microdissection.

Laser Capture Microdissection
We collected frozen sections on membrane-coated slides (PEN ® ; Leica Microsystems, Deerfield, IL, USA) and stained them with cresyl violet from Applied Biosystems (Foster City, CA, USA) [71]. We then microdissected samples of epithelial cells from acini; samples of epithelial cells from intralobular ducts, interlobular ducts, and intralobar ducts; and samples of immune cells from clusters around interlobular ducts using the PixCell II LCM System ® (Arcturus Bioscience, Mountain View, CA, USA) and Cap-Sure HS ® laser capture microdissection caps from Molecular Devices, Sunnyvale, CA, USA. We collected five or six replicate samples of the epithelial segments and five samples of immune cell clusters from each gland; each sample contained approximately 100 cells.

Real-Time Reverse-Transcriptase Polymerase Chain Reaction (Real Time RT-PCR)
We extracted mRNA with RNAqueous Micro ® kits from Ambion (Austin, TX, USA) according to the manufacturer's protocol, eluting with 11-L aliquots of prewarmed elution solutions. After treatment with DNAase and quality control with an ND-1000 spectrophotometer from Nanodrop Technologies (Wilmington, DE, USA), we reversed-transcribed to complementary DNA (cDNA) with High-Capacity cDNA Reverse Transcription Kits and RNase Inhibitor from Applied Biosystems using a DNA Engine ® thermal cycler from Bio-Rad (Hercules, CA, USA). We preamplified cDNA using TaqMan PreAmp Master Mix ® from Applied Biosystems. We performed real-time RT-PCR with a Prism 790HT ® from Applied Biosystems using primer and probe sequences selected with the Primer Express ® program and synthesized by Applied Biosystems. We expressed all transcript abundances relative to the abundance of mRNA for glyceraldehyde 3-phosphate dehydrogenase (GAPDH).

Immunohistochemical Staining
We used the protocol described by Ding et al. to stain frozen sections for PRL immunopositivity [46].

Data Analysis
As in previous studies [49,50], we submitted relative transcript abundances to Partek Genomics Suite ® for principal component analysis, and we accepted all principal components with eigenvalues >1.0 as significant. We used SigmaPlot 14 for nonlinear regression analyses, ANOVA for differences between means of normally distributed values, and the Kruskal-Wallis one-way analysis of variance on ranks for differences between non-normally distributed values.

Conclusions
If immune cell cluster/epithelial cell networks analogous to the networks we described develop in human lacrimal glands, the ability to detect their molecular signatures would have important implications for the diagnosis of incipient Sjögren's foci and for the design and implementation of therapies to prevent progression to clinically significant pathology. Funding: This work was supported by NIH Grants EY05801, EY13720, EY 017731, and DK 048522, and by gifts.