Metagenomic Analysis of the Gut Microbiome in Psoriasis Reveals Three Subgroups with Distinct Host-Microbe Interactions


 Background

Psoriasis is a chronic, inflammatory skin disease that impacts 2–3% of the global population. Besides skin manifestations, psoriasis patients have increased susceptibility to a number of comorbidities including psoriatic arthritis, cardiovascular disease, and inflammatory bowel disease. To understand the systemic component of psoriasis pathogenesis, we examined the fecal metagenome, host colonic transcriptome, and host peripheral blood immune profiles of psoriasis patients and healthy controls.
Results

Our analysis showed increased functional diversity in the gut microbiome of psoriasis patients. In addition, we identified microbial species that preferentially associate with psoriasis patients and which have been previously found to associate with other autoimmune diseases. Intriguingly, our data reveal new insights about psoriasis disease heterogeneity, as we identified three psoriasis subgroups that have distinct microbial and host features. Finally, integrating microbial and host features revealed host-microbe interactions that are specific to psoriasis or particular psoriasis subgroups.
Conclusion

To our knowledge, this is the first study that provides a comprehensive view of gut microbial function and corresponding host response in psoriasis patients. Our findings provide insight into factors that may affect the development of comorbidities in psoriasis patients and may hold diagnostic potential for early identification of psoriasis patients at risk for these comorbidities.


Introduction
Psoriasis is a prevalent immune-mediated disease characterized by in amed skin lesions and epidermal hyperproliferation. 2-3% of the global population is affected by psoriasis and the disease is known to be heterogeneous and multifactorial. It has been well established that psoriasis has different subtypes based on distinct disease characteristics. Plaque psoriasis (psoriasis vulgaris) is the most prevalent psoriasis subtype, with others being guttate, inverse, erythrodermic, and pustular. Psoriasis is driven by the interplay of genetics, the microbiome, diet, and other environmental factors. At the molecular level, aberrant activation of IL-17 signaling pathway is one of the major contributors to the disease. Therapeutically blocking components in IL-17 signaling pathway usually controls skin in ammation. However, the effectiveness of psoriasis treatments varies across patients, further highlighting the heterogeneous nature of psoriasis.
In addition to the skin manifestations, psoriasis patients are at higher risk of developing a number of other comorbidities. Psoriatic arthritis is one of the most prevalent comorbidities among psoriasis patients, as up to one-third of psoriasis patients transition into psoriatic arthritis [1]. Psoriasis patients also have 2.5 fold higher risk of developing Crohn's disease and 1.7 fold higher risk of developing ulcerative colitis linking psoriasis to in ammatory bowel disease (IBD) [2]. Other common psoriasis comorbidities include cardiovascular disease [3] and type 2 diabetes [4]. Although the association between psoriasis and its comorbidities is rmly established, the underlying cause and triggers are yet to be elucidated. Developing these comorbidities in psoriasis patients not only increases disease burden for the patients, it also complicates strategies in treatment and diagnosis. The ability to risk-stratify psoriasis patients for the development of certain comorbidities would greatly enhance strategies towards prevention, early detection, and treatment.
Harnessing the human microbiome for disease management and prevention has gained increasing interest since the initiation of human microbiome project.
The human gut microbiome consists of trillions of micro-organisms including bacteria, fungi and viruses in the gastrointestinal tract. These micro-organisms and their hosts have co-evolved into reciprocal relationships. The gut microbiome impacts many host physiological systems including immune system and neurological system, while host genetics and diet shape the composition and function of the gut microbiome. Dysbiosis in the gut microbiome has been observed in several autoimmune diseases ranging from IBD to rheumatoid arthritis to multiple sclerosis suggesting the systemic impact of these microorganisms. The causal role of gut microbiome in some of autoimmune diseases has been con rmed in animal models, highlighting the importance of the gut microbiome in disease pathogenesis [5].
Together, these ndings have implicated the gut microbiome as an important contributing factor for psoriasis pathogenesis and emergence of psoriasis comorbidities. To better understand the systemic component of pathogenesis associated with psoriasis, we performed a detailed multi-omics analysis with a cohort of psoriasis patients and healthy individuals. We utilized shotgun metagenomic sequencing to pro le both taxonomic composition and functional capacity of the gut microbiome. In addition, we carried out RNA-seq to pro le the host intestinal transcriptome. Host circulating blood was also collected to measure both systemic immune populations and their cytokine-producing capacity. Our study revealed several microbial features associated with psoriasis. Next, using clustering analysis, we identi ed three psoriasis subgroups, each with distinct microbial features that may predispose to certain psoriasis comorbidities. Lastly, we performed multi-omics analysis by integrating all our omics data and revealed disease and subgroup speci c host-microbe interactions. Our work highlights the heterogeneity of psoriasis and the potential role of microbial and host-microbe interactions in psoriasis pathogenesis and comorbidities.

Cohort Recruitment and sample collection
Adult psoriasis patients and healthy volunteers recruited from the San Francisco Bay area were enrolled in the study after providing informed consent. Individuals with abnormal coagulation parameters, positive HIV screening test, history of bleeding disorders, abdominal surgery, gastrointestinal cancer, in ammatory bowel disease, AIDS, immunode ciency or immunosuppressive medications, or concurrent in ammatory skin condition were excluded. All psoriasis patients had a diagnosis of psoriasis from a physician for at least 6 months prior to study enrollment, which was veri ed by study staff. To assess the psoriatic microbiome in an untreated state, subjects were excluded if they had received systemic biologic therapy in the last 6 months, non-biologic systemic medications (methotrexate, cyclosporine, corticosteroids, cyclophosphamide, retinoids, photochemotherapy) or antibiotics in the last month, or phototherapy or topical therapy in the last 2 weeks prior to the clinical visit. Healthy volunteers had no personal or family history of psoriasis. Gut biopsies, stool samples, and blood samples were collected on the same day. Gut biopsies were taken from the sigmoid colon of each subject. Two biopsies were preserved in RNAlater solution (Ambion) and kept at 4 °C overnight prior to storage at -80 °C. Stool samples were collected prior to gut biopsies and stored in -80 °C. Blood samples were taken on the same day as and preceded the gut biopsies. PBMCs were isolated from whole blood and stored in liquid nitrogen. We collected a total of 48 stool samples (15 healthy subjects and 33 psoriasis patients) for shotgun metagenomics and a total of 42 gut biopsies and blood samples (16 healthy subjects and 26 psoriasis patients) were collected for host transcriptomics and immune pro ling.

