Characterization of 3D organotypic epithelial tissues reveals tonsil-specific differences in tonic interferon signaling

Three-dimensional (3D) culturing techniques can recapitulate the stratified nature of multicellular epithelial tissues. Organotypic 3D epithelial tissue culture methods have several applications, including the study of tissue development and function, drug discovery and toxicity testing, host-pathogen interactions, and the development of tissue-engineered constructs for use in regenerative medicine. We grew 3D organotypic epithelial tissues from foreskin, cervix, and tonsil-derived primary cells and characterized the transcriptome of these in vitro tissue equivalents. Using the same 3D culturing method, all three tissues yielded stratified squamous epithelium, validated histologically using basal and superficial epithelial cell markers. The goal of this study was to use RNA-seq to compare gene expression patterns in these three types of epithelial tissues to gain a better understanding of the molecular mechanisms underlying their function and identify potential therapeutic targets for various diseases. Functional profiling by over-representation and gene set enrichment analysis revealed tissue-specific differences: i.e., cutaneous homeostasis and lipid metabolism in foreskin, extracellular matrix remodeling in cervix, and baseline innate immune differences in tonsil. Specifically, tonsillar epithelia may play an active role in shaping the immune microenvironment of the tonsil balancing inflammation and immune responses in the face of constant exposure to microbial insults. Overall, these data serve as a resource, with gene sets made available for the research community to explore, and as a foundation for understanding the epithelial heterogeneity and how it may impact their in vitro use. An online resource is available to investigate these data (https://viz.datascience.arizona.edu/3DEpiEx/).


Introduction
Epithelial tissues are found throughout the body and play important roles in various organ systems.The epithelium is multifunctional: providing basic external barrier protection from pathogens [1], regulation and secretion of bodily fluids [2], absorption of substances internally and externally [3], and sensation [4].Epithelia can be classified based on the number of cell layers and the shape of the cells [5].Tissues involved in absorption, filtration, or diffusion (e.g., the lining of the small intestine and the walls of blood vessels) consist of a single layer of cells [1].The more complex stratified epithelia are made up of multiple layers of cells and are found in locations where protection is the primary function.Stratified columnar epithelia have tall, column-shaped cells in the top layer.Stratified cuboidal epithelia have cube-shaped cells in the top layer.They are found in locations that need to withstand moderate mechanical stress, such as the lining of the bronchi in the respiratory system [6].Stratified squamous epithelia have flattened cells in the top layer.They are found in locations that are subjected to wear and tear, such as the skin, the oral cavity, and the female reproductive tract.The skin comprises multiple cell layers, including a layer of stratified squamous epithelia on the outer surface, which helps protect the body from physical trauma, ultraviolet radiation, and infection.The epidermal cells also help regulate body temperature and hydration.The oral mucosa protects the underlying tissues from the abrasive effects of chewing and the acidic environment created by certain foods and drinks [7], and in the defense against infection [8].The cervix is a crucial part of the female reproductive system, serving as the gateway between the uterus and the vagina.The cervix is lined with stratified squamous epithelium that helps protect the underlying tissues from mechanical and chemical stresses.The stratified squamous epithelium of the cervix also produces mucus that helps keep the cervix moist and helps to prevent infections [9].
Although skin, cervix, and oral epithelia are all classified as epithelial tissues, they have distinct embryological origins and developmental pathways.Skin is derived from the ectoderm [10], the outermost layer of the developing embryo.In contrast, the cervical epithelium is likely derived from the endoderm [11], the innermost layer of the developing embryo.In mice, however, conflicting evidence suggests a potential mesodermal origin [12][13][14].Oral epithelia, including the tonsils of the oropharynx, are derived from both the ectoderm and the endoderm [15,16].
The stratified epithelia of all these tissues are formed through gradual differentiation from stem-cell-like progenitor cells towards fully differentiated cells at the top of the epithelia [17].This stratified differentiation is not recapitulated using traditional 2D cell culture approaches.Therefore, differentiation-dependent gene expression cannot be studied using standard cell culture techniques.On the other hand, direct sequencing of biopsies derived from (human) tissues does not specifically isolate the epithelial cells, leading to contamination with other cell types present in the biopsy.
Organotypic 3D epithelial tissue culture methods are techniques used to grow and maintain epithelial cells in a 3D configuration that more closely resembles the in vivo tissue architecture and function [18][19][20][21][22][23][24].These in vitro methods for epithelial stratification were established in the 1970s and further refined over the following decades [25,26].They involve the use of specialized culture systems and media that support the growth and differentiation of cells in a 3D configuration and the formation of functional tissue-like structures.Specifically, we use a collagen gel matrix to support the growth of cells in 3D [27][28][29].Collagen is a structural protein found in many tissues and organs.It is commonly used in tissue culture because it is biocompatible and can support the growth and differentiation of many cell types [30].Primary, host biopsy-derived epithelial cells are cultured to isolate epithelial cells [31] and homogeneous populations of cells are grown on top of a collagen gel matrix mixed with human fibroblasts.
Cells are cultured at the air-liquid interface to promote differentiation.Organotypic 3D epithelial tissue culture methods have several applications, including the study of tissue development and function [32], drug discovery and toxicity testing [33], host-pathogen interactions [34][35][36][37][38][39], and the development of tissue-engineered constructs for use in regenerative medicine [40].
We chose a reductionist approach to carefully characterize and contrast the tissue-specific transcriptome of primary epithelial cells (foreskin, cervix, and tonsils) differentiated into stratified squamous epithelia.The use of 3D organotypic rafts in conjunction with RNA-seq represents a significant innovation in the study of epithelial tissues, as it enables a more comprehensive and accurate assessment of gene expression patterns [41,42].By characterizing multiple independent donors of distinct-origin epithelial cells, differentiated using an identical 3D raft culture method [43], we focus on these tissues' intrinsic biological similarities and differences.Using this approach, we compare gene expression patterns in these three types of epithelial tissues to better understand the molecular mechanisms underlying their function and identify potential therapeutic targets for various diseases.By characterizing these tissues, we provide a valuable resource for the research community in utilizing origin-specific epithelial tissue for in vitro experiments.
This study compared epithelial cells derived from foreskin, cervical, and tonsillar tissues.Despite their different origins, the epithelial cells within these tissues are involved in responding to infection [44].PAMPs, or pathogen-associated molecular patterns, are molecules associated with pathogens (such as viruses, bacteria, fungi, and parasites) and are recognized by pattern recognition receptors (PRRs) expressed by host cells.The recognition of PAMPs by PRRs triggers the activation of the immune response, including the production of type I interferons.Interferons are produced in response to viral infections or other types of cellular stress.Interferons play a vital role in the immune response by activating antiviral gene expression and inhibiting virus replication [45].Specifically, interferon binds to receptors on the surface of neighboring cells, activating a signaling pathway culminating in the expression of interferon-stimulated genes.
Interestingly, a subset of these interferon-stimulated genes is continuously expressed without an acute interferon signal [46].It is thought that these tonic interferon responses play a role in maintaining immune surveillance and protecting against viral infections.The tonsils play a vital role in the immune system, serving as the first line of defense against infections by trapping and neutralizing pathogens that enter the body through the oral cavity [47].Tonsillar epithelia are equipped with a variety of immune defense mechanisms.Some evidence suggests that the microenvironment in the tonsils may be different from that of other tissues [48][49][50].Tonsils and other lymphoid organs, such as the spleen and lymph nodes, are characterized by a high density of immune cells [51] and a complex network of blood vessels [52].By analyzing the extensive transcriptome data provided as part of this resource, we suggest that the epithelial cells play a critical role in shaping the immune microenvironment to balance tissue homeostasis and immune responses in the tonsils.

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 4 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 5 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 2+ 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 2 O rinsed, mounted with Pro-Long 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 flashfreezing 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.

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 preexisting 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.rproject.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/.

Results
Three-dimensional organotypic raft cultures are the current state-of-the-art method to reconstruct fully stratified 3D epithelial tissue in vitro [43,[73][74][75].Traditionally, most of these 3D tissues are grown from epithelial cells derived from foreskin keratinocytes [73,74,76,77].However, since epithelial tissues have different embryogenic origins, these foreskin-derived tissues may not be representative of other (squamous) epithelia.Therefore, we designed this comparative study to characterize 3D tissues derived from primary cells obtained from foreskin, cervical, or tonsillar biopsies.

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 foreskinderived 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].
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 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 skinderived 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 (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

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.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.
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).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 −5 ).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 tonsilderived 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 (Padj < 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).
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 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

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.
Cervix-derived epithelia overexpress genes involved in ECM remodeling and deplete keratinization compared to foreskin.We first compared foreskin tissue to cervical tissue.This comparison identified 605 differentially expressed genes: 409 up-and 196 down-regulated (Fig 6B).The top significantly up-regulated genes in cervix vs. foreskin were transcription factor FOXE1, putative transcription factor CRIP2, immune regulator THEMIS2, extracellular matrix proteoglycan VCAN, and metalloprotease ADAMTSL3.The top down-regulated genes were the calcium-binder CALML5, transcription factor HOXC10, water channel AQP3, suprabasal differentiation maker KRT10, and transmembrane protease TMPRSS11F.Gene set enrichment analysis on these 605 differentially expressed genes identified 86 biological

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 foreskinderived 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.
Foreskin-derived rafts uniquely up-regulate genes involved in keratinization.The foreskin-identifying epithelial markers consist of 237 differentially expressed genes plus 30 transcripts found in the Euler diagram's core representing differentially expressed genes in all pairwise comparisons.Most foreskin-specific differentially expressed genes are down-regulated in foreskin relative to cervix and tonsil (76.0%, 203/267, S14 File).One of the top down-regulated genes in foreskin was MUC5B, a gel-forming mucin responsible for lubricating oral and cervical mucosa [77].Network analysis of enriched functional terms identified a large interconnected functional module (S6 .These epidermis development genes include keratinization genes such as regulators of cornification (e.g., CALML5, [95]) and ceramide synthases (e.g., CERS3).Ceramides are the principal lipids in the stratum corneum and are an essential component of the skin barrier [96].

