All samples for this study were derived from an epidemiological survey of feline upper respiratory tract disease conducted in Wuhan between 2021 and 2024, encompassing a total of 1,427 cases (Supplementary material 1). From this dataset, we selected 8 FHV-1-positive samples based on stringent inclusion criteria, ensuring all specimens were obtained from adult British Shorthair cats exhibiting upper respiratory symptoms, housed in single-cat household, and having received routine vaccinations. The cohort maintained a balanced gender distribution (4 males and 4 females), with laboratory confirmation via RT-qPCR confirming FHV-1 positivity while excluding co-infections with feline calicivirus (FCV), Mycoplasma felis (Mf), Chlamydia felis (Cf), and Bordetella bronchiseptica (Bb). Additionally, four FHV-1-negative control samples were included, sourced from asymptomatic adult British Shorthairs under identical conditions (single-cat households, vaccinated, 2 males and 2 females), with RT-qPCR verification confirming the absence of FHV-1, FCV, Mf, Cf, and Bb. All samples were collected from animals that had not received antimicrobial treatment within the preceding week (Table 1).

The sterile swabs were inserted transnasally to the nasopharynx (depth adjusted by weight: 4–5 cm for <3 kg, 5–7 cm for ≥3 kg), rotated three times against the posterior pharyngeal wall, and immediately preserved in RNA/DNA stabilization solution at −80°C. All procedures complied with the Guide for the Care and Use of Laboratory Animals (National Research Council, 2011) and were approved by the Institutional Animal Care and Use Committee of HZAU (protocol code: HZAUCA-2025-0038). Owners provided written consent for sample collection and data use, including potential risks and benefits.

Total DNA was extracted using the TGuide S96 Magnetic Soil/Feces DNA Kit (Tiangen Biotech, China) with negative controls to exclude contamination (sterile ultrapure water) according to manufacturer’s instructions. The quality and quantity of the extracted DNA were examined using electrophoresis on a 1.8% agarose gel, and DNA concentration and purity were determined with NanoDrop 2000 UV–Vis spectrophotometer (Thermo Scientific, USA). The full-length 16S rRNA gene were amplified with primer pairs 27F: AGRGTTTGATYNTGGCTCAG and 1492R: TASGGHTACCTTGTTASGACTT. Both the forward and reverse 16S primers were tailed with sample-specific PacBio barcode sequences to allow for multiplexed sequencing. We chose to use barcoded primers because this reduces chimera formation as compared to the alternative protocol in which primers are added in a second PCR reaction. The KOD One PCR Master Mix (TOYOBOLife Science, Japan) was used to perform 25 cycles of PCR amplification, with initial denaturation at 95°C for 2 min, followed by 25 cycles of denaturation at 98°C for 10 s, annealing at 55°C for 30 s, and extension at 72°C for 1 min 30 s, and a final step at 72°C for 2 min. The total of PCR amplicons was purified with VAHTSTM DNA Clean Beads (Vazyme, China) and quantified using the Qubit dsDNA HS Assay Kit and Qubit 3.0 Fluorometer (Thermo Fisher Scientific, USA). After the individual quantification step, amplicons were pooled in equal amounts. SMRTbell libraries were prepared from the amplified DNA by SMRTbell Express Template Prep Kit 2.0 according to the manufacturer’s instructions (Pacific Biosciences, USA). Purified SMRTbell libraries from the pooled and barcoded samples were sequenced on a PacBio Sequel II platform (Beijing Biomarker Technologies, China) using Sequel II binding kit 2.0.

The bioinformatics analysis of this study was performed with the aid of the BMKCloud.¹ The raw reads generated from sequencing were filtered and demultiplexed using the SMRT Link software (version 8.0) with the minPasses ≥5 and minPredicted Accuracy ≥0.9, in order to obtain the circular consensus sequencing (CCS) reads. Subsequently, the lima (version 1.7.0) was employed to assign the CCS sequences to the corresponding samples based on their barcodes. CCS reads containing no primers and those reads beyond the length range (1,200–1,650 bp) were discarded through the recognition of forward and reverse primers and quality filtering using the Cutadapt (version 2.7) quality control process. The UCHIME algorithm (v8.1) was used in detecting and removing chimera sequences to obtain the clean reads. Sequences with similarity > 97% were clustered into the same operational taxonomic unit (OTU) by USEARCH (v10.0), and the OTUs < 2 in all samples were filtered.