Diet habit analysis
We assessed dietary habits of healthy and psoriasis patients via in person distribution of the 2009-2010 National Health and Nutrition Examination Survey (NHANES) dietary screener questionnaire. The questionnaire consists of 25 questions about nutrient intake frequencies of different food groups. The differences in intake frequencies of each food group were compared between disease status and psoriasis subgroups with Fisher's Exact Test, and results with p < 0.05 were considered signi cant. To further assess quality of dietary habits, we summarized diet survey by calculating a Mediterranean diet score that would quantify the overall quality of diet. The traditional Mediterranean diet has been shown to reduce the risk for developing heart disease, metabolic syndrome, diabetes, and depression. The foundation for the Mediterranean diet consists of an abundance of fresh fruits, vegetables, whole grains, legumes, and omega-3 fatty acids and limited food like milk, cheese, red meats, and sweets with added sugars. We adopted Mediterranean score calculation from a previous study [53]. Brie y, we used 6 components of the Mediterranean diet (non-re ned cereals, fruit, vegetables, potatoes, whole grains, and legumes). For the consumption of food groups that aligned with the Mediterranean diet, we assigned scores 0, 1, 2, 3, 4, 5 when a participant reported no consumption, rare (1x/month), frequent (2-3x/month), very frequent (1x/week), weekly (2-6x/week), and daily (1-2x/day), respectively. For non-Mediterranean diet food components, we assigned scores on a reverse scale with scores 0, 1, 2, 3, 4, 5 when a participant reported daily (1-2x/day), weekly (2-6x/week), very frequent (1x/week), rare (1x/month), and no consumption, respectively. All scores were then combined to produce a total score that ranged from 0-50, where the higher the score, the better the diet. To test for statistical signi cance between both participant groups, we utilized the Wilcox test to compare the means and p-values of our data.
Shotgun metagenomics DNA was extracted from fecal samples with the modi ed cetyl trimethylammonium bromide (CTAB) method. The quality of the extracted DNA was assessed by gel electrophoresis and Quant-iT PicoGreen dsDNA Assay (Life Technologies) The extracted DNA was submitted to Vincent J. Coates Genomic Sequencing Laboratory at the California Institute for Quantitative Biosciences (www.qb3.berkeley.edu/gsl) for metagenomic library construction and sequencing. The metagenomic libraries were constructed using PrepX DNA library kit (Takara bio) and sequenced on HiSeq4000 (Illumina) for pair-end 150 bp sequencing.
Average of 22 million pair-end total reads were generated per sample. Reads were rst mapped to human genome (GRCh38, GENCODE release 25) by bowtie2 (v2.2.8) to identi ed human DNA. The non-human reads were extracted and sorted by using samtool (v0. 1.19). This data processing work ow results in average of 16 million non-human reads per sample. The resulting non-human reads were used for the subsequent taxonomic and functional pro ling.
The abundance table of bacterial species were extracted for the subsequent analysis. Functional pro les were generated by HUMAnN2 (v0.11.1), which is an assemble-free method to infer the functional capacities for each microbiome [54]. For each sample, HUMAnN2 maps reads to the pangenomes of species identi ed in the sample and performed additional translated search on unclassi ed reads. The resulting dataset is the abundance table of microbial gene families (UniRef90), which is summarized into the higher-level MetaCyc gene pathways. Our analysis identi ed total of 339 bacterial species, 409 MetaCyc pathways, and 650048 UniRef90 gene families in our dataset. Subsequent metagenomic analyses were performed using R package phyloseq [55]. The abundance of microbial species, UniRef90 gene families, and MetaCyc pathway were normalized to relative abundance to account for different library sizes and input as phyloseq objects and used for calculating microbial community diversities. Due to the high dimensionality and high drop-out rate of the microbial gene families table, we focus our analysis on microbial gene families with at least 3 counts for at least 80% of the total samples which results to 3356 UniRef90 gene families. Alpha diversity was calculated as four indices: Observed (number of observed microbial features), chao1, Shannon and Simpson index which summarize both richness and evenness of the community. The alpha diversities were compared by disease status and psoriasis subgroups with two-sided Wilcox test. Dissimilarities between each microbiome isolated are represented by Bray-Curtis dissimilarity matrix and visualized on a PCoA plot.
Un ltered counts per million (cpm) tables were used for differential abundance analysis. Differential abundant microbial features between psoriasis samples and healthy samples were identi ed using DESeq2 (v1.18.1) [56] with design ~ gender + age.bin + batch + Status. Microbial features with adjusted p-value < 0.05 and absolute log2 fold change > 0.6 are considered statistically signi cant. Differential abundance microbial features with none zero count for at least 10 samples were plotted on the heatmap to exclude features with high drop-out rates. Differential abundant microbial features associated with each subgroup identi ed in this cohort were identi ed by pair-wise comparison of each subgroup using DESeq2(v1.18.1) [56]. Microbial features with adjusted p-value < 0.1 and absolute log2 fold change > 0.6 are considered statistically signi cant.

