Single-cell transcriptomic characterization of mouse skin response to West Nile virus- infected mosquito bite reveals fibroblast-mediated barriers to transmission

Experiment: Single-cell RNA sequencing was applied to skin exposed to infectious mosquito bites in C57bl6/J immuno-competent mouse model to characterize the early molecular and cellular responses to viral transmission: 2-mm skin biopsies corresponding to the areas bitten by WNV-infected mosquitoes were collected at 6 h post biting. Each mouse was bitten by three to five mosquitoes (already placed in 3D-printed cages (3-D mosquito cage plan)). Twenty-two biopsies collected from eight animals were pooled (2 replicates, inf1 and inf2). The same number of skin biopsies were collected from unbitten mice that underwent the same procedures as control (rest skin). Skin biopsies were cell dissociated and dead cells were removed. Live cells were counted using and concentrated to ~1,000 cell / µl. GEM generation, barcoding, post GEM-RT cleanup, cDNA amplification, and cDNA library construction were performed using the Chromium Single-Cell 5′ Reagent version 2 kit (10x Genomics). Libraries were sequenced on an Illumina HiSeq4000 platform (Illumina).

Data processing and quality control for 5'end transcripts: Raw counts were aligned against the Mus musculus (mm39), Ae. aegypti (AaegL5.0), and West Nile Virus genomes using 10X Genomics Cell Ranger v7.1.0. Data was preprocessed and integrated with Seurat. First, cells were removed if total RNA count was low (nCount_RNA < 500), number of expressed features was too low or too high (nFeature_RNA < 500 and nFeature_RNA > 6000), and mitochondrial gene count was high (percent.mt > 10). Second, the samples were normalized, clustered, and projected onto low-dimensional embedding (UMAP) following the standard Seurat workflow with default parameters. The first 20 PCA dimensions were used for computation of neighbors and low-dimensional embedding. Third, the samples (1 rest (resting skin), inf 1 and inf2 (infectious bite) were integrated by selecting integration features (SelectIntegrationFeatures()) and anchors (FindIntegrationAnchors()), resulting in the integrated dataset (IntegrateData()). Finally, the integrated data was rescaled and UMAP representation was computed again. The two infectious bite samples were merged into a single sample for downstream analysis.

Clustering, annotation, and proportion testing of 5'end transcripts: Cell types were identified by computing cell neighbors using the top 15 PCA dimensions, followed by clustering with resolution = 0.02, resulting in n = 8 clusters. This resolution threshold was identified by trial-and-error and based on the clear separation between the cell types in the UMAP space. Clusters were annotated based on the expression of marker genes identified with the FindAllMarkers() Seurat function with parameters (only.pos = TRUE, min.pct = 0.8, logfc.threshold = 0.5) based on Wilcoxon rank sum test and adjusted p-values with Bonferroni correction. Based on annotation at cell type level, fibroblasts, keratinocytes, lymphoid cells and myeloid cells were selected for individual sub-clustering. For each cell type, multiple clustering resolution parameters were tested to identify cell subtypes, and the final parameter value was set based on the interpretability of the resulting clusters (fibroblasts: 0.2; lymphoid: 0.8; myeloid: 0.9; keratinocytes: 0.3), whereby gene markers for each subtype were identified following the same procedure explained above. Proportion tests for cell types and subtypes between Resting and Infectious bite were performed using the scProportionTest R package 12 by applying a permutation test that associates a p-value to the cell fraction fold-change.

Data_submission file includes: barcodes, features and matrix for rest, inf1 and inf2 samples as well as metadata file.