All study procedures conformed to guidelines set by the US National Institutes of Health (Bethesda, MD, USA). The study protocol was approved (2019271A) by the Animal Ethics Committee of West China Hospital within Sichuan University (Chengdu, China).

Details of all regents and antibodies utilized in this study are listed in Table S1 .

Sprague–Dawley rats were housed in a room at constant temperature (25°C) and relative humidity (55%) and exposed to a 12-h light–dark cycle in a specific pathogen-free facility at the Laboratory Animal Center in Chengdu Dashuo Biotechnology (Chengdu, China). All rats had free access to water and chow.

Ten (i.e., 20 ME bullae) male rats (Rattus norvegicus; 8–10 weeks) were divided randomly into the AOM group and control (untreated) group. Otoscopic examinations were undertaken before treatment for all rats to ensure that tympanic membranes were normal and ME effusion was absent. Subsequently, rats in the AOM group were anesthetized with 2% pentobarbital sodium (0.2 mL/100 g, i.p.). Then, lipopolysaccharide (LPS; Sigma–Aldrich, Saint Louis, MO, USA) was injected into the bilateral ME cavity through the tympanic membrane. About 72 h after LPS injection in the AOM group, three rats from each group were sacrificed. In order to reach the requirement of cell number constructing scRNA-seq libraries, the corresponding six ME bullae in each group were isolated carefully and mixed for preparation of single-cell suspensions, respectively. The remaining rats in each group were sacrificed and their ME bullae isolated for hematoxylin and eosin (H&E) staining and immunofluorescence staining.

The MEM was harvested, as described previously (9). The harvested MEM was dissociated with collagenase type I (1 mg/mL; catalog number, SCR103; Sigma–Aldrich) and Dispase^® II (2 mg/mL; neutral protease, grade II; 4942078001; Sigma–Aldrich) in Dulbecco’s phosphate-buffered saline containing Ca²⁺ and Mg²⁺ for 25 min in a 37°C shaking incubator at 100 rpm. The MEM was dissociated further in 1.5 mL of Accumax™ Cell Dissociation Solution (AM105; Innovative Cell Technologies, San Diego, CA, USA) for 8 min in a 37°C shaking incubator at 100 rpm to aid isolation of single cells. Dissociated cells were filtered with 70-μm cell strainers (BD Biosciences, Franklin Lakes, NJ, USA) to eliminate clumps, and incubated in 1 mL of Red Blood Cell Lysing Solution (BD Biosciences) to remove red blood cells. Finally, we utilized 35-μm cell strainers (BD Biosciences) to filter cells and collected them according to manufacturer recommendations. Cell viability >80% was required for subsequent construction of libraries.

Based on the HiSeq X™ platform (Illumina, San Diego, CA, USA), cells were labeled with sample tags from the Human Single-Cell Multiplexing Kit from BD Biosciences. Then, they were counted, multiplexed, and prepared for subsequent single-cell capture. The latter and complementary-DNA synthesis were undertaken by the Rhapsody Single-Cell Analysis System according to manufacturer (BD Biosciences) recommendations. Libraries were sequenced on multiple runs of NextSeq™ (Illumina).