Clustering analysis to identify subgroups in the cohort
The subgroups in this cohort were identi ed by complete linkage hierarchical clustering on the differentially abundant microbial UniRef90 gene families identi ed as descripted above. The clustering was done using R package pheatmap (v.1.0.12), which implement hclust for hierarchical clustering or the gapstat_ord function in phyloseq (v. 1.22.3). Gap statistics on differentially abundant UniRef90 genefamilies estimated three clusters present in our data set. Gap statistics was done by using clusGap function from cluster package (v.2.0.7-1) with bootstrapping 1000 times (B = 1000). The optimal number of clusters was determined as the smallest k that are within one standard error from the local maximal as suggested in Tibshirani et. al. 2001 [56]. The hierarchical dendrogram were cut by cutree function for 3 clusters (k = 3) in stats package (v.3.4.2) to assignment membership of clusters for each sample.
Host RNAseq RNA was extracted from sigmoidal colon biopsies by RNeasy Plus Universal Kits (Qiagen). Total RNA integrity was assessed by using Agilent 2100 Nano and Pico RNA kit, and all samples have RIN score greater than 7. cDNA libraries construction and sequencing were done by BGI (Beijing Genomics Institute). In brief, ribosomal RNA was removed from the total RNA by Ribo-Zero Gold Kit (Illumina) and library was constructed by TruSeq Stranded Total RNA Library Prep kit (Illumina). cDNA libraries were sequenced on HiSeq4000 (Illumina) for paired-end 100 bp sequencing. The average sequencing depth is 66 million per sample. Sequencing reads were mapped to human genome (GRCh38, GENCODE release 25) by using STAR (v. 2.4.2a). Number of counts mapped to each transcript are summarized by HTSeq-count (v. 0.6.0) and used as input for differential analysis in DESeq2 (v1.18.1) [56]. Raw read counts were normalized by median of ratios method using DESeq2. Differentially expressed genes between psoriasis samples and healthy samples were determined by DESeq2 with model design = ~ gender + age + batch + status. Genes with adjusted p-value less than 0.05 and absolute log2 fold change greater than 0.6 are considered signi cantly different between psoriasis samples and healthy samples. Similar analysis was performed to identi ed differentially expressed genes between different subgroups (genecluster.group) in this cohort with DESeq2 model design: ~ gender + age + batch + genecluster.group. Results of pair-wise comparisons were extracted, and signi cant differentially expressed genes were determined with the same signi cant criteria. The expression of the top 5000 most variable genes were transformed by variance stabilizing transformation for principle component analysis.

In silico Cytometry Methods
The colon immune cell composition was imputed from the sigmoid colon bulk RNAseq data by CIBERSORTx (https://cibersortx.stanford.edu/ ) [40] In brief, single cell single cell RNAseq dataset (scRNAseq) associated with sigmoid colon was extracted from colon immune atlas of Gut Cell Atlas project by Seurat [41]. The rst part of the CIBERSORTx work ow computed the signature matrix of immune cell populations in sigmoid colon using sigmoid colon scRNAseq dataset. The second part of the CIBERSORTx work ow imputed cell fractions in sigmoid colon by comparing bulk sigmoid colon RNA-seq to sigmoid colon immune cell signature matrix with the following parameters: relative run mode, 100 permutations, disabled quantile normalization, and S-mode batch correction which was implemented to account for cross-platform deconvolution. The cell compositions were compared by disease status and psoriasis subgroups with two-sided Wilcoxon test or Kruskal-Wallis test, and results with p < 0.05 were considered signi cant. Statistical analyses were performed with R.

Host ow cytometry
Flow cytometry for cell population Frozen PBMC were thawed, counted and 1 million cells were plated in a 96 well V-bottom plate. The cells were centrifuged 5 minutes at 400 x g and surfaced stained with a pre-determined antibody staining cocktail for 15 minutes at room temperature. The cells were washed with sort buffer, centrifuged and resuspended in BD Perm/Wash Buffer and centrifuged for 5 minutes at 700 x g. The supernatant was aspirated, and the cells were stained for FoxP3 by adding diluted uorophore-labelled anti-human FoxP3 antibody in BD Perm/Wash Buffer for 30 minutes at room temperature. Cells were centrifuged 5 minutes at 700 x g, the supernatant aspirated and the wash repeated with BD Perm/Wash buffer. The cells were resuspended a nal time in 50 uL of PBS and held at 4C until analyzed on a BD LSR-II instrument. Brilliant Violet Buffer was added into each staining step. Fluorescent minus one (FMO) controls were used for setting gates. A minimum of 300,000 events were analyzed for each sample.
Flow cytometry for cytokine production.
Frozen PBMC were thawed, counted and 1 million cells were plated in a 96 well U-bottom plate. PBMC were left unstimulated in complete culture media (RPMI + 10% heat-inactivated FBS supplemented with pen/strep/L-glutamine) or stimulated with 5 ng/mL of Phorbol myristate acetate (PMA) and 0.5 ug/mL of ionomycin for 4 hours in the presence of brefeldin A (10 ug/mL). After 4 hours of stimulation, the cells were centrifuged and stained with antibody staining cocktail listed for 15 minutes at room temperature. The cells were washed with sort buffer, centrifuged and resuspended in BD Perm/Wash Buffer and centrifuged for 5 minutes at 700 x g. The supernatant was aspirated, and the cells were stained for intracellular antigens by adding diluted uorophorelabelled antibodies in BD Perm/Wash Buffer for 30 minutes at room temperature. Cells were centrifuged 5 minutes at 700 x g, the supernatant aspirated and the wash repeated with BD Perm/Wash buffer. The cells were resuspended a nal time in 50 uL of PBS and held at 4C until analyzed on a BD LSR-II instrument. Brilliant Violet Buffer was added into each staining step. Fluorescent minus one (FMO) controls were used for setting gates. A minimum of 300,000 events were analyzed for each sample. Reported data represents the stimulated cytokine expression minus the unstimulated cytokine expression referred to as background corrected (BC) cytokine expression. fcs les were analyzed by FlowJo (v.10) and the immune pro les were exported as txt les, which will be import into R for statistical analysis. Wilcox test was used to assess the signi cance of the differences between disease status or disease subgroups.

Multi-omic analysis
We integrated 6 different datasets collected from three different measurement types: shotgun metagenomic sequencing, ow cytometry, and host RNAseq from 14 healthy subject and 26 psoriasis patients. We included the following datasets for data integration: 1.) Bacterial species from shotgun metagenomics with non-zero count in at least 10 samples. 2.) Microbial MetaCyc pathway from shotgun metagenomics with non-zero count in at least 10 samples. 3.) Microbial UniRef90 gene families from shotgun metagenomics with at least 3 counts for 80% of the total samples with metagenomic data. 4.) cell population from ow cytometry 5.) Cytokine production capacity from ow cytometry. 6.) top 500 most variable genes plus all differentially expressed genes identi ed by DEseq2 from host RNAsEq. More information about datasets included in multi-omic analysis are summarized in Table S9. Pair-wised datasets from different measurement type were integrated by Spearman's rank-order correlation. We constructed interaction networks with only strong and robust correlations that t the following criteria: correlations need to have FDR adjusted p-value less than 0.1 and absolute correlation coe cient greater than 0.6. We required microbial or host features with a non-zero count for more than 70% of samples within the group of interest. Multi-omic analyses were done within all subjects (both healthy and psoriasis subjects), healthy subjects alone, psoriasis subjects alone, and within each psoriasis subgroups. The multi-omic networks were visualized using R package, igraph (v.1.2.4.1) The community modules within each network were identi ed using fast greedy modularity optimization algorithm [57].