Clean reads then were conducted on feature classification to output ASVs (amplicon sequence variants) by DADA2, and the ASVs < 2 in all samples were filtered. Taxonomy annotation of the OTUs/ASVs was performed based on the Naive Bayes classifier in QIIME2 using the SILVA database (release 138.1) with a confidence threshold of 70%. The Alpha diversity was calculated and displayed by the QIIME2 and R software, respectively. Beta diversity was determined to evaluate the degree of similarity of microbial communities from different samples using QIIME. The Permutational Multivariate Analysis of Variance (PERMANOVA) analyses were performed using the vegan package in R, with visualizations generated in Python. Prior to the above analysis, rarefaction curves and rank abundance curves were generated to assess the sequencing depth. Furthermore, differential taxa identification between groups was performed using Linear discriminant analysis effect size (LEfSe) analysis with a Linear Discriminant Analysis (LDA) score threshold of 2.0 and a p-value cutoff of < 0.05, combined with Wilcoxon rank-sum tests to validate group-specific differences in individual taxonomic units. For dominant genera with relatively high abundance (relative abundance > 1%), further validation was conducted using ANOVA tests for normally distributed data or Kruskal-Wallis tests for non-normally distributed data. p-values were corrected using the Benjamini-Hochberg method to control the false discovery rate.

Sequencing results showed that the FHV-1-positive group obtained an average of 51,679.375 effective CCS, while the FHV-1-negative group obtained an average of 54,363.25 effective CCS (Table 2). The rarefaction curve results indicated that all curves approached a saturation trend, suggesting that the number of species in this environment would not increase significantly with the increase in sequencing depth (Figures 1A,B). At 97% sequence similarity threshold, microbial community analysis identified 272 shared OTUs between FHV-1-positive group and FHV-1-negative group. FHV-1-positive group harbored 120 unique OTUs, whereas FHV-1-negative group contained 14 unique OTUs (Figures 1C,D). These findings demonstrate significant disparities in microbial diversity between the two groups, with FHV-1-positive group exhibiting remarkably higher microbial diversity than FHV-1-negative group.

Alpha diversity analysis demonstrated significant microbial community differences between groups (Figure 2). FHV-1-positive group exhibited significantly higher Shannon (p < 0.05) and Simpson indices (p < 0.01) compared to FHV-1-negative group, indicating enhanced microbial diversity in infected individuals. However, no significant intergroup differences were observed in Chao1 (p > 0.05) or ACE indices (p > 0.05), suggesting comparable species richness between groups (Table 3).

Principal component analysis (PCA) and Principal coordinates analysis (PCoA) revealed distinct microbial community structures between groups (Figures 3A,B). The FHV-1-negative group showed tight clustering with minimal intra-group variation, while the FHV-1-positive group exhibited significantly greater dispersion, indicating substantial compositional dissimilarities.

The PERMANOVA analysis results demonstrated differential sensitivity of distance metrics in detecting microbial community variation. Analysis using Bray–Curtis dissimilarity showed highly significant differences (R² = 0.285, p = 0.001) (Figure 3C), whereas the binary Jaccard metric yielded more modest but still statistically significant separation (R² = 0.238, p = 0.022) (Figure 3D). A larger R² indicates a higher degree of explanatory power of group partitioning for the observed differences, reflecting greater divergence between groups, while a p-value < 0.05 signifies high reliability of the test results.

The top 10 phyla in terms of species abundance are Bacteroidota, Firmicutes, Proteobacteria, Fusobacteriota, Spirochaetota, Patescibacteria, Desulfobacterota, Synergistota, Campylobacterota, and Actinobacteriota (Figures 4A,B). Among them, the abundance of Actinobacteriota in the FHV-1-positive group is significantly higher than that in the FHV-1-negative group (0.34% ± 0.05% vs. 0.17% ± 0.03%, p = 0.0433) (Table 4).

The top 10 genera in terms of species abundance are Porphyromonas, Bergeyella, Fusobacterium, Fusibacter, unclassified_Peptostreptococcaceae, unclassified_Chlorobi_bacterium_canine_oral_taxon_046, Capnocytophaga, Bacteroides, Treponema, and Frederiksenia (Figures 4C,D). At the genus level, Porphyromonas (16.77% ± 1.66% vs. 25.88% ± 0.86%, p = 0.0043), Bergeyella (4.68% ± 1.45% vs. 15.38% ± 0.76%, p = 0.0006) showed marked reductions in the FHV-1-positive group (Table 5). These reductions may disrupt commensal microbiota-mediated immune surveillance, potentially facilitating secondary infections.

LEfSe analysis revealed that the family Prevotellaceae (LDA score> 4, p < 0.05) was significantly enriched in the FHV-1-positive group compared to FHV-1-negative group (Figure 4E), with its taxonomic hierarchy further illustrated in the cladogram (Figure 4F).