Primary cell culture

Primary human keratinocytes used in these experiments were derived from single donor ex vivo explants of neonatal foreskin [53], hysterectomy-derived ectocervix [54], or tonsillectomy-derived palatine tonsils [43]. Foreskin and tonsil samples were provided by Banner University Medical Center Tucson as approved by The University of Arizona Institutional Review Board. Primary cervical cultures were kindly provided by Dr. Aloysius Klingelhutz [54]. The authors did not have access to information that could identify samples during or after data collection We used three independent donor cultures for each tissue origin for these experiments. Isolated primary keratinocytes were maintained as monolayer co-cultures with irradiated mouse J2 fibroblast feeder cells [31]. To preserve their lifespan, primary epithelial cells were conditionally and reversibly immortalized using the Rho-kinase inhibitor Y-27632 at a 10 μM concentration [55, 56]. Y-27632 was removed prior to growing the cells as 3D organotypic raft cultures. Y-27632 has no impact on the ability of primary cells to differentiate using the raft method [56].

3D organotypic raft culture

Full-thickness 3D organotypic epithelial raft cultures were grown over two weeks to yield stratified squamous epithelium [24, 43]. A pair of raft cultures, one for RNA and one for histology, were grown for each of the nine independent donor cultures: n = 3 for each tissue origin; foreskin, cervix, and tonsil. Rafts were constructed by first creating a dermal equivalent, where hTERT-immortalized neonatal human foreskin fibroblasts (HFF-hTERT cells) were embedded (8 x 10⁴ per raft) into a solidified rat tail type I collagen matrix. Three consecutive batches of dermal equivalents were made. On the next day, primary human epithelial cells were seeded at equal proportions (2.5 x 10⁵ per raft) atop the dermal equivalents, with representation from each of the tissue origins across each dermal batch, to create confluent basal layers. Finally, on the third day, cultures were lifted onto a hydrophilic polytetrafluoroethylene membrane insert (0.4 μm pore size) for an air-liquid interface and fed every two days (seven total media changes) with differentiation media (1.88 mM Ca²⁺ and no epidermal growth factor, EGF) to stratify until fully differentiated raft cultures were harvested on day 15 post-lift.

Histology

Whole organotypic epithelial raft cultures were fixed in 10% neutral-buffered formalin for 24 hrs at room temperature then gently washed and stored in 70% ethanol. Fixed rafts were bisected, paraffin-embedded, sectioned (4 μm thickness), and hematoxylin and eosin (H&E) stained at the University of Arizona Cancer Center’s Acquisition and Cellular/Molecular Analysis Shared Resource (TACMASR). Color micrographs were captured using an Olympus BX50 microscope, DP72 camera, and cellSens Entry software. Whole palatine tonsils were collected from Banner University Medical Center Tucson as described above and histologically processed in the same manner as organotypic raft cultures. Overlapping edge color micrographs were captured using an EVOS M5000 microscope and stitched together using Microsoft Image Composite Editor.

Immunofluorescence

Protein-level characterization was performed as described [43]. Unstained slides with formalin-fixed paraffin-embedded epithelial sections (4 μm thick) were deparaffinized by xylenes and ethanol washes. Slides were steamed for 30 min in 10 mM citrate buffer (pH 6), blocked with 2% bovine serum albumin (BSA) in phosphate-buffered saline (PBS) for 30 min, and incubated overnight at 4°C with primary antibodies diluted in blocking buffer. The following primary antibodies and conditions (dilutions and final concentrations) were used: anti-COL17A1 basal epithelial marker (rabbit polyclonal, 1:500, 0.6 μg/mL, MilliporeSigma, Catalog No. HPA043673-100UL) and anti-CRNN upper epithelial marker (goat polyclonal, 1:200, 1 μg/mL, Novus Biologicals, Catalog No. 103242–512). Following PBS washes, slides were incubated in the dark for 30 min at room temperature with secondary antibodies diluted 1:800 in blocking buffer: donkey anti-rabbit AlexaFluor 488 (Thermo Fisher Scientific, Catalog No. A-21206) and donkey anti-goat AlexaFluor 555 (Thermo Fisher Scientific, Catalog No. A-21432). Slides were PBS washed, DAPI-stained (1:1000), ddH₂O rinsed, mounted with ProLong Diamond Antifade Mountant (Thermo Fisher Scientific, Catalog No. P36965), and imaged by epifluorescence using a Leica DMIL LED microscope, DFC3000G camera, and Leica Application Suite X software. For each set of images taken, the DAPI (blue), COL17A1 (green), and CRNN (red) channels were merged in Adobe Photoshop. Reference images were obtained from The Human Protein Atlas (v21.proteinatlas.org) for skin, cervix, and oral mucosa. Quantification of midzone gap distance, which represents the double-negative space between CRNN and COL17A1 layers, was performed in ImageJ/Fiji by measuring perpendicularly from the inferior portion of the CRNN-positive area to the superior portion of the basal layer. Measurements were taken at 25 μm intervals spanning the length of the tissue within each field of view. To determine relative differences due to tissue origin we calculated midzone gap ratios by dividing each midzone gap thickness measurement by the mean of the donor sample of the smallest gap within each set of tissues (T3 for 3D raft cultures and oral mucosa 1505 for the reference tissues).

RNA-seq