Study design and cohort summary
We recruited a cohort of 33 psoriasis subjects not on systemic therapy and 15 age-and gender-matched healthy individuals (Table 1a) to study the microbial features associated with psoriasis and their potential contribution to psoriasis pathogenesis. All psoriasis patients were clinically diagnosed with psoriasis at the UCSF Psoriasis and Skin Treatment Center and had a mean Psoriasis Area and Severity Index (PASI) of 14.2 representing moderate-to-severe disease. To characterize microbial composition, we collected stool samples and subjected each sample to shotgun metagenomic sequencing that provide both taxonomic composition and functional capacity. For the subsequent analyses, we focused on bacterial species, microbial UniRef90 gene families and microbial MetaCyc pathways. To link the changes of microbial features in psoriasis gut to the host response, we collected biopsies from sigmoid colon and subjected these samples to RNA-sEq. We also isolated PBMC from blood samples to measure immune cell pro les and cytokine production capacity. Together, our study design (Fig. 1) provided a comprehensive survey on both host biology and microbiome capacity in a cohort of psoriasis patients in comparison to healthy controls.
Microbial diversity and community structure between psoriasis and healthy gut microbiome Gut microbiome dysbiosis has been previously associated with decreased microbial diversity. Low microbial diversity has been observed in several human diseases including IBD, obesity, and autism [11,16,17]. It has been hypothesized that losing microbial diversity rise from missing group of "bene cial microbes" in the gut microbiome, which can lead to many detrimental consequences such as loss control in growth of opportunistic pathogens and lack of production of bene cial microbial derived compounds. To compare the microbial diversity in psoriasis patients and healthy subjects, we calculated different alpha diversity indices to estimate overall diversity (Shannon), evenness (Simpson), and richness (chao1) of each community. We observed higher Shannon and Simpson diversity metrics in the gut microbiome of psoriasis patients driven by evenness but not richness of the microbial functional pro le ( Fig. 2A). In contrast, taxonomic diversity metrics at the species level are comparable between the psoriasis and healthy gut microbiome ( Figure S1A). Despite increased functional diversity, the overall microbial community structures were similar between psoriasis samples and healthy samples as no distinct clusters was observed in PCoA plots for taxonomic and functional pro le ( Figure S1B and Figure S1C). All psoriasis patients in our cohort had a normal appearing lower gastrointestinal endoscopic exam, so major differences in diversity and community structure in the psoriasis microbiome as observed with other gastrointestinal diseases might not be expected. Transcriptomic pro les in colonic biopsies were comparable between psoriasis patients and healthy subjects. No gene features were signi cantly different between these groups after correction for multiple hypothesis testing ( Figure S1D).