Tonsil-derived rafts uniquely down-regulate baseline antiviral immune responses.
When analyzing the oral-derived tonsillar epithelia, we identified 271 (242 + 29) tonsil-identifying differentially expressed genes (Fig 7B).116 of these 271 genes (42.8%) are down-regulated in tonsil-derived tissue (S17 File).Analyzing these genes, we identified 115 enriched functional terms involved in innate immune processes, specifically related to defense responses to viruses and type I interferon signaling, which are down-regulated in tonsillar tissues.These data further support earlier data (Fig 4 ) that the 3D organotypic tonsillar epithelia have a baseline innate immune suppression against viruses (S7 Fig and Fig 7E).
In parallel, an analysis of the 155 tonsil up-regulated genes (57.2%, 155/271, S18 File) identified 22 enriched terms (S19 File).Tonsil up-regulated enriched terms related to xenobiotic stimulus and detoxification (e.g., epoxygenase P450 pathway, phenol-containing compound metabolic process, and olefinic compound metabolic process; S7 Fig and Fig 7E) Cervix-derived rafts uniquely up-regulate a subset of Hox genes.The 118 cervix-specific epithelial markers include 101 differentially expressed genes from the shared overlap plus 17 core markers (Fig 7C).45% of these differentially expressed genes are down-regulated in cervix compared to the other tissues (44.9%, 53/118, S20 File).The 65 cervix-specific up-regulated differentially expressed genes (55.1%, 65/118, S21 File) were enriched for 115 terms (S22 File) related to embryogenic and morphogenic biological processes (S8 Fig and Fig 7F).Of note, several HOX cluster genes are up-regulated in cervical tissue.For example, expression of HOXA11 is critical for the development of the cervix and is strongly up-regulated in these tissues.The down-regulated genes do not appear to be enriched in any functional processes.

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 −18 and 3.17 x 10 −44 , relative to cervix and tonsil, respectively), followed by cervix (Wald P-adj = 4.91 x 10 −5 , 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 Padj = 2.9 x 10 −18 ) and tonsil (Wald P-adj = 1 x 10 −3 ), while AQP3 expression is higher in tonsils compared to cervical tissues (Wald P-adj = 9.1 x 10 −6 ).
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 −7 and 9.89 x 10 −20 , relative to foreskin and tonsil, respectively), followed by foreskin (Wald P-adj = 1.40 x 10 −3 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 Padj = 3.4 x 10 −12 and 3.8 x 10 −3 , relative to foreskin and tonsil, respectively), followed by tonsil (Wald P-adj = 2 x 10 −3 relative to foreskin) and foreskin derived tissues.
Compared to foreskin-derived tissues, FOXE1 (Forkhead Box E1) is up-regulated in both tonsil (Wald P-adj = 8.06 x 10 −111 ) and cervix (Wald P-adj = 9.04 x 10 −31 ).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 −25 ) 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 tissuespecific 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 tonsilderived signature positively identified 98.2% of oral mucosa samples (Esophagus and Salivary gland).Notably, all the other samples scored negative (Fig 8C).
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.
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.