RNA was isolated from the epidermal compartment (containing only basal and suprabasal keratinocytes) of organotypic raft cultures using Qiagen’s RNeasy Mini Kit (Catalog No. 74106) after separating them with sterile forceps from their underlying dermal equivalents and flash-freezing with liquid nitrogen [43]. RNA was eluted in nuclease-free water, carry-over DNA was removed with the TURBO DNA-free Kit (Thermo Fisher Scientific, Catalog No. AM1907), then shipped on dry ice to Novogene Corporation (Sacramento, CA). High-quality eukaryotic mRNA libraries were prepared and sequenced on an Illumina platform, yielding >20 million paired-end reads (2x150 bp length) per sample. Raw sequence reads were deposited in the Sequence Read Archive (SRA) as BioProject ID PRJNA916745. High-throughput sequencing data were initially processed using the CyVerse Discovery Environment (http://www.cyverse.org), where FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) was used for quality screening, HISAT2 [57] for mapping sequences to the human reference genome version GRCh38, and featureCounts [58] for generating gene-level read counts. Principal component and differential expression analysis of read counts, including library-size correction and statistical testing, was performed using the DESeq2 Bioconductor package [59] implemented in R (https://www.R-project.org/) via RStudio’s integrated development environment (http://www.rstudio.com/). PcaExplorer was used to plot 95% confidence interval ellipses for the principal components [60]. The 18,925 confidently detected genes were identified using default DESeq2 parameters and were a subset of 31,033 genes with non-zero total read counts among the 58,884 total Ensembl gene IDs used for mapping. Lists of differentially expressed genes, with a stringent adjusted P-value < 0.01, were analyzed for functional enrichment using g:Profiler [61, 62]. Unsupervised clustered heatmaps were constructed in R using the pheatmap package (https://cran.r-project.org/package=pheatmap). Pairwise contrast volcano plots were made using EnhancedVolcano (https://github.com/kevinblighe/EnhancedVolcano). Gene set enrichment analyses (GSEA) and over-representation analysis of pairwise and overlapping genes were done using clusterProfiler [63, 64]. Area-proportional ellipse Euler diagrams were made using eulerr (https://CRAN.R-project.org/package=eulerr). The singscore package [65, 66] was used to calculate tissue-identifying gene expression signature scores. The web-based submission tool for Integrated Motif Activity Response Analysis (ISMARA) was used to identify regulatory motifs and transcription factors that regulate gene expression [67].

Reanalysis of Mostafavi et al., dataset

Mostafavi and colleagues used Affymetrix ST1.0 microarrays to measure the expression of interferon-stimulated genes in mouse B-cells and macrophages [46]. They report the expression of each gene at baseline and 2 hours post interferon stimulation in wildtype (WT) mice. To identify genes that do not require interferon, the experiments were repeated in interferon receptor (IFNAR1) KO mice (KO). Supplemental data S1 provided by Mostafavi and colleagues was downloaded and further analyzed. If a microarray probe was detected in both macrophages and B-cells, the expression values were averaged. Otherwise, the B-cell or macrophage value was used. For some genes, multiple probes were present, for simplicity, only one probe per gene was analyzed. For some genes, the human homologue could not be identified. Similarly, some of the mouse genes were auto-converted to “dates” by Microsoft Excel and could not be recovered [68]. These genes were excluded from our further analyses.

For each gene, a fold change of baseline vs. IFNAR1-KO (WT/KO) expression and the intrinsic responsiveness to interferon (WT+IFN / WT) were calculated. As per Mostafavi and colleagues, highly tonic-sensitive genes are those genes where log₂(foldchange(WT/KO))>log₂(foldchange(WT+IFN / WT))/4. With low tonic-sensitive genes following log₂(foldchange(WT/KO))<log₂(foldchange(WT+IFN / WT))/6. We consider the genes that are larger than log₂(foldchange(WT+IFN / WT))/6 but lower than log₂(foldchange(WT+IFN / WT))/4 to have medium sensitivity.

3DEpiEx application development

The three-dimensional (3D) Epithelial Gene Expression (3DEpiEx) application was developed using the Shiny R package [69]: https://shiny.rstudio.com/. Application development included two main parts: the user interface and the server. The user interface dictates user interaction and the overall appearance of the application, while the server involves handling user inputs and the construction of different outputs. The primary input for 3DEpiEx is the user gene selection. With nearly 19,000 genes from which to choose, it was of the utmost importance to develop an effective tool that would allow users to select specific genes based on Ensembl ID [70] or HGNC gene symbol [71]. Ultimately, the “selectizeInput” search method (https://shiny.rstudio.com/articles/selectize.html), which combines a search bar and a dropdown menu, was utilized. The biomaRt library [72] was used to annotate gene symbols for the pre-existing data which already contained Ensembl IDs. Once the gene selection tool was created, outputs for the box plots and volcano plots were implemented via volcano3D: https://cran.r-project.org/web/packages/volcano3D/vignettes/Vignette.html. The last part of the development process was web hosting and publication via RStudio Connect. This allowed for public access and use of the application, freely available at https://viz.datascience.arizona.edu/3DEpiEx/.

In vitro 3D organotypic epithelia derived from isolated primary cells recapitulate in vivo histology

Primary epithelial cells obtained from foreskin, cervical, or tonsil tissue were cultured to generate 3D organotypic rafts as described [43]. For each tissue, cells from three unique donors were included. These nine 3D organotypic epithelial raft cultures were grown using the identical method and materials as described [43]. All nine 3D organotypic raft cultures were cultured at the same time using the identical reagents. Rafts were processed in sets of three samples, one sample for each tissue. Therefore, any observed differences represent inherent epithelial cell properties (e.g., tissue origin and donor genetics) and are not due to differences in culture methods. Following stratification, we sectioned and stained formalin-fixed tissues with hematoxylin and eosin (H&E; Fig 1A). All cultures (i.e., 3 donors and 3 tissue origins) differentiated into fully stratified 3D organotypic epithelia (S1 Fig). We identified basal keratinocytes at the dermal interface (dotted white line in Fig 1A), suprabasal midzone cells, and differentiated squamous keratinocytes at the apical surface of the tissue. Despite having used identical culturing methods, we observed tissue-specific differences between samples. Compared to foreskin-derived tissues, the suprabasal midzone cells of the mucosal organotypic epithelial tissue (cervix and tonsil) had a greater extent of cytoplasmic clearing. This cytoplasmic clearing indicates high glycogen content [78]. Cytoplasmic glycogen (and other neutral mucins) are rapidly degraded post-harvest and during the formalin fixation process [79] and poorly stained by H&E [80]. Indeed, abundant cytoplasmic glycogen is typical of the normal ectocervical epithelium [81].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Histological and immunofluorescent verification of multi-origin 3D organotypic raft cultures.

(A) The top row represents foreskin epithelia, the middle row represents cervical epithelia, and the bottom row represents oral epithelia. The first column shows H&E micrographs of 3D organotypic epithelial raft cultures. The second column depicts immunofluorescent staining of the 3D raft cultures with COL17A1 (Collagen XVII; green), CRNN (Cornulin; red), and DNA (DAPI-stained; blue). Columns three and four contain immunohistochemical reference staining of COL17A1 and CRNN from The Human Protein Atlas. Colorimetric 3,3’-diaminobenzidine (DAB)-stained COL17A1 micrographs include: normal skin tissue (T-X0500, Female, age 66, Patient id: 4786), normal cervical tissue (T-83000, Female, age 37, Patient id: 4504), and normal oral mucosa tissue (T-51000, Male, age 62, Patient id: 3724) all stained with antibody HPA043673. CRNN micrographs include: normal skin tissue (T-80100, Female, age 66, Patient id: 3357), normal cervical tissue (T-83000, Female, age 55, Patient id: 2005), and normal oral mucosa tissue (T-51000, Male, age 62, Patient id: 1505) all stained with antibody CAB026182. Images are available from v21.proteinatlas.org. Scale bars represent 50 μm. (B) Schematic illustration of the apical, midzone, basal, and dermal layers in stratified squamous epithelia. CRNN-positive cells are typically present in the apical layer (highlighted in red), while COL17A1-positive cells are typically present in the basal layer (highlighted in green). The midzone gap represents the double-negative area from the inferior portion of CRNN to the superior portion of basal epithelial cells. (C) Midzone gap distances were quantified between skin, cervical, and oral epithelia in our 3D rafts and reference Protein Atlas images (patient id indicated for each tissue, v21.proteinatlas.org). Violin plots were made for each independent donor measured while points represent the ratio of each midzone gap distance measurement divided by the lowest donor’s average midzone gap measurement. All measurements were taken at 25 μm intervals spanning the length of visible tissue across non-overlapping fields of view using ImageJ/Fiji. Mean ratios for each donor are indicated by red bars.

https://doi.org/10.1371/journal.pone.0292368.g001

We used immunofluorescence to characterize these tissues further. We stained these tissues with antibodies specific for collagen XVII A1 (COL17A1) and cornulin (CRNN). COL17A1 is expressed in epithelial hemidesmosomes, is a known basal cell marker in stratified squamous epithelial tissue, and is an integral component of the basement membrane [82]. COL17A1 protein was detected in cells localized to the basal layer that contact the basal membrane at the dermal interface (Fig 1A; green). Suprabasal cornulin is a terminal differentiating cell marker in the superficial layers of differentiated epithelial tissue [83, 84]. We detected cornulin in the superficial layers of 3D organotypic epithelial raft cultures (Fig 1A; red). All tissues have a section of double-negative midzone cells (COL17A1^NEG/CRNN^NEG) representing a clear distinction between the basal COL17A1(+) monolayer and apical CRNN^POS layers of these tissues (Fig 1B). Comparing the height of the midzone gap between the tissues (Fig 1C), the gap was most prominent in foreskin-derived 3D epithelia (32 +/- 11 μm, n = 3 tissues), moderate in cervix-derived cultures (17 +/- 4 μm, n = 3 tissues), and smallest in tonsil-derived cultures (10 +/- 2 μm, n = 3 tissues). We identified a similar pattern of COL17A1^NEG/CRNN^NEG midzone gap heights (Fig 1C) in reference micrographs obtained from The Human Protein Atlas (v21.proteinatlas.org) tissue database [85]. The midzone gap measured between the COL17A1^POS basal cells and the CRNN^POS apical cells was most prominent in the skin (69 +/- 10 μm, n = 2), moderate in the cervix (29 μm, n = 1), and smallest in the oral mucosa (18 +/- 6 μm, n = 2). Since 3D organotypic rafts have a reduced total epithelial thickness compared to in vivo tissue, we compare the midzone gap thickness to the oral raft or tissue. The midzone gap in skin-derived tissues is roughly 3 times as wide as those in oral-derived tissue (3.1x 3D rafts and 3.9x for skin). For both cervical rafts and cervical tissue the midzone gap is 1.6x wider than the oral equivalent (Fig 1C).

Tonsillar epithelia are composed of reticular crypt epithelia and surface epithelia (Fig 2) [86]. The reticular crypt epithelium likely plays an important role in regulating the immunological functions of the tonsil [87]. The crypt epithelium also has reduced barrier function, while the surface epithelia, like the surrounding oral mucosa, is thought to serve a typical epithelial barrier role owing to tight junctions [88]. When cultured as 3D organotypic rafts, two of the tissues (T1 and T2) display reticular/crypt epithelial characteristics, including a reduced overall thickness and extent of stratification and a spongey-like appearance (Fig 2). The T3 sample more closely resembles surface epithelium, with more layers of stratified squamous cells and signs of apical keratinization/cornification. Compare the micrograph derived from a tonsillectomy to the H&E-stained 3D organotypic rafts (Fig 2): focusing on the thickness differences between crypt and surface epithelia, as well as the extent of stratification. Also, note the extent of lymphocyte infiltration in the ex vivo crypt epithelia and the reticulated “spongey-like” appearance.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. 3D organotypic raft cultures recapitulate differences between crypt- and surface-like tonsillar epithelia.

(A) Histology of crypt-like and surface-like tonsillar epithelia. The top row contains H&E micrographs of three donor-derived 3D tonsillar epithelial raft cultures (T1, T2, and T3). Below we show micrographs from real tissue: a stitched composite of the human palatine tonsil, which includes different types of tonsillar epithelia: crypt epithelia on the left and surface epithelia on the right. The surface epithelium of tonsils consists primarily of stratified squamous cells and makes up the outer periphery, whereas crypt epithelium serves to increase the overall surface area and is composed of an uneven mixture of stratified squamous cells and non-epithelial cells (e.g., infiltrating lymphocytes), resulting in a spongy appearance. As a result of these epithelial subtypes, there is a variable expression of differentiation-related epithelial markers in tonsillar tissues. Black scale bars represent 50 μm and the red scale bar represents 500 μm. (B) Transcriptional variation between surface-like (T3) and crypt-like (T1, T2) 3D tonsillar epithelia: top 25-varying genes within tonsil samples. (C) Over-representation analysis of surface-like tonsillar epithelia markers (REAC = Reactome pathway).

https://doi.org/10.1371/journal.pone.0292368.g002

In conclusion, 3D organotypic raft cultures recapitulate several histological features of the in vivo tissues from which the cells were initially derived.

Epithelial origin is the primary determinant of transcriptomic variability between foreskin, cervix, and tonsil-derived 3D rafts

We analyzed the transcriptome of 3D organotypic epithelial tissues using RNA-seq (S1 and S2 Files) to compare and contrast foreskin, cervix, and tonsil-derived tissues. We carefully peeled the epidermis off the dermal equivalent and isolated total RNA from these cells. Poly-A enriched cDNA was sequenced by a commercial provider using an Illumina NovaSeq. Table 1 shows that sequencing depth (average ~29 million paired-end reads per sample) and quality (~98% of sequences had at least a Q30 quality score) were sufficient for further analyses [89, 90]. Roughly 97% of the paired-end sequencing reads for each sample were mapped to the human reference (GRCh38, Table 2) using HISAT2 [57]. The library-size correction was performed using DESeq2 [59] (S3 File). Cook’s distances for all samples were inspected via violin plots to ensure no individual samples were overtly contributing to outliers (S2A Fig). From the 31,033 genes with detectable counts, 12,038 were excluded due to low counts (mean count < 3), and 70 were excluded as count outliers (based on Cook’s distance, S2B Fig).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. RNA-seq file quality summary.

https://doi.org/10.1371/journal.pone.0292368.t001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. RNA-seq mapping and gene-level assignment summary.

https://doi.org/10.1371/journal.pone.0292368.t002

As an initial analysis, we used principal components to analyze the top 500 varying genes (based on gene-level counts). Principal component analysis revealed that each tissue had a distinct transcriptional profile (Fig 3). The scree plot in Fig 3A shows that Principal Component 1 (PC1; 41%) and PC2 (28%) explain 69% of the variance. PC1 and PC2 visually distinguish foreskin (green), cervix (orange), and tonsil (blue) derived tissues from each other (Fig 3B), further supporting that 3D organotypic rafts maintain tissue-specific features. The density plots for each tissue along both PC1 and PC2 demonstrate the spatial separation of these transcriptomes along the axes of the individual principal components. We used pcaExplorer [60] to draw 95% confidence interval ellipses around each set of samples belonging to a specific tissue (Fig 3B). The lack of overlap between the ellipses statistically supports that foreskin, cervix, and tonsil-derived tissues have different transcriptional profiles. This separation is further confirmed using a multivariate analysis of variance (MANOVA; P-value = 2.89 x 10⁻⁵).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Transcriptomic characterization of multi-origin 3D organotypic raft cultures.

(A) Scree plot of principal components (PCs) and their variance explained (%), totaling 100%, using the top 500 varying genes (variance-stabilized counts) detected by RNA-seq. PC1 (41%) and PC2 (28%) explained the majority (~70%) of transcriptional variability. (B) Principal component analysis (PCA) plot of 3D organotypic raft culture transcriptomes (top 500 varying genes, variance-stabilized counts), each point represents an individual raft culture from an independent donor, biological sex labeled (♂, male; ♀, female), with 95% confidence ellipses for the tissue origin groups. Density plots for PC1 and PC2 are plotted opposite the respective axis. (C) Verification of biological sex by sex-specific marker XIST. DESeq2 normalized (library-size corrected) counts for each sample confirmed sample identity.

https://doi.org/10.1371/journal.pone.0292368.g003

Importantly, these analyses demonstrate that the tissue origin of keratinocytes, rather than the genetics of the individual donor, is the primary driver of transcriptomic variability of 3D organotypic raft cultures. Of note, biological sex, as determined by the sex-specific expression marker XIST (Fig 3C), is not a significant contributor to transcriptional variability among our samples, as T1 (♂, male) clusters closely to T2 (♀, female) (Fig 3B).

The densities along PC1 separate each of the tissue groups. Therefore, we conclude that PC1 (41% of the variability) separates the three tissues, whereas PC2 has overlapping distributions and does not clearly distinguish foreskin (green) from tonsil (blue) tissue. This increased spread is likely driven by variability in the tonsillar tissue, i.e., T1/T2 (crypt) vs. T3 (surface). To further investigate this, we plotted the top 25 variable genes between the different tonsil tissues (Fig 2B). T3 shows higher expression of genes implicated in tissue keratinization (e.g., CRNN, SBSN, and SPINK7). Indeed, functional profiling of the 19 genes up-regulated in T3 indicates that these genes are implicated in the cornified envelope formation and keratinization (Fig 2C). Likewise, T1 and T2 express the follicular dendritic cell-secreted protein (FDCSP) ~40x higher than T3. FDCSP was previously reported as a marker for crypt epithelia [91, 92]. These data further support the hypothesis that T1/T2 and T3 are derived from different sub-tissues and further confirm that 3D organotypic rafts maintain many features associated with the original, in vivo tissue. Nonetheless, as is clear from Fig 3B, the tonsillar tissues are transcriptionally more similar to each other than to other tissues. Therefore, we will analyze T1, T2, and T3 as biological replicates.

Global gene expression across three different tissue-derived 3D rafts

To characterize the tissue-inherent transcriptional differences, we used a likelihood ratio test (LRT) to detect the globally differentially expressed genes between foreskin, cervix, and tonsil-derived epithelial. Based on our biological question and supported by the exploratory PCA of the primary sources of variation (Fig 3B), we included tissue origin as the grouping factor, with three levels (foreskin, cervix, and tonsil) in the DESeq2 model [59]. The individual donor samples were used as biological replicates (n = 3). Using a strict statistical significance threshold (P-adj < 0.01), we identified 1,238 globally differentially expressed genes distinguishing the three tissues (S4 File). The combination of different donors and strict statistical cut-off increases the likelihood of identifying biologically relevant differentially expressed genes. These data demonstrate important transcriptional differences between foreskin, cervix, and tonsil-derived tissues. In turn, we may miss important differentially expressed genes due to stringent cut-offs.

We used unsupervised clustering to analyze these differentially expressed genes (Fig 4). Unsupervised clustering builds a hierarchy of clusters, which are summarized as a dendrogram. In these heatmaps, individual columns represent the different samples, while differentially expressed genes are shown in rows. The color indicates the relative expression of each differentially expressed gene in the different samples: blue means the gene is expressed less in that sample compared to the normalized gene expression across all samples, whereas red means the gene is overexpressed. Therefore, expression levels should only be compared across samples for the same gene and not between genes. Unsupervised clustering of the samples confirmed that all three tissues are different. Based on this set of genes, the mucosal tissues (cervix and tonsil) are transcriptionally more similar to each other than to foreskin-derived tissues. Hierarchical clustering of the genes identified two main clusters we designate as "HFK-up" and "HFK-down", respectively (see S5 File for human foreskin keratinocyte, HFK, cluster annotations).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Globally differentially expressed genes.

Variance-stabilized, row-centered, counts are shown for all 1,238 genes differentially expressed genes (LRT, P-adj < 0.01). For each cluster, the top five over-represented Gene Ontology (GO) terms for each GO category (BP, biological process; MF, molecular function; and CC, cellular component) are listed based on the lowest g:Profiler P-adj.

https://doi.org/10.1371/journal.pone.0292368.g004

We identified 535 genes up-regulated in foreskin-derived tissues (HFK-up cluster). On the other hand, HFK-down contains 703 genes with reduced steady-state levels compared to the mucosal tissues. Over-representation analysis of the two clusters via g:Profiler (ordered query, g:SCS threshold < 0.05) identified 181 functionally enriched terms for HFK-up (S6 File) and 293 functionally enriched terms for HFK-down (S7 File).

HFK-up genes were enriched for biological processes related to pattern specification and morphogenesis. These embryogenic terms relate to tissue development along an organism’s anterior-to-posterior axis [93]. The top enriched cellular components for HFK-up include differentiation-related structures: the cornified envelope and lamellar bodies common to the upper layers of stratified squamous epithelium that are important for barrier maintenance [94]. The top enriched up-regulated molecular functions in HFKs relate to transcription factor activity. Furthermore, the top enriched cellular components for HFK-up also included chromatin and chromosome subcellular localizations, further supporting the critical role of these DNA-binding transcription factors. While we expected tissue-specific differences in differentiation and transcription factor activity, we were intrigued to find an enrichment of innate immune-related genes and transcription factors supporting the identification of "response to type I interferon" as a top enriched biological process in the HFK-up cluster (Fig 4).

HFK-down genes were enriched for differentiation-related biological processes (e.g., development and morphogenesis), further highlighting differences between the tissues related to differentiation. The analyses also identified differences related to cell adhesion processes. Of note, these cell adhesion differences are also reflected by the observations that the top enriched cellular components included the extracellular matrix and cell periphery. Top molecular functions enriched in the HFK-down group are related to growth factor binding and receptor signaling (Fig 4). Enrichment of these genes suggests that epithelial cells actively contribute to shaping the extracellular matrix and that these tissues differentially organize the extracellular microenvironment. This likely has important consequences for tissue homeostasis.

An online database allows for easy analysis of differential gene expression

This paper describes a careful transcriptional analysis of 3D organotypic rafts derived from different tissues. To maximize the impact of this resource, it is essential to share these data with the community, both in the form of raw data that allows for reanalysis and in the form of a pre-analyzed database that allows for easy reference. We developed an interactive database using the R statistical language and Shiny app, available at https://viz.datascience.arizona.edu/3DEpiEx/. The website allows users to visualize the expression of individual genes of interest across the different tissues (screenshot Fig 5A) as well as in 2D volcano plots (screenshot Fig 5B, similar to the plots in Fig 6) and an interactive 3D volcano plot linking differential expression and statistical support.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. 3DEpiEx interactive database.

Screenshots of our web-based three-dimensional (3D) Epithelial Gene Expression (3DEpiEx) application which allows users to search genes of interest and visualize their expression in foreskin, cervix, and tonsil-derived 3D epithelial raft cultures as (A) boxplots of normalized or variance-stabilized counts and (B) volcano plots in the context of pairwise contrasts.

https://doi.org/10.1371/journal.pone.0292368.g005

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Pairwise contrasts identify the differentially expressed genes between foreskin, cervix, and tonsil-derived 3D epithelia.

(A) Euler diagram with area-proportional ellipse for each pairwise contrast, including up (red) and down (blue) genes. (B-D) Volcano plots of the DESeq2 pairwise contrast results. The top five up- and down-regulated genes are labeled. (E-G) Ridge plots of the aggregated core enriched gene distribution (log₂ fold change) from pairwise gene set enrichment analysis (GSEA): non-redundant functional modules from network analysis of the statistically significant Gene Ontology Biological Process (GO:BP) terms are shown (S3–S5 Figs). The number of core enriched (leading edge) genes are shown beside each distribution. Red = enriched. Blue = depleted.

https://doi.org/10.1371/journal.pone.0292368.g006

Pairwise contrasts between 3D tissues highlight their distinct transcriptional and functional profiles

Due to their availability, primary human foreskin keratinocytes are often used in research [74]. Our global analysis identified that the foreskin samples clustered apart from the cervix and tonsil tissues, indicating that these skin-derived cells are phenotypically and transcriptionally distinct from each of the mucosal-origin tissues studied here. Indeed, the two primary clusters of differentially expressed genes corresponded to down- or up-regulation in the foreskin (Fig 4).

We updated the described DESeq2 model to make direct pairwise comparisons between each of the three epithelial tissue origins [cervix vs. foreskin (S8 File), tonsil vs. foreskin (S9 File), and tonsil vs. cervix (S10 File)]. Among the 18,925 confidently detected genes, we found a total of 1,381 pairwise differentially expressed genes (Wald test P-adj < 0.01). These 1,381 identified genes were plotted using a three-way area-proportional Euler ellipses diagram (Fig 6A).

The Euler diagram in Fig 6A illustrates that many identified pairwise differentially expressed genes are shared across pairwise tissue comparisons. Each pairwise comparison was visualized using individual volcano plots. These plots show the relationship between the magnitude of differential expression (x-axis) and statistical significance (y-axis) for each differentially expressed gene (Fig 6B–6D). Functional processes enriched and depleted for each comparison were identified using gene set enrichment analyses (GSEA). Individual genes often belong to multiple GO terms, leading to overlap between distinct GO terms. To deconvolve this redundancy, we visualized the statistically significant GO terms using enrichment maps and gene-concept networks (S3–S5 Figs). These networks organize enriched terms with edges connecting overlapping gene sets to identify functional modules. For each network, we identified functional modules and plotted the distribution of fold change values of all core enriched (also known as leading edge) genes contained within each module (Fig 6E–6G). Therefore, these ridge plots show the number of aggregated core genes in each module and their differential expression. Red curves show that a set of core genes is enriched, while blue distributions indicate depletion.

Tissue-identifying gene expression profiles highlight differences in embryonic origin and interferon responses

To investigate tissue-specific gene expression, we analyzed the pairwise overlapping regions of the Euler diagram (Fig 6). For example, Fig 7A shows genes differentially expressed in foreskin-derived raft tissues compared to both cervical and tonsillar tissue (i.e., the overlap between the tonsil vs. foreskin and cervix vs. foreskin comparisons, highlighted in green in Fig 7A). We also considered differentially expressed genes in all three comparisons (38 genes in the black center of the Euler diagram in Fig 6A). This set of genes is plotted using a scatterplot. In these plots, genes in the top right quadrant are up-regulated in foreskin tissue, while genes that fall in the lower left quadrant are down-regulated in foreskin compared to the other tissues. Transcripts in the other two quadrants are differentially expressed in the three-way overlap but are incongruent between comparisons. For example, the expression of ESR1 is statistically different in all three comparisons. However, ESR1 is up-regulated in foreskin-derived tissues compared to tonsil, while ESR1 is down-regulated when comparing foreskin to cervical tissues (Fig 7G and 7H). These incongruent transcripts are excluded from further analyses. Genes in the top right (i.e., up-regulated) and bottom left (i.e., down-regulated) quadrants are further analyzed using a combination of g:Profiler and network analyses to identify interconnected functional modules.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Tissue-identifying differentially expressed genes.

(A-C) Scatterplots of pairwise contrast overlaps, axes represent log₂ fold change. (D-F) Over-representation analysis enrichment plots for each tissue-identifying gene set, showing the non-redundant functional modules identified via network analysis (S6–S8 Figs). (G) Heatmap of the core distinguishing epithelial markers (n = 38) of the triple-overlap differentially expressed genes: six groups identified based on tissue expression. The color scale from blue to red represents variance-stabilized, row-centered, counts. (H) Log₁₀-scaled normalized RNA-seq count data were plotted for genes from each of the six groups.

https://doi.org/10.1371/journal.pone.0292368.g007

A core set of differentially expressed genes uniquely distinguish each tissue

The three-way overlap identifies 38 differentially expressed genes in all three pairwise comparisons (Fig 6A). These 38 genes are differentially expressed in all pairwise comparisons and therefore have a unique expression level for each of the three epithelia. Visualizing the expression of these 38 genes in a heatmap (Fig 7G) highlights the different relative expression levels in each tissue. These 38 genes can be grouped into 6 clusters based on relative expression levels. Clusters 1 and 2 contain genes that are up-regulated in tonsils and vary in expression between foreskin and cervix. Clusters 3 and 4 are up-regulated in cervix with differential expression in foreskin and tonsil. Finally, clusters 5 and 6 are up-regulated in foreskin-derived tissues compared to the other 2 tissues.

We highlight 6 genes using plots of the normalized counts to show the magnitude of expression (Fig 7H). KRT10, a keratin-encoding gene, was among this small subset of genes. KRT10 was highest in foreskin (Wald P-adj = 3.14 x 10⁻¹⁸ and 3.17 x 10⁻⁴⁴, relative to cervix and tonsil, respectively), followed by cervix (Wald P-adj = 4.91 x 10⁻⁵, relative to tonsil), and lowest in tonsil (Fig 7H). This keratin is commonly used as an epithelial differentiation marker [43], found in the spinous layer of cutaneous epithelium [97], but clearly has a tissue-specific abundance.

The stratum corneum of the epidermis Is a critical barrier that functions to limit water loss by evaporation. The aquaglyceroporin AQP3 is abundantly expressed in keratinocytes of the mammalian skin epidermis along with the stratified squamous epithelial cells of the oral cavity and esophagus [98, 99]. AQP3 is highest in foreskin compared to both cervix (Wald P-adj = 2.9 x 10⁻¹⁸) and tonsil (Wald P-adj = 1 x 10⁻³), while AQP3 expression is higher in tonsils compared to cervical tissues (Wald P-adj = 9.1 x 10⁻⁶).

The estrogen receptor (ESR1) was up-regulated in tissues with a genital origin. Estrogen receptors are highly expressed in the ectocervix [100] and foreskin [101], and we found it was highest in cervix (Wald P-adj = 5.82 x 10⁻⁷ and 9.89 x 10⁻²⁰, relative to foreskin and tonsil, respectively), followed by foreskin (Wald P-adj = 1.40 x 10⁻³ relative to tonsil), and lowest in tonsil (Fig 7H).

MYEOV (MYEloma Overexpressed gene) was initially identified as a potential oncogene in specific multiple myeloma cell lines [102]. Importantly, it has also been implicated in oral squamous cell carcinomas [102]. We find that it has the highest expression in cervix (Wald P-adj = 3.4 x 10⁻¹² and 3.8 x 10⁻³, relative to foreskin and tonsil, respectively), followed by tonsil (Wald P-adj = 2 x 10⁻³ relative to foreskin) and foreskin derived tissues.

Family with sequence similarity 3 (FAM3) is a cytokine-like gene family that contains four genes: FAM3A, FAM3B, FAM3C, and FAM3D [103]. FAM3B (also called PANDER) regulates the epithelial-mesenchymal transition and has been implicated in the development of oral cancer [104]. In our data, FAM3B (and its family member FAM3C) are highly expressed in tonsillar cells and lowest in cervix (Tonsil vs. Foreskin: Wald P-adj = 1.2 x 10⁻⁷; Tonsil vs. Cervix Wald P-adj = 1.9 x 10⁻¹⁷; Foreskin vs. Cervix: Wald P-adj = 5.8 x 10⁻⁵).

Compared to foreskin-derived tissues, FOXE1 (Forkhead Box E1) is up-regulated in both tonsil (Wald P-adj = 8.06 x 10⁻¹¹¹) and cervix (Wald P-adj = 9.04 x 10⁻³¹). FOXE1 is a member of the forkhead family of transcription factors. FOXE1 is primarily considered a thyroid transcription factor that plays a role in thyroid morphogenesis [105]. Within the mucosal tissues, FOXE1 is further up-regulated in tonsils (Wald P-adj = 1.69 x 10⁻²⁵) relative to the cervix (Fig 7H). During embryogenesis, the Pharyngeal endodermal pouches give rise to tissues responsible for forming the middle ear cavity and eustachian tube, palatine tonsils, thymus, parathyroid glands, and parafollicular cells of the thyroid, potentially explaining the observed expression profile of FOXE1 in tonsils.

In vitro-derived tissue-distinguishing gene signatures are shared with in vivo tissues

In the previous sections, we used our 3D organotypic rafts to characterize tissue-identifying gene expression signatures (Fig 7). If these in vitro tissues (and thus their gene expression signatures) are equivalent to the ex vivo samples, we hypothesized that these in vitro gene expression signatures should allow for the identification of the ex vivo samples.

The Genotype-Tissue Expression (GTEx) project is a resource database containing tissue-specific RNA-seq data primarily obtained from low-post-mortem-interval autopsies [106]. We downloaded gene-level read counts of 2,183 samples (skin, oral, or cervical epithelial) from GTEx. Six hundred four samples were reported to be “suprapubic skin” (not sun exposed), 701 samples were “lower leg” skin (sun-exposed), 9 samples were “ectocervix”, 10 samples were “endocervix”, 142 samples were “uterus”, 555 samples were “esophagus mucosa”, and 162 samples were “minor salivary gland” (from the inner surface of the lower lip). For each of these datasets, we used the “tissue-distinguishing” genes plotted in the top right and bottom left quadrants of Fig 7 to calculate rank-based single-sample signature scores. The higher this score, the more the sample resembles our 3D organotypic tissues (indicated by the dotted line). When using the in vitro foreskin-derived signature, 99.2% of the GTEx samples identified as “skin” have positive scores similar to those observed for the in vitro tissues (dotted line). On the other hand, 84.3% of the other samples scored negatively. Notably, while some of the other tissues score positive on this scale, they do not compare to the magnitude of the “skin” tissue scores (Fig 8A). We observed similar results for the in vitro cervix-derived signature, which positively identified 97.5% of uterine, ectocervical, and endocervical samples in our dataset, while 99.9% of the other tissues scored negative (Fig 8B). Likewise, our in vitro tonsil-derived signature positively identified 98.2% of oral mucosa samples (Esophagus and Salivary gland). Notably, all the other samples scored negative (Fig 8C).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Tissue-distinguishing epithelial signatures identify samples from the GTEx database.

Tissue-derived 3D epithelial signature scores were calculated for n = 2,183 GTEx RNA-seq samples using our tissue-distinguishing overlap gene sets. We accurately score the samples based on their anatomic site (A) foreskin identifies skin samples, (B) cervix identifies female reproductive tract samples, and (C) tonsil identifies oral samples. The dotted line represents the pure epithelial signature score of our 3D raft cultures.

https://doi.org/10.1371/journal.pone.0292368.g008

Overall, these findings indicate that our pure 3D epithelial in vitro signatures can classify complex in vivo tissue samples accurately based on anatomic origin.

Reduced homeostatic interferon-stimulated gene expression in tonsil cells

Intriguingly, baseline interferon-stimulated gene expression is reduced in tonsillar cells compared to other tissues (Fig 7E). We used Integrated Motif Activity Response Analysis (ISMARA) to reconstruct the regulatory networks that control the observed gene expression in these tissues [67]. ISMARA identifies key transcription factors that explain the observed gene expression patterns [67]. The ISMARA analysis identified 4 transcription factor binding motifs with a z-score larger than 2 whose target gene expression does not change between foreskin and cervix but is reduced in tonsil-derived tissue. These transcription factor binding motifs are “IRF2_STAT2_IRF8_IRF1”, “IRF9”, “IRF3”, and “PITX1”. Notably, the baseline expression of these transcription factors is not different between the different tissues (S9 Fig), suggesting that upstream pathways, and not just transcription factor abundance, must be differentially regulated in tonsillar epithelial cells.

A recent study shows an important contribution of STAT2–IRF9 complexes for the constitutive expression of interferon-stimulated genes in resting cells [107]. This subset of interferon-responsive genes is involved in a homeostatic response and is known as ‘tonic interferon-sensitive’ genes. These tonic interferon-sensitive genes were characterized by Mostafavi and colleagues in B-cells and macrophages [46]. We reanalyzed the Mostafavi dataset and classified their set of interferon-stimulated genes (n = 188) as having high (27 genes), medium (15 genes), or low (146 genes) tonic sensitivity (Fig 9A). Genes with high tonic sensitivity have a higher baseline of unstimulated expression compared to medium-sensitive genes. For example, IFIT3 and SLFN5 are "high" tonic-sensitive genes; IFIT3 is differentially expressed in tonsillar epithelia, while SLFN5 does not reach statistical significance.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. Tonsilar epithelia down-regulate tonic-sensitive interferon genes.

(A) Scatterplot showing Interferon Stimulated Gene expression vs. dependence on baseline tonic interferon signaling (data aggregated from [46]). For details on how data was reanlyzed and aggregated, see Materials and methods section. Symbols indicate tonic sensitivity (circle: high, square: medium, and triangle: low), while colors indicate whether the gene is significantly (red: p < 0.01, blue: p > 0.01) differentially expressed in tonsils vs. cervix and foreskin. (B) Heatmap of the Mostafavi gene set with rows (genes) sorted by fold change of tonsils relative to the other two epithelia. The black and grey bars indicate high and medium tonic sensitivity, respectively. (C) Correlation plots of genes with either tonic sensitivity [high (circles) + medium (squares)] or low sensitivity, stratified by whether genes are differentially expressed in tonsil relative to the other two epithelia. R-squared value and p-value for the linear regression line are shown. N indicates the number of genes in each set.

https://doi.org/10.1371/journal.pone.0292368.g009

Similarly, 11/27, 4/15 medium, and 21/146 low sensitivity genes are differentially expressed in tonsils. We used a Chi-square test to determine that differentially expressed genes are enriched among high vs. low sensitivity tonic genes (P < 0.0001). Likewise, high-sensitivity genes are more likely to be differentially expressed than genes with medium sensitivity (P = 0.043). This suggests that tonic-sensitive interferon-responsive genes are specifically down-regulated in tonsils. To further investigate this connection, we plotted the Mostafavi gene set (n = 188) using a heatmap (Fig 9B). The rows (i.e., genes) in this heatmap are sorted by fold change comparing tonsils to the other two epithelia, with the gene with the biggest reduction (CCL2) on the top and the gene with the largest increase in expression (GBP6) on the bottom. Tonic-sensitive groupings are indicated to the left of the heatmap. The top 10% down-regulated genes (n = 18) contain 8 highly sensitive, 2 medium, and 8 lowly sensitive genes. Therefore, the top 10% of down-regulated genes contain 23.8% of the tonic-sensitive (high and medium) genes. Our data suggest tonic-sensitive genes are specifically down-regulated in tonsillar epithelia, potentially limiting interferon responses due to constant exposure to viruses and microbes in the oral cavity [108]. Following that logic, we hypothesized that tonsillar epithelia might specifically reduce those tonic-sensitive genes that have the highest induction following interferon stimulation. For each gene in Fig 9B, we plotted response intensity after interferon treatment (derived from [46]) against baseline expression in tonsils (fold change vs. other epithelial tissues). For statistically differentially expressed genes with high (circles) or medium (squares) tonic interferon sensitivity (Fig 9C), tonic sensitivity is inversely correlated with baseline expression. This suggests that, in tonsillar epithelia, genes that are highly stimulated by exposure to interferon are more repressed at baseline. No correlations were observed for genes with low tonic sensitivity. Taken together, these data suggest that tonsillar epithelia down-regulate a subset of interferon-sensitive genes, and this repression is correlated with the gene’s expression following interferon stimulation.

Donor-derived epithelial tissues retain origin-specific differences and are valid models

Overall, we found that 3D organotypic epithelia retain their origin-specific characteristics, with anatomic origin sites dominant over donor variability. This was evident at the tissue (histology and immunofluorescence) and transcriptional levels. Overall, many of the tissues’ differences are related to homeotic genes. These genes, e.g., HOX and PAX cluster members, are master regulators directing the development of particular body segments or structures [112]. We found homeotic genes among the top differentially expressed genes in each pairwise contrast, with HOXC10 highly expressed in foreskin-derived epithelia, HOXA11 in cervix-derived epithelia, and PAX9 in tonsil-derived epithelia. HOXC10 is expressed during mouse urogenital development [113] in the analogous tissues that would give rise to the human foreskin. HOXA11 is expressed in the developing cervix [114]. PAX9 is essential for tissues developing from the pharyngeal pouches [115], including the palatine tonsil. These in vitro 3D epithelial tissues retain the expression of tissue-specific homeotic regulators, supporting the utility of using epithelial cells from distinct origins in 3D models, implying that these regulatory genes may be valuable biomarkers for distinguishing tissues and verifying their provenance. Importantly, we used gene expression signatures derived from these in vitro 3D tissues to classify samples from the GTEx database. Therefore, these pure 3D tissues maintain many transcriptional and histological features of complex ex vivo tissue samples.

Foreskin-derived cells are imperfect models for mucosal tissues

Due to the relative ease of access to foreskin-derived keratinocytes [116], these cells have become the de facto gold standard of stratified squamous epithelial tissue culture. However, our data clearly demonstrate that, depending on the research question, there may be better models than these cells. Indeed, based on global gene expression analysis, mucosal tissues are transcriptionally distinct from foreskin-derived tissues. This further supports efforts by us and others to transition to more relevant tissues during experimentation [34, 35, 39, 43, 117–119].

Based on gene ontology analyses, many of the differences between tissues are related to "keratinocyte differentiation." Nonetheless, our histology and immunofluorescence data clearly demonstrate that all three tissues support terminal differentiation.

The Gene Ontology Consortium attempts to create a nomenclature describing how a gene product may regulate biology. Almost by definition, GO annotations rely on the source and data to support an asserted association between a gene product and a GO definition [120]. Therefore, since most research on keratinocyte differentiation is based on skin-derived cells, GO terms related to “keratinocyte differentiation” will be biased to skin-derived tissues and will, artefactually, propose that differentiation and associated terms are reduced in these mucosal tissues. Studies like these focused on characterizing these tissues can help update and continue to improve Gene Ontology and related resources.

Additionally, differences within specific epithelia, such as tonsil epithelia, further distinguish epithelia from one another. Clearly, utilizing foreskin-derived cells or 3D epithelia to study host-pathogen interactions or diseases that target specific epithelia elsewhere in the body may not accurately reflect the impact those interactions have on specific genes that are expressed or under-expressed in that particular epithelia.

Tonsil cells have a reduced baseline immune response

Tonsils are a part of the mucosa-associated lymphoid tissue (MALT), protecting the entrance to the respiratory and alimentary tract [121]. As such, tonsils are continuously exposed to microbiota and pathogens. Cells and tissues must balance homeostasis and inflammation. Therefore, tonsils must be able to discriminate between potentially infectious pathogens and innocuous airborne and food antigens. We demonstrate that tonsillar epithelia have reduced baseline expression of interferon-sensitive genes compared to cervical and foreskin tissues. Specifically, this reduction is predicted to be driven through STAT2 and IRF9. Upon type I interferon activation, STAT1, STAT2, and IRF9 form the ISGF3 complex to drive the expression of the full spectrum of interferon-stimulated genes. However, preformed STAT2-IRF9 complexes control the basal expression of interferon-stimulated genes [107]. This basal interferon expression is known as tonic interferon sensitivity [46, 122, 123]. We show that tonsil cells down-regulate these tonic-sensitive genes’ baseline expression. Furthermore, those genes that are strongly induced upon exposure to interferon are the most strongly repressed, presumably increasing the activation energy to lead to transcriptional activation. Of note, STAT2 or IRF9 are not differentially expressed in tonsils, suggesting that other mechanisms must be regulating these pathways. The nature of these control mechanisms remains to be determined but appears to be inherent to explanted tonsil cells as it does not require interaction with other cells or microbes. It is tempting to hypothesize a role for chromatin modifications [124].

These data suggest that tonsillar epithelia may have a reduced sensitivity to interferon. Interestingly, it was demonstrated that the microbiota continuously induces low-level interferon and that this baseline stimulation partially regulates the tonic activity of the interferon signaling pathways [108]. Since tonsils are constantly exposed to microorganisms in the oral cavity, it is tempting to speculate that reduced sensitivity to type I interferon could be beneficial. Interestingly, high concentrations of type I interferons may block B-cell responses [125]. During an adaptive immune response, germinal centers develop within the tonsils. B cells congregate in the germinal centers and differentiate into long-lived antibody-secreting plasma cells and memory B cells. While interferon stimulates the development of germinal centers, dysregulated germinal center responses may lead to the development of auto-immune disease [126–129]. Indeed, interferon signaling and improper B-cell activation may promote lupus pathogenesis [130]. Furthermore, it is well known that excessive interferon signaling has immunosuppressive effects and/or causes inflammation and tissue damage that exacerbate disease. Our data suggest that tonsil epithelia evolved to shape the local immune environment to avoid excessive immune activation.

Indeed, mice without a microbiome due to antibiotic treatment do not mount appropriate innate and adaptive antiviral immune responses following exposure to different viruses [108, 131]. Furthermore, these mice had severe bronchiole epithelial degeneration after influenza virus infection [108, 132], implicating an important role for tonic interferon signaling in balancing inflammation and immune responses. Indeed, it was previously demonstrated that tonsillar dendritic cells demonstrate a reduced type I interferon response to influenza infection compared to dendritic cells isolated from blood [49]. Remarkably, different bacteria (S. aureus, H. influenzae, and N. meningitidis) activate a robust type I interferon response in freshly purified blood dendritic cells [48]. However, when the same dendritic cells were exposed to bacteria in media conditioned by tonsillar epithelia, the interferon response decreased 6-fold compared to control media [48]. This suggests that tonsillar epithelia produce paracrine (and potentially autocrine) factors that regulate the interferon response in dendritic cells. Importantly, these tonsillar epithelial “educated” dendritic cells still support T-cell differentiation and T-cell-derived cytokine production [48]. This suggests that tonsillar epithelia signal dendritic cells to become less inflammatory in response to constant microbiota exposure while supporting robust adaptive immunity. However, unlike what was proposed by the authors [48, 50], we do not think that reprogramming the local immune microenvironment is a shared feature of all mucosal epithelia. Our data suggest that tonsillar epithelia may be unique in this by reducing the baseline expression of specific interferon-stimulated genes.

This differential regulation of interferon-stimulated genes may partly explain a puzzling result. Chatterjee and colleagues used microarrays to compare the gene expression profile of different tissues infected with HPV16 [34]. The authors noted that unlike in foreskin-derived tissue, HPV16 did not reduce the expression of immune-related genes in tonsil-derived tissue. However, this study was not designed to directly compare uninfected foreskin and uninfected tonsil tissue [34]. Our data suggest that the expression of these immune-related genes is already reduced in tonsils. Importantly, this raises the possibility that papillomaviruses (and other pathogens) may differentially interact with tonsils vs. other tissues. This is especially interesting in light of the increasing incidence of papillomavirus induced oropharyngeal (i.e., tonsil and base of tongue) cancer [133].

Origin-specific epithelial models to study host-pathogen interactions

These differences in tissue-specific gene expression may be important to consider when studying host-pathogen interactions, as they may influence the susceptibility and response of these tissues to infection. Epithelial tissues are a niche for pathogenic microbes, as well as commensals, including a broad range of epitheliotropic fungi, parasites, bacteria, and viruses. Many of these microorganisms have a very specific tropism for particular epithelia. For example, the dermatophyte fungi that cause ringworm thrive in the terminally differentiated apical layers of cornified epithelium and do not colonize host mucosa [134]. In contrast, other microorganisms, such as specific human papillomaviruses, exhibit tropism toward mucosal rather than cutaneous sites (and vice versa). Importantly, microbes know the difference between specific mucosal stratified epithelia. For example, Neisseria gonorrhoeae specifically infects the cervical epithelia, while a related pathogen (Neisseria meningitidis), is specific for the oral cavity [135]. Similarly, herpes simplex virus type 1 and type 2, primarily infect the genital and oral mucosa, respectively [136, 137]. Specific papillomavirus types may be even more fastidious, and evolved to infect specific niches on a host epithelia [138, 139]. What causes the exquisite tissue tropism seen in these infections is not always apparent, but these specializations do point to biologically critical differences between tissues. The data presented in this resource will allow researchers to generate testable hypotheses explaining these differences in tissue tropism.