Identi cation of microbial features associated with the psoriasis gut microbiome
We hypothesized that even though the gut microbiome from psoriasis patients have seemingly normal overall microbial community structure and diversity, the differences between psoriasis and healthy microbiome may be in speci c microbial features. To identify microbial features that are differentially abundant between gut microbiome from psoriasis patients and healthy individuals, we performed differentially abundance analysis using DEseq2 [18] which estimates differential abundant features using negative binomial model after controlling for known confounding factors for gut microbiome such as gender, age and experimental batch. Our analysis revealed bacterial species and microbial gene families that are differentially abundant between microbiome in psoriasis patients and healthy individuals ( Fig. 2C and 3A, Table S2 and S3). No microbial MetaCyc pathways were identi ed as differentially abundant. Among bacterial species with higher abundance in the psoriasis gut microbiome, Bacteroides vulgatus and Parasutterella excrementihominis (Fig. 2D) have been associated with other immune-mediated diseases. Increased intestinal colonization of Bacteroides vulgatus and elevated Bacteroides vulgatus reactive serum antibodies have been reported in ulcerative colitis patients [19]. Moreover, colonization of Bacteroides vulgatus is su cient to promote colitis in several animal models, which further supports the active role of Bacteroides vulgatus in development of colitis [20][21][22]. Parasutterella excrementihominis belongs to a newly studied genera Parasutterella with limited literature. An increase in Parasutterella excrementihominis has been found in ileal submucosal tissue of Crohn's disease patients [23] and in fecal samples of irritable bowel syndrome (IBS) patients [24]. Our analysis also revealed a decreased abundance of Phascolarctobacterium succinatutens in psoriasis compared to the controls (Fig. 2D). As suggested by its name, Phascolarctobacterium succinatutens utilizes succinate and converts it into acetate or propionate [25]. Accumulation of succinate has been observed in fecal samples of IBD patients [26,27] and DSS induced colitis model in mice [28]. Moreover, reduced succinate utilizing Phascolarctobacterium was observed in both patients with Crohn's disease and ulcerative colitis [14]. It is intriguing that the gut microbiome of psoriasis patients in our cohort shares some microbial features associated with IBD or IBS despite no intestinal symptoms reported in our patients. At the systemic level, we also detected a trend of elevated active CD4 + effector T cells in PBMC samples in psoriasis patients compared to healthy controls ( Figure S2A), but no difference in cytokine production capacities ( Figure S2B). Similar IBDassociated intestinal microbial signatures and elevated activated CD4 + effector T cells in the circulation of psoriasis patients suggest these patients may be more susceptible to systemic in ammation.
Microbial gene family analysis reveals three psoriasis subgroups with distinct microbial and host features Although most microbial gene families that were differentially abundant between psoriasis patients and healthy controls had no functional annotations, hierarchical clustering analysis on these microbial gene families identi ed three distinct groups in our cohort ( Fig. 3A and 3B). To con rm these subgroups were not observed by chance, we performed bootstrapped gap statistics on both the Euclidean distance and Bray-Curtis distance, which con rmed the presence of the three distinct groups in the cohort (Fig. 3C, Figure S3). Among the three groups identi ed by clustering, Group 1 consists of a mixture of healthy and psoriasis samples (14 healthy samples and 15 psoriasis samples), Group 2 consists of all psoriasis patients (9 psoriasis samples), and Group3 consists of almost all psoriasis patients except for 1 healthy control (1 healthy sample and 9 psoriasis samples) (Fig. 3D). The clustering was not confounded by BMI, age, and gender or diet ( Figure S4A, S4B, S4C, S4G, Table S1). The clustering analysis revealed three distinct psoriasis subgroups, which highlight the heterogeneity of the disease. For the subsequent investigations of these psoriasis subgroups, we term the psoriasis subjects from these subgroups PSO1, PSO2, and PSO3. Speci cally, psoriasis samples in group PSO1 clustered closely with most of healthy samples while group PSO2 and PSO3 represent more psoriasis dominant subgroups (Fig. 3A). All psoriasis subgroups had similar disease severity, disease duration, and age of disease onset ( Figure S4D, S4E, and S4F). In addition, the microbial diversity at the level of species, gene families, and pathways were seemingly similar among the three psoriasis subgroups ( Figure S4H, S4I, S4J). However, we found that each psoriasis subgroup had a number of distinct microbial and host features (Fig. 4A), which together suggest the biologic relevance of these psoriasis subgroups. Some of the interesting features associated with each psoriasis subgroup are highlighted below.
Interestingly, both PSO2 and PSO3 are dominated by psoriasis subjects, suggesting these are two distinct psoriasis speci c subgroups. We identi ed several common microbial features shared by these psoriasis subgroups, as well as some microbial features that are distinct to each subgroup. While psoriasis samples have lower abundance in Phascolarctobacterium succinatutens compared to healthy controls (Fig. 2D), the reduced abundance of Phascolarctobacterium succinatutens is speci c to PSO2 and PSO3, but not PSO1 (Fig. 4A). In addition, samples in PSO2 and PSO3 are less abundant with Turicibacter sanguinis and unclassi ed Turicibacter species (Fig. 4A). Less abundant Turicibacter species has been implicated in pediatric Crohn's disease in several cohorts [29,30]. Samples in PSO2 are more abundant with Bacteroides xylanisolvens and less abundant with Prevotella copri, Streptococcus thermophilus and Coprococcus sp ART55 1 (Fig. 4A). On the contrary, samples in PSO3 have lower abundance in Ruminococcaceae bacterium D16 and Lachnospiraceae bacterium 1 1 57FAA (Fig. 4A). Reduced abundance of Ruminococcaceae and Lachnospiraceae have been associated with the gut microbiome of psoriatic arthritis [6]. It is worth noting that all patients in PSO3 reported having joint pain or swelling, whereas only a fraction of patients in PSO1 or PSO2 did ( Figure S4H). It is intriguing that the gut microbiome from patients in PSO2 and PSO3 share some microbial features observed in patients with IBD and psoriatic arthritis suggesting patients in these psoriasis subgroups may have increased association with these psoriasis comorbidities compared to patients in PSO1.
In addition to distinct taxonomic features, PSO2 has a distinct pro le in microbial functions, especially in microbial pathways. Both abundant and depleted pathways were found in PSO2 relative to other psoriasis subgroups suggesting a shift of microbial functions in PSO2 (Fig. 4A). Microbial communities in PSO2 have lower abundance in the arginine and polyamine biosynthesis superfamily, suggesting lower capacity for polyamine production (Fig. 4A). Polyamines are small cationic molecules that modulate essential cellular process through binding of anionic DNA, RNA, and proteins. Polyamines are derived from arginine and intestinal polyamines can be generated by dietary supplementation, host metabolism, or microbial metabolism. Polyamines play a crucial role in intestinal mucosa maintenance and resident immune cell development [31]. Interestingly, it has been demonstrated that spermine, a class of polyamine, reduces IL-18 cytokine secretion by inhibiting NLRP6 in ammasome activation [32]. Samples in PSO2 also have increased capacity in pyridoxal 5-phosphate biosynthesis, which is a biologically active form of vitamin B6. Vitamin B6 can served as co-factor for various metabolic enzymes that are involved in wide range of cell metabolism pathways [33]. Intestinal commensal bacteria and dietary supplements are crucial source of vitamin B6 since mammalian cells lack the ability to de novo produce these micronutrients. Vitamin B6 de ciency has been linked to several immune-mediated diseases including rheumatoid arthritis and IBD [31,34]. Our data suggest that gut microbial communities in PSO2 have increased microbe-derived vitamin B6 and reduced polyamine production ( Fig. 4A). In contrast to PSO2, we did not identify microbial pathways that are uniquely associated with PSO1 and PSO3.
From the perspective of host response, PSO2 and PSO3 patients have distinct intestinal transcriptomic signatures (Fig. 4A). Gut biopsies from PSO2 patients have increased expression in ATP13A5 and PDE9A and reduced expression in sulfotransferases, SULT1C2 and SULT1C3 (Fig. 4A). In contrast, most of the transcriptomic signatures associated with PSO3 have increased expression compared to other psoriasis samples (Fig. 4A). It is interesting that some of the PSO3 transcriptomic signatures are similar to the ones reported in IBD patients or IBD animal models. FOLH1 encodes folate hydrolase, which is essential for absorption of dietary folate. Elevated FOLH1 expression has been reported in intestinal biopsies of patients with IBD and inhibiting FOLH1 activity ameliorates IBD associated abnormalities in mouse models [35,36]. HTR3A encodes one of the serotonin receptors and increased HTR3A colonic expression has been observed in patients with Crohn's disease [37]. UCP2 encodes for mitochondrial Uncoupling Protein 2 and has been implicated in several autoimmune diseases [38]. Expression of UCP2 is elevated in a DSS-induced mouse model of IBD and the severity of IBD can be ameliorated by knocking down UCP2 expression by siRNA [39]. Together, our data link PSO3 to IBD through intestinal transcriptomic pro les.
To gain more insight into the host intestinal immune response, we deconvoluted the immune cell composition in sigmoid colon bulk RNAseq using the digital cytometry framework CIBERSORTx [40]. As part of the CIBERSORTx framework, we rst de ned gene signatures of various immune cell populations in sigmoid colon using single cell RNAseq (scRNAseq) generated from healthy sigmoid colon [41]. We then used CIBERSORTx to apply the gene signature matrix to the bulk RNAseq to infer the composition of cell populations. Our analysis revealed differences in intestinal immune cell composition in different psoriasis subgroups which may represent distinct intestinal immune responses. Compared to other psoriasis subgroups, sigmoid colon from patients in PSO3 has more abundant CD8 + T cells and less abundant NK cells and activated CD4 + T cells (Fig. 4B). Although PSO3 patients possess some similarities in gene expression signature to IBD, PSO3 patients have less abundant colonic NK cells and activated CD4 + T cells, which are thought to be key contributors to IBD through cytokine production [42]. Interestingly, increased mucosal CD8 + T cell frequency has been reported in patients with Crohn's disease [43] and CD8 + T cells can trigger intestinal in ammation in mice colitis model [44]. Our data suggest PSO3 patients may be linked to IBD and particularly Crohn's disease through a CD8 + T cell mediated pathogenesis.
In addition to transcriptomic signatures, we observed distinct features in circulatory immune pro les and in the gut transcriptome between the two psoriasis dominant groups. Patients in PSO3 have higher activated memory CD4 + effector T cells compared to PSO2 while the frequency of memory CD4 + effector T cells and total CD4 + effector T cells are comparable between these two psoriasis subgroups (Fig. 4C). Similarly, memory CD8 + T cells have more activated population in PSO3 patients (Fig. 4C). Despite low T cell activation, patients in PSO2 have a higher capacity to produce IL-22 by circulating CD8 + T cells and TNF-alpha by circulating gamma-delta T cells (Fig. 4C) compared to healthy controls or PSO1. Together, our data reveal differential underlying circulatory immune responses associated with PSO2 and PSO3. PSO2 patients have higher capacity for pro-in ammatory cytokine production, while PSO3 patients have higher baseline T cell activation.

