In this prospective case-control study, dogs with CR (n = 6) and healthy control dogs (n = 6) presented to Colorado State University Veterinary Teaching Hospital (CSU-VTH) for medical workup (CR) or routine dentistry (healthy controls) were screened for study enrolment. The protocols for this study were approved by the CSU Institutional Animal Care and Use Committee and the CSU-VTH Clinical Review Board (protocol #2858). Inclusion criteria for CR dogs included at least 1-month history of unilateral or bilateral nasal discharge, congestion, and/or sneezing. Histologic documentation of lymphoplasmacytic rhinitis with or without concurrent neutrophilic or eosinophilic rhinitis based on nasal mucosal biopsies, without evidence of neoplasia, fungal infection or other identified causes of CR, based on CT and rhinoscopy, was also required for study entry.

Exclusion criteria included any animals receiving immune suppressive medications, including glucocorticoids, leflunomide, mycophenolate, azathioprine, or cytotoxic chemotherapy in less than a month prior to the enrolment. Animals with previously diagnosed endocrine diseases (Cushing’s, diabetes mellitus) that might alter nasal immunity were also excluded. Dogs with pronounced epistaxis were also excluded. A diagnosis of other diseases that may cause the nasal signs (neoplasia, aspergillosis, foreign body, oronasal fistula, nasopharyngeal stenosis, chronic regurgitation) were also criteria for exclusion.

A control group of n = 6 dogs was selected from animals presented to the CSU-VTH for routine dental cleaning, with no reported nasal signs or other respiratory signs, gastrointestinal signs or dental diseases that may affect the nasal cavity (e.g., oronasal fistula) were included as the healthy control group. All animals enrolled had no antibiotic usage for 30 days prior to swab collection.

Nasal swabs were collected under general anesthesia by gently swabbing the distal third of the nasal cavity using 6-inch PurFlock Ultra, sterile flocked collection devices (Puritan Medical Products, Guilford, ME). Swabs (2 swabs obtained per each nostril) were placed in 15 mL conical tubes (Corning Inc. Corning, NY) containing RNA later (Thermo Fisher Scientific, Waltham, MA) and stored at 4C prior to RNA and DNA extraction. To extract RNA from nasal cells, swabs were vortexed for 1 min at high speed to dislodge cells, PBS was then added to RNA later at a ratio of 1:5. Cells were pelleted, then processed for RNA extraction using Qiagen RNeasy micro kit (Qiagen, Hilden, Germany) following manufactures instructions or DNA extraction as described below.

RNA concentrations were verified on Nanodrop 1000 Spectrophotometer (Thermo Fisher Scientific), and then sent to Novogene Corp. Inc. (Sacramento, CA) for RNA sequencing. RNA quality was determined using an Agilent 2100 Bioanalyzer system to generate RIN numbers (RNA integrity number), which ranged from 6.9 to 10 for all RNA samples submitted. At Novogene Corp., mRNA was purified using poly-T oligo-attached magnetic beads. After fragmentation, the first strand cDNA was synthesized using random hexamer primers followed by the second strand cDNA synthesis. The library was completed following end repair, A-tailing, adapter ligation, size selection, amplification, and purification. Quantified libraries were pooled and sequenced on an Illumina NovaSeq 6000 (Illumina, San Diego, CA). 150 bp paired end reads were generated, and files were delivered as de-multiplexed fastq files.

Sequence data were analyzed on Partek Flow software, version 10.0 (Partek Inc. Chesterfield, MO). Raw data were filtered by removing reads containing adapters and reads containing N > 10% and for Phred scores >30. Filtered reads were aligned with STAR 2.7.3a, using CanFam3.1 genome assembly. Aligned reads were annotated and counted using HT-seq (31) with Ensembl 107, and differentially expressed genes were identified using DEseq2 (32). (Differential gene expression analysis based on negative binomial distribution). Biological interpretations included gene ontology and gene set enrichment analysis (GSEA).¹ Gene sets Hallmarks v2022.1, biocarta v2022.1, KEGG v2022.1, Gene Ontology go.bp v2022.1, and ImmuneSigDB v2022.1 were used for comparisons. Significant pathways were filtered using false discovery rate (FDR) q-value of≤0.25 and NOM p-value ≤0.05. Single cell gene expression data from the Human Protein Atlas² was used to identify olfactory and epithelial cell specific gene sets that were not expressed by immune cells.

Previously published single cell sequencing data from human nasal epithelium (33) was used to identify unique genes in nasal cell types (secretory cells, goblet cells, club cells, ionocytes, ciliary cells and rapidly dividing cells RDC) that were not found in previously published RNAseq data fromcanine peripheral blood mononuclear cells (34).

Additional swabs were collected for microbial analysis. Swabs were cut and immersed in extraction tubes following Qiagen DNeasy PowerSoil Pro Kit instructions (Qiagen, Hilden, Germany). Microbial DNA was frozen at −80°C and sent to Anschutz Center for Microbiome Excellence (ACME), at University of Colorado Anschutz Medical Campus, Aurora, CO for microbial sequencing. The library was prepared according to Earth Microbiome project protocol,³ with 35 PCR cycles, using 515F and 806R primers. Samples were run on MiniSeq cartridges on Illumina Miseq sequencing instruments. Microbial sequence analyses were performed with QIIME2 (35). Microbial community similarity was displayed with principal coordinate analysis (PCoA) plots. Alpha diversity was determined using Shannon, Faith, and pielou indices. Beta diversity using weighted and unweighted UniFrac, as well as Bray Curtis. Alpha diversity indices were compared using a paired t-test, and beta diversity metrics were compared with PERMANOVA. Analysis of composition of microbiomes (ANCOM) was employed to determine the sequence variants that differed significantly between treatment groups (36). In addition, Linear discriminant analysis Effect Size (LEfSE) was also used to calculate the taxa that best discriminated between rhinitis or healthy group⁴ (37). Because of the low overall bacterial abundance in the swab samples, additional PCR cycles were needed to reach the required number of reads for each sample. Microbial features were filtered for a minimum frequency of 50 were removed and features not present in >2 samples were also removed, resulting in a total of 106 total features. Median frequency for n = 12 samples was 14,455.5 and total frequency was 228,553.

Single cell RNA sequence data from n = 17 dogs (38) was used to generate a signature matrix file for Cibersort analysis of immune cells in the nasal swabs (39). Using the data from 74,067 immune cells in blood that were obtained previously from scRNAseq analysis from 17 dogs, 41 transcriptionally distinct leukocyte populations were identified. For analysis of leukocyte populations in dog nasal samples, each of the 41 cell populations were downsampled to obtain equal representation across each cell type, then normalized count data for each population was summed up to obtain a cell type by gene matrix. After collapsing count data by cell type, the dataset was filtered to only include the most highly variable features that define each cell population. Filtering was completed using the output of the FindAllMarkers function, with features being selected if they had an average log2 fold change (relative to all other cells) greater than 0.85 and a percent expression greater than 75% within the cell population. After initial filtering, any overlapping features were then removed from the dataset. This data set was then uploaded to⁵ along with the median ratio normalized counts from the current study to determine immune cell type enrichment (similar datasets were not available for epithelial cell populations in dogs).

For this analysis, significant differentially expressed genes (DEGs) were extrapolated from the RNA sequencing data. These median ratio normalized reads from individual samples were then correlated to the percent relative abundance of 160 unique bacterial taxon found in n = 12 samples (healthy and rhinitis) using rcorr (Rstudio) (40). p values for significance and r values for correlation were generated for each gene to taxon pair. Protein coding genes with correlations p-values ≤0.05 were then entered into string protein database⁶ (41) for categorizing the protein sets.

The RNA sequencing data including fastq files and read counts presented in this manuscript have been deposited in NCBI’s Gene Expression Omnibus and are accessible through GEO Series accession number GSE255529. All other data is available from the corresponding author on reasonable request.

To address the question of how the nasal transcriptome and microbiomes may be altered in dogs with CR, we designed a prospective study to sequence deep nasal swab samples from n = 6 healthy dogs and n = 6 dogs with CR of at least 1 month duration. The demographic data of the CR dogs and heathy dogs are listed in Tables 1, 2. This study was approximately matched for age and gender, though the healthy dogs were overall younger than CR dogs. Within the CR group, the severity of clinical signs ranged from mild to erosive and more severe CR (Table 1).

The first question addressed in this study was whether the transcriptomes of nasal mucosal cells, including epithelial cells, immune cells, and stromal cells, were significantly different in CR dogs compared to healthy. A total of 17,975 protein coding genes were identified in the 12 analyzed samples. The samples from the n = 6 CR dogs showed distinct clustering compared to the n = 6 healthy dogs (Figure 1A) with a high degree of transcriptome heterogeneity present in the CR samples. By contrast (Figure 1B) the heathy nasal transcriptomes displayed a relatively high level of within-sample homogeneity. Next, differential gene expression (DEseq) analysis was done comparing CR transcriptomes to healthy transcriptomes (based on analysis of 14,504 protein coding genes with known annotation). We found 234 significantly upregulated genes with a fold-change (FC) ≥ 2 Log₂ and FDR adjusted p-value of ≤ 0.05. The analysis also revealed 385 genes whose expression was significantly downregulated, with FC ≤ −2 Log₂ and FDR ≤ 0.05 (Figure 1C). The most upregulated genes in the CR samples included genes encoding OXER1, OSM, PROK2, MMP8 and TEX48, while the genes with the most downregulated expression included OR7G6, CDKL4, SNTN, DZIP1L and LMNTD1 (Figure 1D). Pathway analysis of the differentially expressed immune-related genes (DEG) revealed upregulation of inflammatory pathways such as the TNF-α and IFN-γ and IFN-α pathways in the CR animals, along with significant downregulation of cell cycle (CTCF) and activin receptor kinase (ALK) and TGF-β pathways in CR dogs (Figure 2A).

Epithelial cell dysfunction may also be an important contributor to the pathogenesis of CR, so we also analyzed the epithelial cell transcriptomes in study dogs. To conduct this analysis, we interrogated human single cell databases (33) for epithelial cell specific genes, including those associated with mucous and serous glandular cells, salivary duct cells, ciliated cells, squamous epithelial cells. Our analysis also excluded genes expressed by canine leukocytes (34). From this analysis, we found that there were numerous genes (117 in total) associated with ciliated airway cell function whose expression was significantly downregulated in dogs with CR (Supplementary Figures S2A,B). In addition, genes associated with cilia and flagella associated proteins, mucin and olfactory receptors were also downregulated (Supplementary Figure S2C), along with genes expressed by other nasal epithelium specific cell types (Figure 2B). Pathways related to cilium assembly, movement and regulation were also found to be downregulated in samples from dogs with CR (Figure 2C). These findings suggest that ciliary dysfunction and defective clearance of airway mucus and debris may play an important role in CR pathogenesis in dogs. However, the current study was not able to determine whether ciliary dysfunction was a primary or secondary occurrence in dogs with CR.

In addition to bulk transcriptome analysis, a Cibersort matrix generated from canine single cell leukocyte sequencing data (38) was used to deconvolute the bulk RNAseq data from the nasal swabs in the present study into gene sets derived from specific immune cell populations. Using Cibersort analysis, we identified a total of 23 different immune cell populations in the 12 RNAseq samples from healthy dogs and dogs with CR. We found that the percentage of polymorphonuclear-myeloid derived suppressor cells (PMN-MDSC) was significantly higher in nasal swabs from dogs with CR compared to healthy dogs, while unexpectedly the percentage of eosinophils was significantly lower in dogs with CR compared to healthy dogs (Figure 2D).

Microbial 16S rRNA sequencing revealed 160 unique taxa in the nasal swab samples. To assess the overall diversity in bacterial populations between groups, the alpha and beta diversity were compared between CR and healthy samples, which revealed that there were no significant differences in alpha diversity, either using faith or Shannon metrics (Figures 3A,B). However, beta diversity was significantly different between CR and healthy samples, with distinctly separated clusters of the CR and healthy samples for Bray Curtis beta diversity and for weighed UniFrac distance measurement (Figures 3C,D). These analyses therefore indicated the presence of a dysbiotic bacterial flora in the nose of dogs with CR, which has been reported previously (24).

Of the 11 phyla identified in the nose of dogs with CR and in healthy dogs, the most abundant phylum was Proteobacteria in the healthy dogs and Firmicutes in dogs with CR (Figure 4A). Four phyla were found to be significantly different in terms of relative abundances between dogs with CR and healthy dogs, as assessed by ANCOM analysis (Figure 4C). These included Proteobacteria and GN02 which had a significantly higher abundance in dogs with CR, whereas Firmicutes and Actinobacteria were significantly less abundant in dogs with CR compared to healthy dogs (Figure 4C).

Sixteen class level taxon were identified, with the most abundant class being Gammaproteobacteria in healthy dogs and Bacillus in dogs with CR (Figure 4B). Significant differences were found by ANCOM analysis between CR and healthy dogs for both of these most abundant bacterial classes (Figure 4D). Linear discriminant analysis Effect Size (LEfSe) was also used to interrogate the differences in bacterial populations, and this analysis revealed additional family and genus level differences which significantly contributed to the altered mucosal bacterial populations present in CR samples. For example, Arthrobacter and Kingella family Chitinophagaceae had higher relative abundance in the CR group, while the family Mycoplasmataceae had a greater relative abundance in the healthy control group (Figure 4E).

Finally, we used the nasal transcriptome and microbiome data to develop interaction mapping to better understand potential interactions between host nasal responses and microbial populations. First, the DEGs (see Figure 1) were selected from the transcriptome data, and the normalized expression values for 619 differentially expressed genes were then correlated to the relative abundance percentage of 160 bacterial taxa. This analysis revealed that 502 host genes in the nose exhibited a significant correlation (p ≤ 0.05) with at least a single bacterial taxon.

Five bacterial genera, including Enterococcus, Streptococcus, Helcococcus, Porphyromonas, and Arthrobacter, exhibited significant correlations with 15 differentially expressed genes (Figure 5A). These potentially interacting genes were then entered into the String protein database to find meaningful connections or interactions. Protein network association showed enrichment for host gene sets amongst 3 of the 5 bacterial genera. Interestingly, the gene sets that correlated with Porphyromonas and Arthrobacter mapped to multiple immune pathways (Figures 5B,C). In contrast, the gene sets correlated with Streptococcus involved primarily motility and cytoskeletal structure related genes, including many of the ciliary pathways found to be downregulated in CR dogs (Figure 5D). Thus, our interactome analysis revealed several important host and bacterial interactions, involving both immune and nasal epithelial cell responses and particular populations of bacteria, including some (e.g., Porphyromonas) previously identified as pathobionts.