Discussion
For this study, we used a reductionist approach, focusing on the epithelial cells comprising stratified squamous epithelial barriers.While modifications can be made to these 3D models depending on the experimental context and complexity desired, we used an organotypic 3D raft model with the same fibroblasts, collagen, media, growth conditions, and duration to maximize our ability to compare between tissues.For example, 3D organotypic skin and cervical models have been embedded with additional cell types, such as Langerhans cells [109,110], and others have tested tissue-specific hormone supplementation (estrogen and progesterone) when culturing 3D organotypic ectocervical epithelium [100].A primary goal of our study was to provide evidence that these different in vitro tissues maintain similarities with their in vivo tissue equivalents, independent of other cell types or microbiota normally present in these tissues.Our study used a strict statistical cutoff for identifying differentially expressed genes.Therefore, we may have missed minor (but biologically relevant) differences between the tissues.However, the identified transcripts are likely highly biologically relevant.Aggregating the pairwise contrast data yielded slightly more differentially expressed genes than the global likelihood ratio test: 1381 pairwise vs. 1238 global (see S23 File for the combined global and pairwise results).Comparing the pairwise vs. global differentially expressed genes, 1,194 (83.8%) were pairwise and global differentially expressed.44 (3.1%) were only global differentially expressed genes, and 187 (13.1%) were only pairwise differentially expressed genes.These differences between the sets of pairwise and global differentially expressed genes are due to the different statistical methods [59] underlying each test (global likelihood ratio test across all three sample groups vs. pairwise Wald tests between two sample groups) and the high stringency cutoffs (P-adj < 0.01) resulting in variable statistical power for detecting differentially expressed genes [111].

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][118][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.
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 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][127][128][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.

Conclusions
We characterized biological differences between foreskin, cervical, and tonsillar-derived 3D organotypic epithelial raft cultures using histological, immunofluorescent, and RNA-seq techniques.We confirmed that the identical 3D raft culturing technique yields distinct, origin-specific, stratified squamous epithelia.Epithelial origin is the primary differentiator between these differentiated cultures, yielding distinct transcriptional profiles in the form of thousands of differentially expressed genes.Interestingly, we demonstrate that the transcriptional profile of organotypic tonsillar epithelia is distinguished from foreskin and cervix by a down-regulated baseline expression of important (anti-viral) innate immune genes.Overall, these data are an important resource for those using 3D epithelial models in their research, epithelial biologists, and those studying host-pathogen relationships in these distinct epithelial sites.

Fig 1 .
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 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 surfacelike (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

Fig 4 .
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 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.

Fig 5 .
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

Fig 6 .
Fig 6.Pairwise contrasts identify the differentially expressed genes between foreskin, cervix, and tonsil-derived 3D epithelia.(A) Euler diagram with areaproportional 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 2 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 Fig) related to extracellular matrix organization (Fig 7D).Functional enrichment on the 64 up-regulated differentially expressed genes (S15 File) identified 4 terms (S16 File, S6 Fig) and is defined by keratinocyte differentiation genes (Fig 7D)

Fig 7 .
Fig 7. Tissue-identifying differentially expressed genes.(A-C) Scatterplots of pairwise contrast overlaps, axes represent log 2 fold change.(D-F) Overrepresentation 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 10 -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

Fig 8 .
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

Fig 9 .
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,

Table 1 . RNA-seq file quality summary.
Filename, Size (in bytes), and MD5 checksum were used for verifying file identity and integrity.FastQC was run on each sequencing read file to verify quality.The percentage of sequences with quality scores of at least 30 (�Q30%, 1 in 1000 probability of an incorrect base call, 99.9% inferred base call accuracy) was calculated from the raw "fastqc_data.txt"output for each file.https://doi.org/10.1371/journal.pone.0292368.t001

Table 2 . RNA-seq mapping and gene-level assignment summary.
Homo sapiens (GRCh38) mapped reads comprised 96.72 +/-0.20% of the total reads via HISAT2 (mapped + unmapped).FeatureCounts was used for gene-level counting of the HISAT2 alignment output through a series of category filters (left to right, starting with human gene "assigned", which is used for downstream differential expression analyses). https://doi.org/10.1371/journal.pone.0292368.t002