Distinct correlations between psoriasis severity to microbial features in each psoriasis subgroups
To further understand the relationship between the psoriasis subgroups and psoriasis heterogeneity, we tested if microbial features in these psoriasis subgroups are signi cantly correlated with psoriasis parameters such as severity, disease onset, and duration. We detected no signi cant correlations between the microbial features and any of the disease information when we combined all psoriasis patients in the cohort (data not shown). However, the microbial features-disease parameter correlations were observed in individual psoriasis subgroups. In PSO1 patients, disease severity is correlated with microbial pathways involved in purine degradation (Fig. 4D). Several microbial pathways are positively correlated with disease severity in PSO2 patients including pentose phosphate biosynthesis and L-arginine biosynthesis (Fig. 4D). The microbial L-rhamnose degradation pathway is positively correlated with disease duration in PSO2 (Fig. 4D). The microbial L-isoleucine biosynthesis pathway is correlated with disease duration positively in PSO1 patients and negatively in PSO2 patients (Fig. 4D). No signi cant correlations between microbial features and disease information were observed in PSO3. The different correlations in psoriasis subgroups might re ect the differential impact of gut microbiome in modulating psoriasis progression.

Multi-omic analysis in psoriasis patients and subgroups reveals distinct host-microbe interactions
In order to gain a comprehensive view of the crosstalk between the gut microbiome and host biology, we constructed multi-omic correlation networks by integrating microbial and host features across different measurements (Table S9). A subset of 40 subjects from our cohort with complete measurements from shotgun metagenomic sequencing, gut RNAseq, and circulating immune pro ling were included in the multi-omic analysis. This multi-omic cohort consists of 26 psoriasis subjects and 14 age-and gender-matched healthy subjects (Table 1b). Multi-omic networks were constructed within each disease status and within each psoriasis subgroup to reveal disease speci c and subgroup speci c host-microbe relationships. We applied Spearman's rank correlations to features from different measurement types and identi ed strong and robust correlations using the following criteria: (1) the correlations had adjusted p-value smaller than 0.1, (2) the correlations had absolute correlation coe cient greater than 0.6, (3) the correlations had features with at least 70% non-zero counts to avoid correlations that are driven by few extreme data points. The resulting multi-omic networks consist of nodes that represent host or microbial features and edges that represent signi cant correlations between nodes. Compared to the multi-omic network in healthy subjects, the network in psoriasis patients are denser and consists more edges (Fig. 5A). The nal multi-omic network consists of 73 signi cant correlations within healthy subjects and 97 signi cant correlations within psoriasis patients (Table 2). Intriguingly, the multi-omic network in each psoriasis subgroup displayed a distinct network structure. Multiomic networks in PSO1 and PSO2 were more interconnected compared to the network in PSO3 (Fig. 5A). It is interesting that PSO3 had fewer edges and nodes compared to other psoriasis subgroups and all psoriasis samples combined (Table 2). Overall, this analysis revealed distinct multi-omic network structures in psoriasis patients and healthy controls, as well as among different psoriasis subgroups.
The different multi-omic networks might re ect the different host-microbe interactions in psoriasis and healthy subjects. We hypothesize that the unique hostmicrobe interactions might contribute to pathogenesis of psoriasis and its comorbidities. Indeed, we identi ed in ammation associated modules in the multiomic network of psoriasis patients. One of the modules in the psoriasis associated multi-omic network positively links microbial pathways involved in microbial ppGpp biosynthesis and pyrimidine ribonucleosides salvage pathways with circulating IL-17 production in both CD4 + effector T cells and regulatory T cells (Fig. 5B). IL-17 cytokines are key disease drivers in psoriasis and blocking IL-17 can effectively ameliorate in ammation in psoriatic skin [45]. The associations between these microbial pathways and IL-17 production were only observed in psoriasis samples despite these microbial pathways are also being abundant in healthy samples ( Figure S5A). Our analysis also revealed psoriasis speci c associations between colonic expression of two proin ammatory chemokines, CXCL1 and CXLC3, and three unannotated microbial gene families (Fig. 5C, Figure S5B). Increased expression of CXCL1 and CXCL3 is associated with an in ammatory state of the intestine [46]. Our result suggests potential microbial controls of colonic expression of CXLC1 and CXCL3 in psoriasis patients prior to developing intestinal in ammation.
Besides psoriasis speci c host-microbe interactions, our analysis also identi ed host-microbe interactions that are speci c in psoriasis subgroups. A multiomic network identi ed in PSO1 consists of a module that associates microbial gene families with activated non-memory CD4 + T cells and IL-17 production capacity in CD4 + effector T cells (Fig. 5D), suggesting potential microbial controls in T cell activation and effector function in this psoriasis subgroup. Indeed, 95 out of 124 host-microbe associations in PSO1 link features of host circulating immunity to microbial features, which is higher compared to networks in PSO2 and PSO3 patients ( Table 2). In PSO2, we identi ed many associations linking microbial features to colonic transcriptome (Fig. 5A) including genes associated with IBD. CXCL8/IL8 is a chemokine that attracts neutrophils to the site of in ammation and elevated expression of CXCL8/IL8 has been found in in amed ulcerative colitis sigmoid colon [47]. Our data show that the expression of CXCL8 is positively correlated with microbial putrescine biosynthesis pathway in PSO2 patients (Fig. 5E). In addition, microbial sugar metabolism pathways are negatively correlated with colonic expression of FOS and PDE4C (Fig. 5F). FOS encodes c-fos which is a key subunit of AP-1, which is an important transcriptional factor regulates immune function. It has been shown that cfos plays a protective role in suppressing in ammation in both LPS and DSS induced colitis models [48,49]. PDE4C belongs to PDE4 family, which is a family of phosphodiesterase that degrades cyclic AMP (cAMP), which is an anti-in ammatory signaling molecule. Inhibiting PDE4 activity has been demonstrated as an effective strategy to treat in ammatory disease including psoriasis and IBD [50]. Together, our results suggest potential microbial controls of host colonic immune response through microbial metabolism in PSO2 patients. Compared to multi-omic networks in PSO1 and PSO2, PSO3 associated network has fewer associations (Fig. 5A, Table 2). Despite the sparse network, we identi ed a module in PSO3 network negatively correlates microbial tetrapyrrole biosynthesis pathway with several genes in B-cell biology including CD79B, PAX5, TCL1A, and MYBL1 (Fig. 5G). Tetrapyrroles are metal binding compounds consist of four pyrrole rings. With different modi cations and metal bound, tetrapyrroles serve as different cofactors, such as heme, cobalamin (Vitamin B12), and coenzyme F430 and have crucial roles in regulating diverse cell functions. The role of tetrapyrroles in B-cell functions has not been described. However, IBD patients are found to have depleted tetrapyrroles in their fecal samples [51] and increased colonic humoral immune response [52]. The negative association between microbial tetrapyrrole pathways and B cell genes found in PSO3 patients reminiscent of changes in tetrapyrroles level and humoral immune response observed in IBD patients may suggest an increased risk of developing IBD for PSO3 patients.

