RNA Sequencing and Analysis

RNA from sorted BMDCs was extracted using the Macherey-Nagel™ NucleoSpin™ RNA Plus kit (Macherey-Nagel™, Cat# 740984.250) and quantified using the Qubit RNA HS Assay Kit (Thermo Fisher Scientific, Cat# {"type":"entrez-protein","attrs":{"text":"Q32852","term_id":"75319324","term_text":"Q32852"}}Q32852). Quality was measured by capillary electrophoresis using the Fragment Analyzer and the ‘Total RNA Standard Sensitivity Assay’ (Agilent Technologies, Inc. Santa Clara, USA, Cat# DNF-471-0500). All samples in this study showed high quality RNA Quality Numbers (RQN > 9.4). Library preparation was performed according to the manufacturer’s protocol using the ‘VAHTS™ Stranded mRNA-Seq Library Prep Kit’ for Illumina^®. Briefly, 500 ng total RNA were used for mRNA capturing, fragmentation, the synthesis of cDNA, adapter ligation and library amplification. Bead purified libraries were normalized and finally sequenced on the NextSeq550 system (Illumina Inc. San Diego, USA) with a read setup of 1x76 bp. The bcl2fastq2 tool was used to convert the bcl files to fastq files.

The reads of all probes were adapter trimmed (Illumina TruSeq) and the clean reads were analyzed using FastQC software to identify potential issues with data quality. The clean reads were then mapped to the mouse reference genome (Mus musculus, GRCm39/mm39) using STAR software. The percentage of uniquely mapped reads were greater than 80%. The uniquely mapped reads to each gene were counted using featureCounts. In order to assess the sample quality, we performed the principal component analysis (PCA) and hierarchical clustering for all samples. No batch effect was detected. The differently expressed genes (DEGs) (|log2FC| >= 1, FDR < 0.05) between non-stimulation and BEA, LPS or BEA with LPS stimulation following the previously described methods (23) were identified using DEseq2 package. DEGs expression was visualized as clustered heat maps using pheatmap package. The functional enrichment analysis (KEGG pathways and GO terms) of DEGs was carried out using enrichR package. Gene Set Enrichment Analysis (GSEA, Version 4.0.3) was used to identify enriched functional gene sets based upon the definitions of the Molecular Signatures Database (24, 25). The included gene set collections were “C2 curated gene sets”, “C5 ontology gene sets” and “C7 immunologic signature gene sets”. An enrichment map of significantly enriched gene sets was produced via Cytoscape (Version 3.8.0) (26) and the GSEA Enrichment Map plugin (27). Since Cytoscape defines a FDR of 0.25 as significant, this value was used as a cut off for inclusion into the network. In the enrichment networks, nodes represent gene sets, while edges represent mutual overlap between gene sets. Genes with overlapping genes and functional annotations were clustered manually to highlight the functional results. These clusters were encircled and labeled with an encompassing terminology. To achieve a simplified and more precise figure all clusters with less than three signatures were discarded from the network.