Raw scRNA-seq data were processed according to a bioinformatics pipeline (BD Rhapsody™ WTA 1.9.1; https://bitbucket.org/CRSwDev/cwl/src/master/). This process included demultiplexing FASTQ reads, mapping to the rat genome (Rnor_6.0, STAR v2.5.2b), and generating gene/barcode matrices.

R package Seurat v4.0 was used for downstream data analyses (21). Following a standard workflow, we filtered cells that had a unique feature count >4000 or <200 and cells with a mitochondrial count >20% in the inflamed-MEM group. We also filtered cells that had a unique feature counts >3000 or <200 and cells with a mitochondrial count >40% in the normal-MEM group. Then, the default parameters of the “NormalizeData” function of Seurat were used to normalize the feature-expression measurements for each cell by the total expression. Finally, 7024 cells of the case group and 5258 cells of the control group were introduced into a combined Seurat object via “FindIntegrationAnchors” and “IntegrateData” functions. Then, variable genes were carried forward into scaling and principal component (PC) analysis. Significant PCs (top 30) were used for t-SNE analysis and clustering through “RunTSNE” and “FindClusters” functions (resolution = 0.6). To identify cell types, we checked whether the well-studied marker genes were the top differentially expressed genes (DEGs), and annotated the most likely identity for each cell cluster. The remaining cell types were identified by manual retrieval from a database.

Differentially expressed genes (DEGs) among clusters were detected by the Seurat function “FindAllMarkers”. Selected functional DEGs of each cluster were visualized via stacked violin plots. Volcano plots were applied to show the genes with upregulated or downregulated expression. Heatmaps showing the expression distribution of marker genes of cell clusters were created by pheatmap 1.0.12 (R package). Genes with Benjamini–Hochberg-adjusted P < 0.05 and absolute log₂fold change between two groups >2 were used for analysis of fuctional enrichment using the Gene Ontology (GO) database (clusterProfiler 3.16.1). In addition, ranked gene set enrichment analysis (GSEA) was undertaken: genes were ranked based on their phenotypes, and the GSEA algorithm proposed by Subramanian and coworkers (Subramanian et al., 2005) was used to calculate the enrichment score of each gene.

We used Monocle 2 v2.5.4 (R package) (22) to order single cells in “pseudotime”, and placed them along an inferred trajectory. After passing quality control, genes were placed into the Reversed Graph Embedding algorithm of Monocle to shape the trajectory. Then, Monocle applied a dimensionality reduction to the data and ordered the cells in pseudotime.

Cell-to-cell communication (CellChat 0.0.2; R package) was ascertained by evaluating expression of pairs of ligands and receptors within cell populations (23, 24). We examined the interaction between different cell types, and gene expression of 0.2 was set as the valid cutoff point.

The CS denotes the mean expression of specific gene functions (e.g., scores for the cell cycle, M1 macrophages, M2 macrophages, epithelia) in single cells. The CS was calculated, as described previously (25, 26). Briefly, a control gene set (Gc) was defined and genes were partitioned into 25 “bins” according to their mean expression across all cells. Then, for each gene from the target gene set (Gt), 100 genes were selected randomly from the bin to which this gene belonged to constitute the Gc. Therefore, the CS was calculated as mean (Gt) – mean (Gc). Relative expression was used to calculate the CS.

Two representative bullae from AOM group and control group were carefully isolated, and then fixed and embedded them with formalin and paraffin, respectively. The formalin fixed paraffin embedded bullae were sectioned and processed by using Opal Polaris™ 7-color Manual IHC Kit (Akoya Biosciences) as the manufacturer’s recommendation (27). Specifically, the opal panel included DAPI (Abcam), anti-CD68 (Abcam), anti-EPCAM (Abcam), anti-MKI67 (Abcam), anti-COL1A1 (Abcam) and anti-SFRP4 (Abcam).

To systematically investigate the cellular diversity of the inflamed MEM in rats as well as the sophisticated intercellular crosstalk within its multicellular ecosystem. A model of AOM in rats was established ( Figure 1A ), which was examined by histology using H&E staining ( Figures 1B, C ). Subsequently, we undertook scRNA-seq on 7023 cells from six samples of the inflamed MEM and 5258 cells from six healthy samples via BD WTA Rhapsody Analysis Pipeline 1.9.1 after quality control and removal of doublets ( Figure 1A and Table S2 ).

Then, sequencing data from the 12 samples were integrated and analyzed using Seurat. Overall, the transcriptomes of eight major cell types were captured based on expression of known canonical gene markers ( Table S2 ). Unsupervised clustering analysis categorized all cells further into 20 distinct clusters ( Figure 1D and Table S3 ). Each of the 20 clusters possessed specifically expressed marker genes, which were consistent to their distinct cell identities, as described previously ( Figure 1E , Table S4 and Supplementary Information 1 ), including neutrophils (C0, C2, C3 and C9), macrophages (C1 and C13), fibroblasts (C4, C6 and C8), epithelial cells (C5, C7, C11, C18 and C19), T cells (C10), dendritic cells (C12), endothelial cells (C14, C15 and C16) and B cells (C17).

Overall, scRNA-seq of the MEM of rats revealed considerable cell diversity that was underestimated previously. Detailed marker gene lists were provided in Table S2 .

To demonstrate the conversion of cell types in the MEM of rats in response to proinflammatory stimuli. We compared the cell compositions between the normal MEM and inflamed MEM, of which obvious differences were identified. Epithelial cells (C0, C3, C12, C14, C15 and C17) and fibroblasts (C1, C2 and C4) were the major cell types of the normal MEM, which accounted for approximately two-thirds of all single cells. However, after exposure to LPS, rapid infiltration of innate immune cells, including neutrophils (C0, C1, C2 and C4) and macrophages (C3, C5 and C6) occurred, and overwhelmed the normally predominant cell types ( Figures 2A, B ). More specifically, compared with their counterparts from the normal MEM, an obvious increase in the number of neutrophils and macrophages as well as an obvious reduction in the number of fibroblasts and epithelial cells, were observed in inflamed MEM samples ( Figure 2C ).

Taken together, these data suggested that, after exposure to LPS, the MEM in rats underwent significant conversion of cell types, which was characterized by the rapid infiltration of macrophages and neutrophils.

Expression of various ligand–receptor pairs was investigated to decipher the cell–cell interactions in the multicellular system of the normal MEM and inflamed MEM. In the inflamed MEM, the macro_2 cluster (C1), which specifically expressed Mrc1 (canonical marker gene of M2 macrophages) and was polarized towards the M2 phenotype ( Table S4 ), displayed vast communication with nearly all other cell clusters (especially those mediated by other macrophage clusters and neutrophils) and seemed to play a key part in inflammatory intercellular crosstalk ( Figures 3A, B ). We further investigated the specific ligand–receptor interactions among different cell clusters. In the inflamed MEM, we identified strong somatotrophic interactions among macrophages as well as between macrophages (especially macro_2) and other clusters, including colony-stimulating factor 1–colony stimulating factor 1 receptor (CSF1–CSF1R) and interleukin 34 (IL34)–CSF1R, which might mediate the maintenance and proliferation of macrophages after exposure to external stimuli (28). Besides, recruitment of neutrophils, natural killer cells, and T cells by macrophages (especially macro_2) via chemoattraction was observed (e.g., C-C motif chemokine ligand 3–C-C motif chemokine receptor 5 (CCL3–CCR5) and C-X-C motif chemokine ligand 1–C-X-C motif chemokine receptor 2 (CXCL1–CXCR2)). Specifically, interleukin 1 beta–interleukin 1 receptor 2 (IL1β–IL1R2) interactions between macrophages and neutrophils (which have been shown to be of vital importance in regulating the proinflammatory activities of IL-1 and contributing to the resolution of acute inflammation) (29) were detected ( Figure 3C ).

In the normal MEM, CSF1–CSF1R and IL34–CSF1R interactions were enhanced significantly between epithelial cells and macrophages (especially macro_2) as well as fibroblasts and macrophages (especially macro_2). These data suggested that epithelial cells and fibroblasts might be the main source of CSF1/IL34 under normal conditions, which has a crucial role in macrophage survival (28). Besides, inhibition of expression of proinflammatory phenotypes of potential non-professional phagocytes, including epithelial cells, fibroblasts, and endothelial cells, by macrophages (especially macro_2) via insulin like growth factor 1–insulin like growth factor 1 receptor (IGF1–IGF1R) interactions, was also observed, which might help to maintain the microenvironment of the normal MEM ( Figure 3D ) (30). Overall, these data suggested that M2 macrophages might play an important part in the MEM under normal conditions and inflammatory conditions.

Ligand–receptor interactions within certain signaling pathways were also predicted. Specific signaling pathways, including transforming growth factor beta (TGFβ), CCL, IGF and IL1, were identified. They displayed highly distinct cell-communication networks between the normal MEM and inflamed MEM, and might play an important part in the development and progression of AOM ( Figure 3E ) (20, 31–37). More detailed cell–cell interactions are presented in Figure S2 .

Overall, M2 macrophages seemed to have a key role in inflammatory intercellular crosstalk, which facilitated the maintenance and proliferation of macrophages, cell chemotaxis, as well as regulation of the proinflammatory activities of cytokines.

Studies have revealed the coexistence of professional phagocytes and non-professional phagocytes in multiple tissues, between which crosstalk disruption has been shown to be responsible for various diseases and might play an important part in the inflammatory response (38, 39). That hypothesis was bolstered in the present study. Specifically, besides professional phagocytes (neutrophils, macrophages and dendritic cells), other clusters from the inflamed MEM and normal MEM also displayed high expression of phagocytosis-related genes, such as Cd14 and Cd68 ( Tables S4 and S6 ), indicating that these cell clusters might have dual features. Therefore, to further confirm the existence of phagocytosis-related dual-feature cells in the MEM, we specifically extracted all CD68 ⁺ single cells for subsequent analyses. Doublets were identified and excluded by using “DoubletFinder” 2.0.3 (R package). Clustering tree 0.4.3 revealed that a stable clustering strategy was achieved at a resolution of 0.6, which generated 12 clusters belonging to seven major cell types ( Figure S3 , Figure 4A , Table S6–4 and Supplementary Information 2 ).

Apart from well-known professional phagocytes such as macrophages (C0, C1, C5 and C10), monocytes (C2), dendritic cells (C4, C6 and C11), and granulocytes (C3), three rare cell clusters (C7–C9) with dual features were also identified. Cluster C7 displayed high expression of cell cycle-related genes, including Ccna2, Ccnb1 and Cdk1, as well as stem-like genes Mki67 (marker gene of proliferation) and Aspm (involved in regulation of the mitotic spindle and coordination of mitotic processes), so cluster C7 comprised Mki67 ⁺ progenitor cells ( Figure 4B and Tables S6–4 ) (40). Cluster C8 specifically expressed Epcam (a gene associated with epithelial cells) and Krt18 (cytokeratin gene), which strongly indicated that it was an epithelial cell ( Figure 4B and Table S6–4 ). Cluster C9 showed high expression of genes related to the development and maintenance of bone, including Bmp5 and Cdh11, as well as genes related to matrix homeostasis, including Col5a2, Col12a1 and Mmp2, which is consistent with osteoblastic stromal cells ( Figure 4B and Tables S6–4 ). Additionally, the majority of CD68⁺ cells derived from inflamed MEM and the cell compositions as well as cell proportions showed obvious differences between normal MEM and inflamed MEM ( Figures 4C, D ).

To further quantify the dual features of these non-professional phagocytes, we calculated the CS of target gene sets of selected cell clusters (Mki67 ⁺ progenitor cells, epithelial cells, and osteoblastic stromal cells; expression of their selected marker genes is shown in Figure S4A ). Some Mki67 ⁺ progenitor cells had a high phagocytosis score and proliferation score ( Figure 4E ), which further demonstrated their phagocytic-proliferating dual-feature ( Tables S8-1 , S8–2 ). To a certain extent, some epithelial cells and osteoblastic stromal cells also showed trends of a dual feature related to phagocytosis ( Figures 4F, G and Tables S8–1 , S8–3 , S8–4 ). The M1 score and M2 score were calculated to clarify the polarization status of the four subpopulations of macrophages. M2 was the main phenotype of these four clusters, whereas the gene signatures of M1 and M2 were not mutually exclusive. For example, many C1qc ⁺ macrophages also presented higher expression of M1 marker genes than the other three clusters ( Figures S4B, C ).

To further confirm the existence of phagocytosis-related dual-feature cells, multiplex immunofluorescence was conducted in both normal MEM ( Figures 5A, B, E ) and inflamed MEM ( Figures 5C–E ). And the results showed that certain cells in both normal MEM and inflamed MEM co-expressed phagocytosis marker CD68 and proliferation marker MKI67, consistent with the definition of phagocytic-proliferating dual-feature cells ( Figure 5F ). Besides, the colocalization of CD68 and epithelial marker EPCAM were observed in both normal and inflamed MEM, suggesting of cells with phagocytic-epithelial dual feature ( Figure 5F ). Similarly, the existence of phagocytic-fibroblastic dual-feature cells were also confirmed in both normal and inflamed MEM based on the co-expression of CD68 and fibroblast marker SFRP4 or COL1A1 ( Figure 5F ).

Taken together, we identified three rare cell clusters with a phagocytosis-related dual feature. These cells are consistent with the definition of previously described non-professional phagocytes, and coexist with professional phagocytes in the MEM of rats.

To investigate the roles these dual-feature cells might have in inflammation progression, analyses of enriched genes expressed specifically by each of the three clusters was undertaken using the GO database. Mki67 ⁺ progenitor cells were characterized by enrichment of the pathways associated with “cell cycle”, “production and response to cytokines”, and “phagocytosis” ( Figure 6A ). Epithelial cells showed enrichment in the signaling pathways associated with “endocytosis”, “cell chemotaxis”, “cell differentiation”, and “response to inflammatory cytokines” ( Figure 6B ). Osteoblastic stromal cells showed significant enrichment in the signaling pathways associated with “bone development”, “endocytosis”, “cell chemotaxis” and, to a lesser extent, “negative regulation of inflammatory response” ( Figure 6C ). Hence, all of these dual-feature cells showed enrichment in signaling pathways associated with immune regulation and phagocytosis-related functions.

We ordered these dual-feature cells along a pseudotime trajectory to aid understanding of which developmental stages these cells might have in their immunoregulatory functions. This ordering resulted in a bifurcated early-to-late trajectory which represented the dynamic changes of the normal MEM after exposure to LPS ( Figure 6D ). These results showed that most of the epithelial cells and osteoblastic stromal cells were derived from the normal MEM, whereas Mki67 ⁺ progenitor cells were scattered along the whole trajectory ( Figure 6E ).

Taken together, it is postulated that epithelial cells and osteoblastic stromal cells exist mainly in the normal MEM, which might help to maintain the homeostasis of the immune microenvironment through phagocytosis-related function and negative regulation of the inflammatory response. However, Mki67 ⁺ progenitor cells might be activated and proliferate in response to proinflammatory stimuli and undertake their immunoregulatory functions throughout the whole process of inflammatory progression.