Discussion
The pathogenesis of psoriasis is highly heterogeneous, which poses challenges in diagnosis and disease control. In addition, psoriasis patients are more susceptible to developing systemic comorbidities which increases disease burden in these patients and complicates treatment plans. Therefore, it is important to understand the heterogeneous and systemic nature of psoriasis pathogenesis. Our study is one of the rst efforts to study the heterogeneity of psoriasis pathogenesis through the lens of multi-omic datasets using shotgun metagenomic sequencing, host colonic RNAseq, and host circulating immune pro ling. In this cohort, we identi ed microbial features that are differentially abundant between psoriasis and healthy samples (Fig. 2C, Fig. 3A). In addition, clustering analysis of the differentially abundant microbial gene families revealed three psoriasis subgroups, suggesting that the gut microbiome contributes to psoriasis heterogeneity. Each psoriasis subgroup has distinct microbial and host features (Table 3). Interestingly, both PSO2 and PSO3 share some features with IBD patients and PSO3 shared additional microbial features with psoriatic arthritis patients. On a self-reported questionnaire, all patients in PSO3 reported to have symptoms of joint pain or swelling ( Figure S4H) while only one PSO3 subject was formally diagnosed with psoriatic arthritis (data not shown). Both IBD and psoriatic arthritis are common psoriatic comorbidities and our study suggests an intriguing possibility that these psoriasis subgroups might represent psoriasis populations with differential risk for developing comorbidities such as IBD and psoriatic psoriasis.
The human gut microbiome has been implicated in many human diseases through microbial pro ling. These efforts have revealed microbial species or functions that are enriched in the disease of interest but have provided little mechanistic insight about how these microbial species and functions interact with host cells and how they might contribute to host biology. To overcome this challenge, we used a multi-omic approach to combine shotgun metagenomics with gut RNA-seq and circulating immune pro les to gain better insights on the potential roles of the microbiome in modulating psoriasis pathogenesis and psoriasis comorbidities. Our multi-omic analysis revealed interesting associations between microbial pathways with circulating IL-17A production and expression of proin ammatory chemokines, CXCL1 and CXCL3, in colon. These host-microbe interactions are only observed in psoriasis patients but not in healthy subjects, suggesting psoriasis speci c host-microbe interactions might be crucial drivers of the pro-in ammatory changes in psoriasis patients. In addition to psoriasis speci c host-microbe interactions, our analysis also revealed speci c host-microbe interactions in each psoriasis subgroup in our cohort.
The overall structures of host-microbe interaction network are different in each psoriasis subgroup (Fig. 5A). The network in PSO1 is highly interconnected and dominated by associations linking microbial features with circulating immune responses. The PSO2 network is also interconnected and dominated by associations linking microbial features to colonic pro-in ammatory gene expression. PSO3 has a relatively sparse network compared to PSO1 and PSO2. The fewer associations in PSO3 might due to small sample size of this subgroup. Alternatively, it may re ect less microbial impact in this subgroup. Together, our data suggest the host-microbe interaction can be context dependent, whereby the same microbial species or function may have differing impact based on the disease state or disease subgroup.
In this study, we provided a comprehensive survey of host and microbial changes in a cohort of psoriasis patients and healthy controls. Our analysis revealed distinct psoriasis subgroups within our cohort and identi ed unique features and associations in psoriasis patients and psoriasis subgroups. We hypothesize that psoriasis subgroups identi ed in our study may have diagnostic potential for early identi cation of psoriasis patients with higher risk of developing comorbidities. We are aware of several limitations of our study that could be improved by future work. Our analysis is limited by modest cohort size, so it would be important to validate our ndings in a larger independent psoriasis cohort. In addition to validation, a larger cohort will allow for the application of machine learning approaches to better identify psoriasis subgroups. Due to the scope of the study, we focused our efforts on cross-sectional observations of our cohort. Host-microbe interaction can be relatively dynamic over time, and some of the features or associations identi ed may be relatively transient.
Collecting longitudinal data may assist in the identi cation stable associations compared to transient associations. In addition, a longitudinal dataset would allow better tracking of host-microbe interactions associated with the development of comorbidities. Even with these limitations, the ndings of this study provide valuable insights about complexity of psoriasis biology that can help us better identify psoriasis patients with higher risk of comorbidities and unique biological pathways.

Conclusion
Using a metagenomics approach, we described characteristics of gut microbial diversity and composition in psoriasis patients and healthy controls. This study identi ed three psoriasis subgroups with distinct enrichment of microbial gene families. Some of these psoriasis subgroups had similar host and microbial characteristics as patients with psoriatic arthritis and in ammatory bowel disease, which are common psoriasis comorbidities. Intriguingly, we revealed distinct host-microbe interactions within psoriasis patients as well as within each psoriasis subgroup. Our results suggest that the pathogenesis of psoriasis and its comorbidities is highly individualized. The ndings of this study are relevant as they may have implications for early diagnosis or prevention of psoriasis comorbidities, with the possibility of personalized disease management strategies. Our analysis framework can be applied to study disease heterogeneity in other immune-mediated diseases. Dr. Liao has received research grant support from Abbvie, Amgen, Janssen, Novartis, P zer, Regeneron, and TRex Bio. All other authors declare no con ict of interest.   Figure 1 Multi-omic study design In this study we collected 6 different datasets for multi-omic analysis: Shotgun metagenomic sequencing from stool samples generated pro les of (1) microbial species, (2) microbial gene families, and (3) microbial gene pathways. RNAseq from sigmoid colon biopsies generated (4) host colonic transcriptome data. Flow cytometry analyses from PBMCs generated (5)     Red box depicts Group1, blue box depicts Group2, and green box depicts Group3. (C) Gap statistics to determine the optimal number of clusters. The y-axis denotes the K statistics which represents deviations of the cluster in (B) from a simulated cluster generated from a reference distribution. The optimal cluster number is determined by the smallest cluster number within 1 standard deviation of the maximal K. The gap statistics was done via bootstrapping 1000 times. (D) Stacked bar plot represents distribution of disease status of the three groups identi ed by cluster analysis. The height of each bar represents the size of each group and the color represents the disease status with red for healthy subjects and blue for psoriasis samples.