extfieldsofstudy: Medicine Biology

s2fieldsofstudy: Biology Medicine

text:
Axial spondyloarthritis patients have altered mucosal IgA response to oral and fecal microbiota

Axial spondyloarthritis (axSpA) is an inflammatory arthritis involving the spine and the sacroiliac joint with extra-articular manifestations in the eye, gut, and skin. The intestinal microbiota has been implicated as a central environmental component in the pathogenesis of various types of spondyloarthritis including axSpA. Additionally, alterations in the oral microbiota have been shown in various rheumatological conditions, such as rheumatoid arthritis (RA). Therefore, the aim of this study was to investigate whether axSpA patients have an altered immunoglobulin A (IgA) response in the gut and oral microbial communities. We performed 16S rRNA gene (16S) sequencing on IgA positive (IgA+) and IgA negative (IgA-) fractions (IgA-SEQ) from feces (n=17 axSpA; n=14 healthy) and saliva (n=14 axSpA; n=12 healthy), as well as on IgA-unsorted fecal and salivary samples. PICRUSt2 was used to predict microbial metabolic potential in axSpA patients and healthy controls (HCs). IgA-SEQ analyses revealed enrichment of several microbes in the fecal (Akkermansia, Ruminococcaceae, Lachnospira) and salivary (Prevotellaceae, Actinobacillus) microbiome in axSpA patients as compared with HCs. Fecal microbiome from axSpA patients showed a tendency towards increased alpha diversity in IgA+ fraction and decreased diversity in IgA- fraction in comparison with HCs, while the salivary microbiome exhibits a significant decrease in alpha diversity in both IgA+ and IgA- fractions. Increased IgA coating of Clostridiales Family XIII in feces correlated with disease severity. Inferred metagenomic analysis suggests perturbation of metabolites and metabolic pathways for inflammation (oxidative stress, amino acid degradation) and metabolism (propanoate and butanoate) in axSpA patients. Analyses of fecal and salivary microbes from axSpA patients reveal distinct populations of immunoreactive microbes compared to HCs using the IgA-SEQ approach. These bacteria were not identified by comparing their relative abundance alone. Predictive metagenomic analysis revealed perturbation of metabolites/metabolic pathways in axSpA patients. Future studies on these immunoreactive microbes may lead to better understanding of the functional role of IgA in maintaining microbial structure and human health.

Axial spondyloarthritis (axSpA) is an inflammatory arthritis involving the spine and the sacroiliac joint with extra-articular manifestations in the eye, gut, and skin. The intestinal microbiota has been implicated as a central environmental component in the pathogenesis of various types of spondyloarthritis including axSpA. Additionally, alterations in the oral microbiota have been shown in various rheumatological conditions, such as rheumatoid arthritis (RA). Therefore, the aim of this study was to investigate whether axSpA patients have an altered immunoglobulin A (IgA) response in the gut and oral microbial communities. We performed 16S rRNA gene (16S) sequencing on IgA positive (IgA + ) and IgA negative (IgA -) fractions (IgA-SEQ) from feces (n=17 axSpA; n=14 healthy) and saliva (n=14 axSpA; n=12 healthy), as well as on IgA-unsorted fecal and salivary samples. PICRUSt2 was used to predict microbial metabolic potential in axSpA patients and healthy controls (HCs). IgA-SEQ analyses revealed enrichment of several microbes in the fecal (Akkermansia, Ruminococcaceae, Lachnospira) and salivary (Prevotellaceae, Actinobacillus) microbiome in axSpA patients as compared with HCs. Fecal microbiome from axSpA patients showed a tendency towards increased alpha diversity in IgA + fraction and decreased diversity in IgAfraction in comparison with HCs, while the salivary microbiome exhibits a significant decrease in alpha diversity in both IgA + and IgAfractions. Increased IgA coating of Clostridiales Family XIII in feces correlated with disease severity. Inferred metagenomic analysis suggests

Introduction
Axial spondyloarthritis (axSpA) is an immune-mediated inflammatory arthritis, which affects the sacroiliac and spinal joints, and is associated with microscopic lesions or inflammation in the gut (1). It includes ankylosing spondylitis (AS), in which imaging evidence of inflammation in spine or sacroiliac joints is observed, and non-radiographic axSpA in which radiographic changes are not observed. While the genetic association of AS with the major histocompatibility molecule HLA-B27 has been known for almost five decades (2), the mechanistic link remains elusive. Like several other complex polygenic diseases such as inflammatory bowel disease (IBD), diabetes mellitus and multiple sclerosis, gene-environment interactions and particularly host-microbiota interactions have been centrally implicated in pathogenesis (reviewed in (3)). The latter is well supported by the observations that patients with AS have a dysbiotic gut microbiota and over half exhibit subclinical bowel inflammation (4,5). Furthermore, a recent study found that gut microbial dysbiosis is associated with worst disease activity in patients with axSpA (6).
Previously, we have shown that HLA-B27 expression perturbs the gut microbiota in an experimental model of spondyloarthritis (SpA) (7)(8)(9). In this model, transgenic rats expressing multiple copies of human HLA-B27 in conjunction with human b2-microglobulin (HLA-B27 TG) develop spontaneous bowel and joint disease (10). Enhanced mucosal immune responses to the gut microbiota strongly correlate with disease severity in these animals. For instance, HLA-B27 TG rats with arthritis elicit a stronger microbiota-specific IgA response (increased frequency of IgA-coated microbes) than control HLA-B7 TG animals without joint disease (8). Recently, we have shown that healthy individuals carrying the HLA-B27 allele also have altered gut microbial composition (11). While multiple studies employing 16S rRNA (16S) and metagenomic sequencing have revealed gut microbial dysbiosis associated with various spondyloarthropathies (4,5,(12)(13)(14)(15), identification of a causal microbe has been elusive. Compositional analysis of the microbial community may not be sufficient to determine the interactions between microbes and the host immune response. To identify such bacteria that preferentially affect disease susceptibility and severity, we employed IgA-sequencing (IgA-SEQ), which couples flow cytometry sorting of IgA-coated bacteria and 16S sequencing, to identify IgA coated microbes (16). This technique has been used to determine pathologically relevant IgA coated microbes in Crohn's disease (CD)-associated SpA (17). Colonization of germ-free mice with IgA coated microbes, particularly with adherent-invasive E. coli (AIEC) pathotype, induced Th17 immunity through microbial metabolic enzyme propanediol dehydratase (17), which induced IL-1b production by mononuclear phagocytes, thus perpetuating inflammation (18). These findings identified immune cells as metabolic sensors and highlight the importance of IgA coated microbes and their metabolites as potential therapeutic targets in CD treatment. AS and CD share several aspects of clinical overlap including bowel, eye, and joint disease (19). However, this type of investigation of the mucosal immune response and metabolic function of IgA coated microbes and has not yet been reported in axSpA patients.
The gut is the largest immune organ, facilitating interaction between the gut microbiota and local mucosal immune response. We hypothesize that other mucosal sites that are home to microbes contribute significantly to systemic rheumatic disease, particularly the oral cavity. It has been shown that patients with rheumatoid arthritis (RA) have a higher prevalence of severe periodontitis as compared to healthy controls (20). Disease activity in RA also tracks with periodontal disease activity with a corresponding partial resolution of microbial dysbiosis upon treatment (21,22). A recent study reported increased incidence of periodontitis in axSpA patients in comparison with the healthy controls; however, they did not observe differences between the bacterial richness and diversity or community structure (23). Furthermore, it has not been elucidated if perturbation of the oral cavity mucosal barrier is also altered in those with axSpA. Our hypothesis thus extends beyond the gut and investigates if immune modulation of the mucosal barrier occurs in the oral cavity of patients with axSpA.
In this study, we evaluated IgA coating (IgA-SEQ) and traditional 16S sequencing of microbes in the saliva and feces from axSpA patients and healthy controls (HCs). IgA-SEQ analyses identified potentially immune-targeted microbes that were not observed to be differentially abundant using traditional 16S sequencing alone. In addition, differential IgA coating on some fecal and salivary microbes was associated with disease activity. Using a computational approach to predict the microbial function, we found that the IgA + fraction in the fecal and salivary samples showed perturbation of many predictive metabolites and metabolic pathways in axSpA patients. These potentially immune reactive microbes may reveal novel host-microbial interactions in health and disease. Taken together, this study adds to our understanding of HLA-B27-associated host immune response and its interaction with the gut microbes in axSpA.

Study design and participants
We performed a prospective cohort study at the Oregon Health and Science University (OHSU). Subjects were excluded from all cohorts if they were younger than age 18 years, pregnant, had a history of prior intestinal surgery or colon cancer, or had antibiotic use 6 months prior to sample collection. This study was approved by the Institutional Review Board at OHSU and written informed consent was obtained from all participants.
AxSpA cohort: Subjects with a diagnosis of axSpA were recruited from either the OHSU rheumatology clinic or via the Spondylitis Association of America (SAA). The SAA is a patient advocacy group with physicians on its board, and members with AS are highly engaged and informed. Supplementary Table 1 lists the detailed location and recruitment status of axSpA patients and HCs (local vs SAA) for the stool and saliva cohort. Medical records were reviewed by a rheumatologist (JTR) to confirm a diagnosis of AS/axSpA based on the modified NY (New York) Criteria and the ASAS (Assessment of Spondyloarthritis International Society) classification respectively (24,25). A total of 26 subjects with AS/axSpA were enrolled (12 subjects provided fecal samples, 9 provided salivary samples and 5 subjects provided both fecal and salivary samples). Of these, 22 patients had a diagnosis of AS in their rheumatological medical records (11 of those were confirmed with imaging review), while 4 had limited external medical records but were sufficient to verify axSpA, thus our patient cohort is classified as axSpA in this study even though AS is highly represented. All of the patient medical information was reviewed again to confirm axSpA diagnosis during study design. While none of the axSpA subjects had psoriasis or peripheral arthritis, one axSpA patient in the salivary cohort was also diagnosed with IBD. The Bath ankylosing spondylitis disease activity index (BASDAI) was administered to subjects with axSpA. Their current medications, including biologics, nonsteroidal anti-inflammatory drugs (NSAIDs) as well as clinical HLA-B27 testing, if known, was also documented.
Samples: All subjects enrolled at OHSU provided a blood sample for genomic DNA extraction and subsequent HLA-B typing. Subjects recruited through the SAA shipped fecal and/or salivary samples to our lab using next-day shipping packets. For microbiome analysis, fecal and salivary samples were snap frozen upon collection and stored at -80°C until further analyses. All processing thereafter was performed blinded to patient phenotype. Supplementary Figure 1 illustrates the study experimental overview.
HLA-B typing: Genomic DNA was extracted from anticoagulated blood by the salting out method performed in a shared core resource at the Casey Eye Institute, OHSU. The LABType XR HLA-B SSO typing kit from One Lambda Thermo Fisher (RSSOX1B) was used according to manufacturer's instructions in the OHSU Laboratory of Immunogenetics and Transplantation.

Microbial community analysis
To determine the microbial community structure, genomic DNA from the salivary and fecal samples was isolated using the methods described previously (9). Briefly, genomic DNA was isolated using the Power Soil DNA Isolation Kit (MoBio Laboratories Inc., Carlsbad, CA) according to manufacturer's instructions. Amplification of the 16S ribosomal RNA genes was performed using the 515-806 primers specified by the Earth Microbiome Project (https://www.earthmicrobiome.org/ ). Sequencing of the V4 region was accomplished with Illumina MiSeq using barcoded primers (26). For IgA-SEQ, fecal samples were homogenized in phosphate buffered saline (PBS), supernatant was collected and washed in PBS with 1% fetal bovine serum (FBS) and incubated with blocking buffer (PBS with 20% rat serum) for 20 min. The IgA + fraction was stained with phycoerythrin (PE) labeled mouse anti-human IgA (Miltenyi 130-093-128), and enriched using anti-PE magnetic activated cell sorting (MACS) beads, and the negative fraction was collected for IgAmicrobes. IgA + fraction was further purified through a FACS Aria (BD Biosciences). Genomic DNA from IgA + and IgAmicrobial fractions was isolated and the V4 region of the16S rRNA gene was sequenced as described above.
Sequencing data were demultiplexed and processed using the DADA2 package in R (27). This included quality filtering, denoising using the DADA2 algorithm, chimera detection and removal, and tabulation into amplicon sequence variants (ASVs). Taxonomy assignment was performed using the RDP classifier (for Phylum to Genus annotations) and exact matching (for Species level annotations) with the Silva database (v. 132). Phylogenetic trees were generated using phangorn and decipher. Further analyses were performed using microbiome data analysis packages in R (phyloseq, vegan, MicrobiotaProcess) (28)(29)(30). IgA-SEQ data were subjected to decontam (31), a software tool to identify and remove contaminating sequence features in 16S sequencing data, and analyses were performed using in-house written scripts (32), available on Github (https://github.com/KarstensLab/igaseq). Samples that had 7,500 reads per sample in both the IgAand IgA + fractions were used. The IgA index was calculated as the log ratio of the difference between IgA coated and IgA uncoated bacteria over the sum of IgA coated and IgA uncoated bacteria (IgA index= -(log(IgA + taxon)-log(IgAtaxon))/(log(IgA + taxon) + log(IgAtaxon))) (33). Alpha diversity was evaluated by using multiple indices to determine richness and evenness in the community. While the number of observed genera/ASVs evaluates richness only using counts, the Shannon Index contains a logarithmic values including both richness and evenness, and Inverse Simpson Index focuses mostly on evenness using sum of squared proportions (34). Beta diversity was evaluated using the phylogenetic tree-based metric, UniFrac metric (35) and visualized with Principal Coordinate Analysis (PCoA).

Predictive metagenomic analysis
To infer the metagenome of oral and gut microbiota from 16S rRNA amplicon sequences obtained, we employed PICRUSt2 (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States 2) (36). PICRUSt2 estimates the metagenomic contribution of gene families in bacteria and archaea using 16S rRNA amplicon sequencing data, while taking into account the copy number of 16S rRNA gene. Functional potential was estimated by mapping 16S rRNA sequencing data to predicted functional annotation via alignments and hidden state prediction models (31), using default parameters. We used the quantitative abundance counts that mapped the Kyoto Encyclopedia of Genes and Genomes (KEGG) ortholog (KO) functional categories and focused our analysis on pathways using the multi-organism database MetaCyc (37) and individual metabolites (KO), which are both default functional annotations of PICRUSt2. Overlap between KO is determined between various groups, and area proportional Euler plots were created using Eulerr package in R (38).

Statistical analysis
All statistical analyses were performed in R. For 16S sequencing and IgA-SEQ analyses, non-parametric tests (Wilcoxon Rank Sum) were used to assess within and between group differences for the alpha diversity measures. Differences in individual microbial relative abundances at the genus and ASV level were assessed with non-parametric tests (Wilcoxon Rank Sum) and corrected for multiple tests using the Benjamini-Hochberg method (39). Differences in overall microbial community structure were evaluated using Permutational Multivariate Analysis of Variance (PERMANOVA). Correlation between relative abundance and disease activity (BASDAI) was performed using Spearman correlation with adjusted p-values after correcting for multiple comparisons using the Benjamini-Hochberg method (39). LEfSe was used to determine the presence of significant discriminant taxa by pairwise comparison (40). For PICRUSt2 analyses, LEfSe was used to identify candidate biomarker features among both individual metabolites (KO) and pathways (MetaCyc) (40). The LEfSe method performs a series of three tests: (1) at the class level (representing diagnosis axSpA or HC) an ANOVA (Kruskal-Wallis sum-rank test) with an alpha of <0.05 was applied (2) at the subclass level (representing patient sex) a Wilcoxon test with an alpha of <0.05 was applied and (3) a linear discriminant analysis (LDA) is finally applied to features to estimate their effect size and provides a ranking order of features most likely to explain group differences. Comparisons were made between IgA + and IgAfractions from salivary and fecal microbiome data from axSpA patients and HCs.

Clinical characteristics defining axSpA patients
To analyze microbial differences between axSpA subjects and HCs, we performed 16S sequencing on the unsorted microbial community and IgA-SEQ to identify IgA + and IgAsorted fractions from the fecal and salivary samples. The fecal cohort had 31 subjects (17 axSpA patients and 14 HCs), and the salivary cohort had 26 subjects (14 axSpA patients and 12 HCs). While there were no significant differences in age or gender (P> 0.2 and 0.6 respectively) between the axSpA patients and HCs for the fecal samples, the HCs were younger (P= 0.018), and with less female representation (P= 0.038) than the axSpA patients in the salivary cohort. However, both axSpA patients and HCs had similar BMI in the fecal (P = 0.5) and salivary (P = 0.4) cohort. AxSpA patients had a mean disease activity (BASDAI) of 3.5 in the fecal cohort, and 2.0 in the salivary cohort. The fecal cohort was comprised of 47% axSpA patients on biologics such as adalimumab, etanercept or secukinumab, while 36% of the salivary axSpA cohort were on biologics (Table 1). Furthermore, 47% of the fecal cohort and 14% of the salivary cohort of axSpA patients were on NSAIDs. In comparison, 21% of fecal cohort and 17% of salivary cohort healthy individuals reported NSAIDs use (Table 1).

AxSpA patients have altered microbial composition in saliva
To characterize microbial diversity in the fecal and salivary samples from patients with axSpA in comparison with the HCs, we performed 16S sequencing on both fecal (Figures 1A-C) and salivary ( Figures 1D-F) samples. We did not observe any differences between the microbial diversity (alpha diversity, Figure 1A) and community structure (beta diversity, Figure 1B) in the fecal samples from axSpA patients in comparison to HCs. In contrast, salivary samples from axSpA   patients showed a decrease in all (Shannon P = 0.016, P FDR = 0.048, observed P=0.047, P FDR = 0.094, and Inverse Simpson P= 0.048, P FDR = 0.094) measures of alpha diversity ( Figure 1D), however the microbial composition (beta diversity, Figure 1E) of axSpA patients was not significantly different than HCs at the genus level. We also compared individual microbial differences between axSpA patients and HCs at the genus level for fecal and salivary samples. We observed many microbes trended towards altered relative abundance in the fecal (e.g., Anaerofilum, Fecalitalea, Fournierella; P < 0.05, P FDR < 0.79; Supplementary Table 2) ( Figure 1C) as well as in the salivary samples (e.g., Megasphaera, Selenomonas, Oribacterium, Aggregatibacter, Campylobacter, Ruminococcaceae_UCG-014, Treponema, Atopobium, Corynebacterium and F0058) ( Figure 1F) of axSpA patients in comparison to HCs (P < 0.05, P FDR < 0.52; Supplementary Table 2). These differences no longer reached statistical significance after multiple testing correction. Similar results were seen in both fecal and salivary samples at the species level as measured using ASVs ( Supplementary Figures 2A-F).
Notably, there was a significant difference between the salivary alpha diversity in axSpA patients and HCs (Shannon's Index, P= 0.0006, P FDR = 0.002), however microbial composition as assessed by beta diversity was not significantly different between axSpA subjects and HCs (permanova, P = 0.06) (Supplementary Figures 2D, E). We did not observe these differences in the fecal microbial diversity or community structure in axSpA patients in comparison with HCs ( Supplementary Figures 2A, B). Despite this, relative abundance comparisons of individual microbes at the ASV level showed a trend towards increased abundance (P<0.05, P FDR = 0.74; Supplementary Table 2) Figure 2F).

Enrichment of IgA-coated microbes in fecal samples of axSpA patients
Although we did not observe differences in the overall microbial composition in fecal samples as measured by 16S sequencing, we suspected that there may be differences in how the host immune system interacted with these microbes. To determine whether microbiota from the fecal and saliva samples can elicit an immune response, we performed IgA-SEQ on IgA + and IgAmicrobes sorted from saliva and fecal microbial communities of the axSpA patients and HCs. This analysis revealed a tendency towards increased microbial diversity of IgA + fraction in the fecal microbiome of axSpA patients as compared to the HCs, whereas the IgAfraction showed a tendency towards loss of diverse in axSpA patients as compared with HCs ( Figure 2A). On the other hand, salivary microbiome of axSpA patients showed a significant decrease in microbial richness as measured by alpha diversity in both IgAas well as IgA + fractions as compared to the HCs ( Figure 2E). The beta diversity comparison (IgA + and IgAfractions) between axSpA and HCs did not yield any significant differences (data not shown).
Further analyses of IgA-enriched (Increased relative abundance in the IgA + fraction as compared with the IgAfraction) and depleted (decreased relative abundance of that microbe in IgA + fraction as compared with IgAfraction) bacteria from patients with axSpA identified a greater number of unique IgA-coated microbial genera in fecal ( Figure 2C) and salivary ( Figure 2G Figure 2D). In the salivary microbiome, IgA coating of Actinobacillus was increased while that of Butyricimonas was decreased when compared to HCs ( Figure 2H). These analyses were also performed at the ASV level and similar results were observed (Supplementary

IgA-SEQ reveals immunoreactive microbes correlating with disease activity
To determine whether IgA enriched microbes in feces from axSpA patients were associated with disease, we correlated genus and ASV level microbes with the disease activity score (measured by BASDAI). The analysis showed a significant negative correlation between Ruminococcaceae UCG-003 (ASV level, Figure 3B) in the fecal community, and a significant negative correlation between genus Actinobacillus and Alcaligenes in the salivary microbes ( Figure 3C). In addition, we found many microbes trending towards a correlation between IgA enrichment and disease activity. Of these, genus Clostridiales Family XIII correlated positively, while the genera Catenibacterium, Intestinibacter and members of family Lachnospiraceae (UCG010) and Ruminococcaceae (UCG-002) trended towards a negative correlation [ Figure 3A, (P < 0.05, P FDR < 0.20)]. At the ASV level ( Figure 3B), we observed a trend towards negative correlation with [Bacteroides spp.,

Perturbation in microbial metabolites in axSpA patients and healthy subjects
To infer the microbial metagenome and metabolic pathways, we employed PICRUSt2, which uses marker gene data (16S sequencing) to predict microbial metabolic potential. Previous inferred microbial metabolic profiling of SpA patients (41) and HLA-B27 TG rats (9) using PICRUST has revealed perturbation of inflammatory pathways associated with disease and/or HLA-B27 status. Therefore, in this manuscript, based on our findings from HLA-B27 TG rats, we further investigated the contribution of IgA + or IgAfractions in the perturbation of microbial metabolic function. We found multiple metabolic pathways (MetaCyc) significantly increased (alpha<0.05 for class and subclass analysis) in the fecal IgA + fraction of axSpA patients ( Figure 4A) including biosynthesis of unsaturated fatty acids, isopropanol biosynthesis, polyamine biosynthesis, and pathways for amino acids biosynthesis (phenylalanine, tyrosine), and colonic acid biosynthesis with a decrease in superpathway of UDP-N acetylglucosamine derived O antigen biosynthesis and pyruvate fermentation to acetone in the IgAfractions ( Figure 4B). On the contrary, salivary IgA + fractions ( Figures 4C, D) showed decreased pathways for L glutamate degradation VIII (to propanoate), acetyl CoA fermentation to butanoate II, and degradation of amino acid pathways (leucine, histidine) in axSpA patients as compared with HCs. We also observed alteration in nucleotide biosynthesis pathways (pyrimidine deoxyribonucleotides de novo biosynthesis I, 5 aminoimidazole ribonucleotide biosynthesis I, guanosine ribonucleotides de novo biosynthesis) in both fecal and salivary fractions (Supplementary Figures 4A, C). In the salivary IgAfraction of axSpA patients, we observed an increase in enterobacterial common antigen biosynthesis, with a decrease in pathways for L glutamate and L glutamine biosynthesis, L glutamate degradation VIII (to propanoate), acetyl CoA fermentation to butanoate II, purine nucleobases degradation I (anaerobic), guanosine nucleotides degradation ( Figures 4C, D). Additional exploratory analyses of metabolic pathways were also performed (alpha <0.1 for class and alpha 0.05 for subclass analysis), which revealed alteration in pathways for various amino acids biosynthesis, fatty acid biosynthesis, propanoate degradation etc. (Supplementary Figures 4A-D).

Discussion
In this study, we examined the IgA coating of microbes in the fecal and salivary microbiome from axSpA patients and HCs. To our knowledge, this is the first study to focus on IgA-enriched and depleted microbes and their predictive metabolic contributions in axSpA patients. We found an enrichment of IgA coated microbes in the fecal samples from axSpA patients. Also, the salivary microbiome from axSpA patients had a lower number of IgA coated microbes despite showing significant microbial dysbiosis. Some of these IgA coated microbes showed an inverse correlation with the BASDAI. In addition, predictive microbial function revealed various inflammatory and metabolic pathways perturbed in fecal and salivary samples from axSpA patients. Taken together, our data suggest that IgA selectively marks inflammatory or immune reactive microbes, resulting in an altered microbial community structure and function in axSpA patients.
The present study combining 16S sequencing of the oral and fecal microbiome and the IgA coated microbes proved useful in deciphering whether dysbiotic microbes are differentially coated with IgA. Despite minor changes in the overall microbial diversity and composition in axSpA patients and HCs, we observed a significant change in the IgA coating of oral and fecal microbes in these patients as compared to HCs. Notably, IgA-SEQ revealed many microbes that were not differentially abundant using the 16S sequencing alone, e.g., Akkermansia, and members of family Ruminococcaceae and Lachnospiraceae. These microbes have been shown to be associated with clinical studies in SpA and IBD (12,42) as well as in experimental models of SpA (7)(8)(9). Interestingly, Akkermansia has been implicated as a pathogenic microbe in SpA and IBD (12,42); and a protective microbe in metabolic diseases and obesity studies (43). We also observed other IgA coated microbes in the fecal (e.g., Lachnospiraceae, Ruminococcaceae) and salivary (e.g., Prevotellaceae, Actinobacillus) samples from axSpA patients, which have been shown to be associated with various spondyloarthropathies (12,14). Furthermore, in a recent study on IgA coated microbes in CD-SpA patients also reported increased IgA coating of Escherichia/Shigella, Lachnospiraceae (17), as seen in axSpA patients from our study. Almost 50-70% of axSpA patients have reported to have subclinical gut inflammation (44), which is further supported by increased calprotectin levels reported in SpA patients (45). This could explain the enrichment of Escherichia/Shigella similar to that seen in CD-SpA patients (17). Therefore, our study highlights Akkermansia, Escherichia/Shigella and members of Lachnospiraceae, Ruminococcaceae and Prevotellaceae for eliciting an IgA response which suggests immunologic significance. Since these organisms are able to activate the host immune response as indicated by an IgA response, the targeted micro-organisms should be studied further for functional analysis to determine host-microbe interactions in disease development.
Recent studies have focused on specific microbes as a disease biomarker, patients with AS and SpA have shown correlation of disease activity (BASDAI) with increased relative abundance of Ruminococcus gnavus and Dialister, respectively in two different studies (12,13). Other studies have reported gut dysbiosis correlating with the fecal calprotectin levels (6,46), suggesting that using multiple markers of disease activity (e.g., BASDAI, fecal calprotectin) in future studies may improve correlation between disease and dysbiotic microbiota. Our current study correlates various IgA coated microbes with the disease activity (BASDAI), and has revealed many microbes with negative correlation between IgA coating and BASDAI score in the fecal and saliva fractions. However, a member of class Clostridiales Family XIII showed positive correlation between IgA coating and increased BASDAI score, without a significant alteration in its relative abundance in axSpA patients (data not shown). While IgA coated microbes have shown to exacerbate gut inflammation (16,17), IgA coated microbes from healthy individuals have been shown to protect from gut inflammation (47) and regulates microbial composition (48). Additionally, a gut commensal microbe Bacteroides fragilis can alter its cell wall to enhance IgA coating, allowing it to adhere to the gut epithelial cells, and discouraging colonization by potential pathobionts (49). These opposing roles for IgA may be due to the T cell dependent IgA production affected by Th-17 cells or Regulatory T cells (17,47). Furthermore, our studies correlating microbial relative abundance with dysregulated inflammatory genes in HLA-B27 TG rats revealed that the relative abundance of only a few microbes correlated with disease severity, and most microbes may increase or decrease in response to gut inflammation, and therefore may not be the causal microbes (50). Together, these studies suggest a potential for using IgA coating as a microbial biomarker for disease activity, since many of the microbial perturbations in axSpA patients might be in response to altered gut microenvironment.
Our results demonstrate alterations in the salivary microbial community associated with axSpA patients. While the salivary microbes in axSpA patients showed significant alterations in the microbial community and structure, we observed a smaller number of IgA coated microbes in comparison with the fecal samples. However, there are significant changes in the predictive microbial function in both IgA + and IgAfractions from saliva samples of axSpA patients. A possible explanation is that perturbation of oral mucosa may be a result of disease instead of driving the disease in axSpA patients.
To evaluate the functional potential of IgA enriched microbes in axSpA patients, we employed PICRUSt2 to predict microbial metabolites altered in IgA + and IgAfractions in feces and saliva. Despite secretory IgA being the most abundant immunoglobulin produced in our body (51), only 5-10% of total microbes are IgA coated (16,52). Despite the difference in abundance, most of the microbial metabolic perturbations in the fecal samples were associated with the IgA + fractions, which further highlight their importance in axSpA pathogenesis. We found that most microbial metabolites perturbed in axSpA patients belong to inflammatory pathways [e.g., isopropanol biosynthesis, degradation of amino acid pathways (leucine, histidine, tyrosine), degradation of various sugars (mannan, glucose, xylose)], nucleotide biosynthesis pathways (50, 53) and pathways for short chain fatty acid (SCFA) metabolism (e.g., glutamate degradation VIII (to propanoate), acetyl CoA fermentation to butanoate II) (8,50,53). We also found alterations in the biosynthesis of tetrahydrofolate, a cofactor required for the synthesis of amino acids and nucleic acids, which is interesting since we also saw altered biosynthesis of both amino acids and nucleotides in axSpA patients as compared to the HCs. Additionally, alteration in biotin [vitamin B7] biosynthesis pathway is seen in IgA + fractions of fecal samples in axSpA patients. Biotin deficiency was recently reported in IBD patients (54) and mice deficient in biotin have exacerbation in colitis, which is alleviated by biotin supplementation (55). These results show a perturbation in the metabolic biosynthesis and degradation pathways, which could be involved in disease development and progression.
While the pathways for tryptophan biosynthesis were not significantly different, various genes/metabolites in that pathways were significantly altered in axSpA patients in comparison with HCs. Recent studies from Berlinberg and colleagues also show increased gene abundance for various metabolites in the tryptophan pathways in patients with axSpA and axSpA-CD (42). Alterations in tryptophan pathway are also reported by Stoll and colleagues to be associated with pediatric SpA (41). Metabolomic analysis of the plasma samples from AS patients also found a decrease in tryptophan metabolites (56), whereas another study found increase of tryptophan metabolites in AS patients (57). These studies indicate that decrease in tryptophan or increase in tryptophan metabolism might play an important role in AS pathogenesis. Also, multiple genes belonging to SCFA metabolism were perturbed in the saliva and fecal microbiome in axSpA patients. Previous studies by Shao and colleagues (58) have shown that AS and RA patients have decreased production of SCFAs -butanoate and propanoate. In addition, Asquith and colleagues performed fecal metabolomic analysis in HLA-B27 TG rats and found many inflammatory and SCFA metabolism pathways were altered in these animals (59). Gut inflammation in these HLA-B27 TG rats was ameliorated by adding propanoate in their drinking water (59), further emphasizing the protective role of SCFAs. These studies suggest the need for an in-depth analysis of fecal and blood metabolites and their association with AS. Furthermore, metagenomic analysis of the IgA coated microbes will allow studying host-microbe interactions at the species or strain level.
Our study has some limitations. First, we did not observe loss of fecal microbial diversity in axSpA patients in comparison with HCs, which could be due to our small sample size and increased variability due to our cohort being recruited from across the US. Additionally, this may be indicative of axSpA being a true dysbiotic model, instead of a being associated with loss in microbial diversity as reported in IBD (60). For example, a large cohort of patients with juvenile idiopathic arthritis (JIA) (75 JIA patients and 32 controls) did not find significant differences in the alpha or beta diversity of gut microbes when compared to controls (61), despite having significant changes at the level of individual microbes. Another study performed metagenomic sequencing on fecal samples from patients with CD, axSpA, axSpA-CD and healthy controls did not observe any differences in the microbial alpha diversity between these groups (42). Second, our observations were limited to fecal samples, and a more severe microbial dysbiosis may be present in the intestinal mucosal microbiome as shown in our previous study in HLA-B27 TG rats (9). Despite the lack of/smaller changes observed in the microbial composition, we were intrigued by the ability of IgA-SEQ to reflect the immunological impact of these microbes.
One more limitation is the inherent microbial diversity which complicates microbiome comparisons between patients and healthy individuals. We are still in the phase of exploring the microbiome's relationship to human health. Therefore, microbiome analyses test several hypotheses to identify disease specific alterations, which inflates the chance of false positive findings. To account for this, while also not excluding interesting observations for researchers in the similar field, we report both pvalues (P) and adjusted p-values (P FDR ). Similarly, for the inferred microbial metabolic pathway analysis, we performed a stringent and relaxed cutoff as done in our previous microbiome studies (9).
Another limitation is the use of biologics in almost 50% of the axSpA patients in our stool and saliva cohorts. Treatment with various biologics have shown to be associated with alteration in the microbial community. In SpA patients, treatment with IL-17 inhibitors correlated with features of sub-clinical gut inflammation, associated with changes in the relative abundance of certain microbial and fungal taxa (62). Conversely, TNF inhibitors have shown to shift the microbial community similar to the healthy controls (63), which may also account for the lack of alpha and beta diversity measures in axSpA patients by using 16S sequencing alone. Additionally, 47% of axSpA patients were on some form of NSAIDs as compared to 21% NSAID use in HCs. Studies on AS patients have shown association of NSAID use with alteration of the gut microbial community (64), increased levels of fecal calprotectin (46), as well as with increased number of bacterial species (65). However, many of these patients still had an active disease as measured by the BASDAI scores, and our studies revealed distinct IgA coated microbes associated with axSpA despite the limitation of biologics and NSAIDs used in our cohort. Further studies with increased sample size will allow us to dissect the effect of biologic and NSAIDs use on the microbial community structure in axSpA patients. However, our essential observation that IgA-coating of bacteria identifies a potentially pathogenic subset of the microbiome is supported by multiple observations reported herein as well as by prior publications (16,17).
The etiology of axSpA is complex involving an interplay of genetic, environmental and microbial factors in disease development and severity. Our results highlight the value of identifying the IgA coated microbes, which are immune-targeted microbes, based on 16S sequencing of the IgA + and IgAfractions. Further in-depth examination of IgA coated microbes using shotgun metagenomic sequencing of IgA + and IgAfractions to determine potentially inflammatory microbes at a strain specific level is warranted. Our results suggest that IgA coated microbes in axSpA patients indicate a stimulated immune response to gut commensal microbes, which might contribute to the development of axSpA. Future studies to determine the mechanistic links between IgA coated microbiota and disease will further highlight the functional implications of IgA coated bacteria in axSpA.

Data availability statement
The data presented in the study are deposited in Sequence Read Archive (SRA), accession number under bioproject PRJNA873789.

Ethics statement
This study was reviewed and approved by OHSU Institutional Review Board. The patients/participants provided their written informed consent to participate in this study.
Author contributions MA ( †) was involved in the study design and data generation for the current manuscript. TG, TM, JR and LK had access to all the data and patient records and share the responsibly for data integrity and its accurate analysis. All authors contributed to the article and approved the submitted version.


Hemogenic and aortic endothelium arise from a common hemogenic angioblast precursor and are specified by the Etv2 dosage

Significance Hematopoietic stem cells (HSCs) are generated from specialized endothelial cells, called hemogenic endothelial cells (HECs). It has been debated whether HECs and non–HSC-forming conventional endothelial cells (cECs) arise from a common precursor or represent distinct lineages. Moreover, the molecular basis underlying their distinct fate determination is poorly understood. We use photoconvertible labeling, time-lapse imaging, and single-cell RNA-sequencing analysis to trace the lineage of HECs. We discovered that HECs and cECs arise from a common hemogenic angioblast precursor, and their distinct fate is determined by high or low dosage of Etv2, respectively. Our results illuminate the lineage origin and a mechanism on the fate determination of HECs, which may enhance the understanding on the ontogeny of HECs in vertebrates.

Hematopoietic stem cells (HSCs) are generated from hemogenic endothelial cells (HECs) in the floor of the dorsal aorta (DA) via endothelial-to-hematopoietic transition (EHT). Yet whether HECs and conventional endothelial cells (cECs) in the DA share a common precursor is controversial, and the molecular mechanisms governing their fate specification remain incompletely defined. Using a combination of fate mapping, timelapse imaging, genetic manipulation, and single-cell RNA sequencing, here we show that HECs and cECs display strictly spatial separation in the DA where nearly all the endothelial cells in the floor and roof are HECs and cECs, respectively. We further show that HECs and cECs in the DA arise from a common hemogenic angioblast precursor, which differentiates into HECs and cECs during axial migration prior to the DA formation. The specification of HECs and cECs from hemogenic angioblasts is governed by the Etv2 dosage by differentially regulating Fli1a, Notch, and Sclβ. Finally, we document that pan-endothelial overexpression of transcriptional factor runx1 is sufficient to promote HEC fate in the DA roof. Our study reveals the lineage origin of HECs and cECs in the DA and uncovers the molecular network controlling their fate specification.
endothelium | hemogenic endothelium | lineage origin | fate determination During vertebrate embryogenesis, hematopoietic stem cells (HSCs) are generated from a subset of specialized endothelial cells (ECs), named hemogenic endothelial cells (HECs), via endothelial-to-hematopoietic transition (EHT) in the dorsal aorta (DA) (1)(2)(3)(4). Timelapse imaging and in situ analysis of hematopoietic markers indicate that HECs are located exclusively in the ventral domain/floor of the DA in the aorta-gonad-mesonephros (AGM), budding toward intraaortic lumen in avians (5,6) and mice (2,7,8) or budding toward subaortic space between the DA and posterior cardinal vein in zebrafish (3,4). Although HSC-forming HECs and non-blood-forming conventional endothelial cells (cECs) are both located in the DA, whether HECs and cECs arise from a common bipotent precursor remains debatable. Several studies using human pluripotent stem cells or murine model have suggested that HECs and cECs represent distinct lineages indicated by the expression of distinct cell surface markers, such as CD73 and CD184 (9), or by their different requirements of Notch signaling strengths (10). Studies in the avian model especially suggested that different mesoderm origins, somatopleural and splanchnic mesoderm, give rise to cECs and HECs, respectively (11,12). In contrast to the distinct lineages model, other studies suggested a common precursor model that endothelial precursors appear to adopt arterial endothelial fate before committing into HECs (13,14). Consistent with this model, recent single-cell RNA-sequencing (scRNA-seq) studies (15)(16)(17) constructed the continuous developmental trajectory from arterial cECs to HECs, showing that this transition is accompanied with gradual loss of arterial signatures and acquirement of hematopoietic features. To clarify this controversial issue, delicate in vivo fate mapping studies are required. In parallel, although multiple signaling pathways including Hedgehog, Wnt, Vegf, Notch, and BMP are reported critical for HEC and cEC fate specification (18)(19)(20)(21), it seems that these signaling pathways eventually converge on the ECs-autonomous requirement of Notch signaling to activate key downstream transcriptional factors, including Gata2, Runx1, and Gfi1 (22)(23)(24). Besides, these key factors are also regulated by the upstream factors, including Etv2 and Scl (25)(26)(27). Precisely how these upstream factors interplay with Notch signaling and form a genetic network to control cEC and HEC fate specification remains unexplored.
Zebrafish is a prominent model for studying many aspects of developmental biology including angiogenesis and hematopoiesis (28,29). In zebrafish, the first wave of endothelial progenitors, named angioblasts, in the trunk region arises from the posterior lateral

Significance
Hematopoietic stem cells (HSCs) are generated from specialized endothelial cells, called hemogenic endothelial cells (HECs). It has been debated whether HECs and non-HSC-forming conventional endothelial cells (cECs) arise from a common precursor or represent distinct lineages. Moreover, the molecular basis underlying their distinct fate determination is poorly understood. We use photoconvertible labeling, time-lapse imaging, and single-cell RNA-sequencing analysis to trace the lineage of HECs. We discovered that HECs and cECs arise from a common hemogenic angioblast precursor, and their distinct fate is determined by high or low dosage of Etv2, respectively. Our results illuminate the lineage origin and a mechanism on the fate determination of HECs, which may enhance the understanding on the ontogeny of HECs in vertebrates.
plate mesoderm at ∼11 h postfertilization (hpf ), migrates toward the midline during 14 to 18 hpf, and finally coalesces to form the DA at 21 hpf (30,31). Similar to mammals, HECs in zebrafish are located in the floor of the DA and express several early hematopoiesis-related genes, including gata2b, sclβ, and runx1 (13,23,26). During 30 to 60 hpf, HECs undergo EHT and directly bud from the floor of the DA toward the subaortic space (3,4). Because of the unique advantages of in vivo imaging and genetic amenability, the zebrafish model system provides a powerful tool to dissect the lineage origin of HECs and cECs and the mechanism controlling their fate specification.
In this report, we combined in vivo photoconversion lineage tracing, time-lapse imaging, and scRNA-seq analysis to study the ontogeny of HECs and cECs. We showed that the ventral and dorsal endothelia of the zebrafish DA are largely HECs and cECs, respectively. We further documented that HECs and cECs arise from a common precursor, termed hemogenic angioblast, which undergoes HEC and cEC fate specification prior to the DA formation. Genetic studies further unraveled that at early stage, high and low dosages of Etv2 determine HEC and cEC fates through differentially regulating the expression and activity of fli1a, sclβ, and Notch. Finally, we showed that ectopic expression of transcription factor Runx1 in the roof of the DA can convert cECs into HECs.

Results
The Endothelia in the Floor of the DA Are HECs. To have a better understanding of the ontogeny of HECs and cECs, we generated a double transgenic line Tg(flk1:eGFP;flk1:NLS-Eos), in which all ECs in the DA express enhanced green fluorescent protein (eGFP) in the cytoplasm and photoconvertible Eos protein in the nucleus (32). This transgenic line allows us to label ECs in the roof or the floor of the DA and follow their behavior by timelapse imaging (Fig. 1A) during EHT from 28 to 72 hpf (3,4). We photoconverted 5 to 10 ECs per embryo in the DA floor at 26 hpf and immediately tracked their behavior by time-lapse imaging until 72 hpf (Fig. 1B). Unexpectedly, we found that ∼96% (211 out of 220) of the photoconverted floor ECs underwent EHT (Fig. 1D), of which ∼49% (104 out of 211) of the cells underwent EHT directly, ∼45% (96 out of 211) of the cells divided at least once before initiating EHT, and ∼5% (11 out of 211) of the cells burst into pieces, a phenotype similar to abortive EHT observed in runx1 knockdown embryos (3,26) (Fig. 1C). Whether the abortive EHT is a natural phenomenon or a consequence of experimental artifacts is unclear. Nonetheless, these data indicate that ∼96% of the ECs in DA floor undergo EHT, and most of them eventually become hematopoietic stem/progenitor cells (HSPCs). Since there are only ∼34 cells in the DA floor at  27 hpf per embryo, this result predicts that 33 of them are HECs, and only 1 cell is cEC (Fig. 1E). Considering that ECs are patterned around the DA circumference unevenly (31), we believe that the rare cECs which we observed in the DA floor in our photoconversion assay might result from the mislabeling of the lateral ECs of the DA. Based on these findings, we conclude that nearly all the ECs in the DA floor at around 26 to 28 hpf are HECs.
The Endothelium in the DA Roof Are cECs Lacking Hemogenic Potential. In avian and mammalian embryos, HECs are reported to undergo EHT toward the DA lumen to form intraaortic hematopoietic clusters (IAHCs) (2,5,33,34). These IAHCs are not only restricted to the floor but also observed in the DA roof of murine embryos (7,33,35). However, whether these IAHCs in the DA roof are generated in situ or arise from the floor remains elusive. Owing to the advantages of time-lapse imaging and high resolution of ECs labeling, the Tg(flk1:eGFP;flk1:NLS-Eos) fish provides a useful tool to address this issue. Therefore, we photoconverted the ECs in the DA roof at 26 hpf and monitored their behaviors ( Fig. 2A). Time-lapse imaging showed that, while the majority (∼67%; n = 103/154; fish = 21) of the ECs remained in the DA roof/lateral during 28 to 72 hpf, a substantial portion (∼33%; n = 51/154; fish = 21) of the ECs moved ventrally and repositioned to the floor of the DA, some of which underwent cell division before and during the ventral movement ( Fig. 2 B and C). Remarkably, the ventral movement of the roof ECs appeared to be temporally correlated with the EHT in the floor (Fig. 2D), indicating that the loss of HECs due to EHT in the DA floor is largely replenished by the ventral movement of roof ECs, which is consistent with the previous observation (3,36,37). However, direct budding of ECs from the DA roof was never observed, although two EHT events (∼1%; n = 2/154; fish = 21) were captured after the roof ECs moved ventrally and repositioned to the floor (Fig. 2C). Given that there are ∼35 ECs in the DA roof at 27 hpf per embryo, this result means that less than 1 cell from the DA roof undergoes EHT per embryo (Fig. 2E). Taken together, these data demonstrate that the ECs in the DA roof are all cECs lacking hemogenic potential.

HECs and cECs Arise from a Common Precursor and Separate
prior to the DA Lumen Formation. To delineate the mechanism governing the fate determination of HECs and cECs, we isolated single ECs from the DA roof and the floor respectively using the Tg(flk1:eGFP;flk1:NLS-Eos) zebrafish after photoconversion and then performed scRNA sequencing (Fig. 3A). These single ECs were isolated at two different time points, 21 and 28 hpf, when the lumen of the DA is established (31) and the EHT begins (3,4), respectively. Uniform manifold approximation and projection (UMAP) clustering analysis identified two clusters, Cluster roof and Cluster floor, enriched in the ECs from the roof and floor, respectively ( Fig. 3A and SI Appendix, Fig. S1A). Differentially expressed gene and gene ontology (GO) analyses showed that the feature genes, including flk1, cdh5, tbx20, sox18, efnb2a, hey2, and dll4, enriched in roof ECs (Fig. 3B) were associated with vasculogenesis, angiogenesis, and vessel development (SI Appendix, Fig. S1B), showing the vascular identity of these cells. In contrast, the enriched GO terms of floor ECs were associated with catabolic and metabolic processes (SI Appendix, Fig. S1B) and were highly enriched in HEC markers and definitive hematopoiesis-related genes, such as gfi1aa, scl, myb, dnmt3ba, ikzf1, gata2b, and runx1 ( Fig. 3B), showing that they are in an active metabolic state with hematopoietic potential. These data indicate that Cluster roof and Cluster floor represent cECs and HECs, respectively, and that their fates have already been separated by 21 hpf. To   further define the time window during which HECs and cECs are specified, we reanalyzed the Wagner's scRNA-seq database from early developmental stages of zebrafish embryos (38) at 10-, 14-, and 18-hpf stages. By integration analysis with our single-cell RNA-seq dataset (21-hpf stage), we found that the 21hpf cells from our dataset were predominantly (33/36; ∼92%) clustered with mesoderm cells defined in Wagner's study (38) (SI Appendix, Fig. S1C). To gain insights into these clusters, we further extracted the information of etv2 + or flk1 + cells, which presumably include HECs, cECs, and their precursors, angioblasts (30,39), from 10-, 14-, and 18-hpf stages. We then integrated the extracted etv2 + and flk1 + scRNA-seq dataset (referred to as Wagner's dataset in Fig. 3C) with our scRNAseq dataset (21-hpf stage) and reperformed clustering analysis. UMAP results identified three clusters, Cluster 1, Cluster 2, and Cluster 3 ( Fig. 3C and SI Appendix, Fig. S2A). We found that the roof ECs and floor ECs from the 21-hpf embryos were predominantly clustered in Cluster 2 and Cluster 3, respectively (Fig. 3C), suggesting that Cluster 2 and Cluster 3 represent ECs which have already committed to cEC and HEC lineage, while Cluster 1 likely represents noncommitted precursors. Indeed, this notion was further supported by the findings that the Cluster 1 contained a high percentage of etv2 single positive cells with feature genes involved in early developmental process, e.g., embryonic organ morphogenesis and anterior/posterior pattern specification (SI Appendix, Fig. S2 B-D). Moreover, the cEC marker dll4 and HEC marker gata2b were highly expressed in Cluster 2 and Cluster 3, respectively, but absent in Cluster 1 (SI Appendix, Fig. S2E). Next, we examined the distribution percentage of the three clusters at different time points of zebrafish development. The results showed that Cluster 1 consists of cells from 10-, 14-, and 18-hpf stages exclusively, whereas Cluster 2 and Cluster 3 consist of cells predominantly from 18-and 21hpf stages, with a small fraction of cells from 14 hpf (∼0.8%) in Cluster 2 ( Fig. 3D), suggesting that cECs and HECs in the DA are likely to arise from a common angioblast precursor, and their fates are largely determined during the period of 14 to 18 hpf, which correlates with the time window when somite-derived Notch signals were shown to specify HEC formation (40,41).
To further support our hypothesis, we performed single-cell fate mapping analysis using the Tg(flk1:eGFP;flk1:NLS-Eos) zebrafish. In this assay, we photoconverted a single flk1-Eos + cell in each embryo at 15 or 18 hpf and defined the identity of the photoconverted cell by examining the distribution pattern of its progenies (if the converted single cell produced two or more progenies) in the DA roof and floor at 26 hpf (Fig. 3E). We reasoned that if cECs and HECs indeed originated from a common precursor and their fates were specified between 14 and 18 hpf, we would expect to observe two different distribution patterns, bilocalization and unilocalization. The bilocalization indicates successful labeling of a common precursor capable of giving rise to both cECs and HECs, whereas the unilocalization indicates successful labeling of a lineage committed cell capable of giving rise to either cECs or HECs (Fig. 3E). In agreement with our hypothesis, single-cell fate mapping revealed that ∼55.5 and ∼44.5% of the cells displayed bilocalization and unilocalization potential at the 15-hpf stage, respectively (Fig. 3F ). As fish developed to the 18-hpf stage, the percentage of the cells with bilocalization potential dropped to ∼22.2%, accompanied by an increase of the percentage (∼77.8%) of the cells with unilocalization potential (Fig. 3F ), indicating a gradual fate commitment process from bipotent progenitors to unipotent lineage-committed cells during this developmental time window. These results, together with the scRNA-seq data analysis, demonstrate that cECs and HECs in the DA indeed originate from a common flk1 + precursor, which we refer to as hemogenic angioblast, and their fates are largely determined during the period of 14 to 18 hpf, well before the lumen formation of the DA.

Etv2 Dosage Determines HEC and cEC Fate through Differen-
tially Regulating the Activities of Fli1a, Notch, and Sclβ. Having defined the lineage origin of HECs and cECs, we next investigated the mechanism controlling their fate specification. The Ets variant 2 (etv2) has been identified as a master regulator for HEC and cEC development (25,42), yet the underlying mechanism remains incompletely clear. To dissect the role of Etv2 in HEC and cEC development, we generated a transgenic fish Tg(etv2:NLS-d2eGFP), in which a fast-degraded d2eGFP protein is under control of the etv2 promoter (43) (SI Appendix, Fig. S3A), and performed fluorescence-activated cell sorting (FACS) analysis at 16 hpf, when HECs and cECs are known to be specified (Fig. 3 C, D, and F ). Intriguingly, two etv2 + cell fractions, etv2-d2eGFP high and etv2-d2eGFP low , were detected (Fig. 4A). qPCR analysis showed that the expression of the HEC marker, gata2b (23), was significantly enriched in etv2-d2eGFP high fraction (Fig. 4A). This observation prompted us to hypothesize that Etv2 dosage might govern the fate choice between HECs and cECs. To test this hypothesis, we generated a double transgenic line Tg(etv2:mCherry-T2a-CreER T 2 ;flk1:loxP-DsRedx-loxP-eGFP) (SI Appendix, Fig. S3A) and performed lineage tracing analysis by pulsing embryos with a range of 4-Hydroxytamoxifen (4-OHT) doses between 15 and 17 hpf when HECs and cECs are undergoing fate specification (Fig. 4B). The resulting flk1-GFP + cells were examined at 32 hpf, which is the latest timing point used in the study to address the presence and localization of HECs in the aortic floor. Results showed that low (2.5 μM) and median (5 μM) dose treatment produced limited numbers of GFP + cells exhibiting bias distribution pattern toward the floor of the DA (Fig. 4B). On the other hand, high-dose (15 μM) treatment led to a significant increase of GFP + cells displaying comparable number in the floor and roof (Fig. 4B). In parallel, we also tested different lengths of time (1, 1.5, and 2 h) of 4-OHT treatment with median dose (5 μM). Results showed that three groups exhibited labeling bias toward floor of the DA, but the bias toward floor was decreased in the group of 2-h treatment (SI Appendix, Fig. S3B). These data suggest that the progenitors with high and low Etv2 dosages preferentially adopt the HEC and cEC fates, respectively. To further support this notion, we generated another transgenic line Tg(hsp70l:etv2-P2a-mCherry), in which Etv2-P2a-mCherry protein is under control of the hsp70l promoter (SI Appendix, Fig. S3A), and asked whether overexpression of Etv2 during 14 to 18 hpf would enhance HEC formation. Indeed, whole-mount in situ hybridization (WISH) showed a robust increase of gata2b and runx1 expression in the DA (Fig. 4C), suggesting that high Etv2 dosage promotes HEC fate. Consistent with the overexpression results, injection of lowdose etv2 MO (25,42) into fertilized embryos at the one-cell stage (SI Appendix, Fig. S3C) dramatically reduced gata2b and runx1 expression in the DA (Fig. 4C), showing that reduced etv2 expression blocks HEC formation. Notably, the DA formation remained intact in both Etv2 overexpression and lowdose etv2 MO knockdown embryos (Fig. 4C), indicating that low Etv2 dosage is sufficient for cEC specification. Collectively, these data indicate that etv2 expression level, high and low dosage, determines HEC and cEC fate specification from hemogenic angioblasts, respectively.
Previous studies have documented that the etv2-fli1a axis is sufficient for cEC specification (25,44) (SI Appendix, Fig. S4A). Consistent with this notion, fli1a expression was detected, although at a lower level, in low-dose etv2 MO knockdown embryos (SI Appendix, Fig. S4B), which explained the intact vascular structure in low-dose etv2 MO knockdown embryos (Fig. 4C). Yet the mechanism controlling the HEC specification is less clear. To dissect the genetic network governing HEC fate specification behind Etv2 high dosage, we monitored the Notch signaling, which has been shown as a key regulator in   HEC fate determination (20,45). We found that the expressions of notch1 receptors (notch1a/notch1b) and Notch downstream targets (hey/hey2/her12) were significantly enriched in the etv2-d2eGFP high cell fraction (Fig. 4A), suggesting that high Etv2 dosage may function through regulating Notch activity. To test this possibility, we examined the expression of notch receptors and Notch activity under etv2 overexpression and low-dose etv2 MO knockdown conditions. Results showed that the expressions of arterial notch receptors (notch1a/notch1b/notch3) (46), Notch downstream targets (hey/hey2/her12), and Tp1-eGFP reporter in the DA were significantly enhanced in etv2 overexpression embryos and dramatically reduced in low-dose etv2 MO knockdown embryos ( Fig. 4D and SI Appendix, Fig. S4 C and D). In contrast, etv2 expression, as indicated by the formation of the bilateral stripes of etv2 + cells, was largely unaffected in Delta-Notch signalingdeficient mib mutants (47) (SI Appendix, Fig. S4E), indicating that Etv2 regulates HEC specification, at least in part, through modulating Notch activity. Notably, fli1a expression was also highly enriched in etv2-d2eGFP high population (Fig. 4A), which is consistent with previous findings that fli1a is a downstream target of Etv2 (44), suggesting that Fli1a may be involved in the HEC specification. Indeed, we showed that injection of fli1a messenger RNA (mRNA) into low-dose etv2 MO knockdown embryos could partially rescue the expression of Tp1-eGFP reporter (SI Appendix, Fig. S4C), suggesting that Notch signaling is regulated directly by Etv2 as well as indirectly through Fli1a. We therefore wondered whether restoration of Notch signaling is sufficient to rescue HEC formation in low-dose etv2 MO knockdown embryos, in which fli1a was expressed at low level (SI Appendix, Fig. S4B) by overexpressing Notch intracellular domain (NICD), the dominant activator of Notch pathway. Results showed that NICD overexpression only marginally rescued HECs (Fig. 4E), suggesting that additional downstream factor(s) are required for HEC specification. We then turned to sclβ, an etv2 downstream target which has been shown indispensable for runx1 expression in the DA (25,26). We found that overexpressing Sclβ together with NICD, but not Sclβ or NICD alone, robustly rescued the HEC formation (Fig. 4E), suggesting that Sclβ and Notch are two key Etv2 downstream targets which act in parallel and collaboratively to activate gata2b-runx1 axis for HEC specification. Indeed, epistatic analysis showed that Sclβ and Notch were independent of each other as indicated by the normal expression of sclβ and Tp1-eGFP reporter in mib mutants and sclβ morphants (48), respectively (SI Appendix, Fig. S4 E and F ). Consistent with this notion, gata2b and runx1 expression were diminished in sclβ MO knockdown embryos and overexpressing sclβ drastically enhanced gata2b and runx1 expression (Fig. 4F ), placing sclβ upstream of the gata2b-runx1 axis. Taken together, these data demonstrate that the fate commitment of cECs and HECs from hemogenic angioblasts is determined by the Etv2 dosage, in which low dosage of Etv2 triggers low levels of Fli1a and Notch to direct cEC formation (25), whereas high dosage of Etv2 induces high levels of Fli1a, Notch, and Sclβ expression, leading to the promotion of HEC specification (SI Appendix, Fig. S5).

DA.
Our above analysis placed the activation of the gata2b-runx1 axis as the final effector for HEC specification. Yet whether induction of gata2b-runx1 expression is sufficient to promote HEC formation in vivo remains unclear. To address this issue, we generated a transgenic line Tg(flk1:runx1-P2a-GFP), in which Runx1-P2a-GFP expression was under control of the pan-endothelial flk1 promoter (31) (Fig. 5A), and asked whether ectopically expressing runx1 in the DA roof endothelium could promote the HEC formation. To facilitate our analysis, we outcrossed the Tg(flk1:runx1-P2a-GFP) line with the Tg(flk1:NLS-Eos) fish. The resulting double transgenic fish Tg(flk1:runx1-P2a-GFP;flk1:NLS-Eos) enabled us to label the ECs in the DA roof at 26 to 28 hpf and followed their behavior by time-lapse imaging. Consistent with previous results (Fig. 2C), only 3 out of 352 roof ECs (3/352, ∼0.85%) underwent EHT in control siblings (Tg − ; n = 32) (Fig. 5 B and C). Strikingly, in runx1 overexpression embryos (Tg + ; n = 40), 57 out of 408 roof ECs analyzed underwent EHT, a 17-fold increase compared to control siblings (∼14% versus ∼0.85%) (Fig. 5C). We further showed that the DA roofborn Eos + HECs were able to colonize the caudal hematopoietic tissue (CHT) and thymus after EHT (Fig. 5 D and E), suggesting that the roof-born HECs in runx1 overexpression fish are indeed capable of becoming hematopoietic cells. This conclusion was further supported by the findings that the Tg(flk1:runx1-P2a-GFP) fish produced significantly more myb + HSPCs than control siblings (Fig. 5 F and G). Taken together, these results demonstrate that runx1 overexpression is sufficient for the conversion, at least in part, of HECs' fate in the DA roof. Interestingly, we noticed that all 57 EHT events from the DA roof in runx1 overexpression embryos took place after they had moved ventrally and repositioned in the floor of the DA (Fig. 5D), supporting the notion that EHT occurs exclusively in the floor of the DA.

Discussion
In this study, we reported a systematic investigation on the ontogeny of cECs and HECs from early lineage fate mapping, genetic fate determination, and spatial distribution in the DA. By combining photoconversion with time-lapse imaging, we unexpectedly found that nearly the entire aortic floor endothelium are HECs, whose loss during EHT is replenished by the ventral movement and proliferation of the roof cECs in the DA. Interestingly, our above finding is similar to the previous study in the avian model that the aortic floor appears to be a transitory structure and would be completely replaced by roof and side ECs in the aorta (11), suggesting that a similar model would also be utilized in mammals.
Our study here provides in vivo evidence that cECs and HECs in the DA arise from a bipotent common progenitor, which undergoes fate separation during the time window of 14 to 18 hpf, and finally colonize the DA roof and floor accordingly. This common progenitor, as shown in our fate mapping data, is in a previously defined flk1 + angioblasts population (30,31) and distinct from the previous conception, hemangioblasts, which refer to preprimitive wave mesodermal cells giving rise to yolk sac blood islands in mouse (49) and primitive erythrocytes and ECs in zebrafish (50). We therefore term the common precursors of cECs and HECs as hemogenic angioblasts, which represent a branch point between definitive hemogenic endothelium and conventional endothelium in the DA. Notably, our common progenitor model appears different from the previous findings in avians, which showed that distinct mesoderm regions, the somatopleural and the splanchnopleural mesoderm, give rise to aortic roof and floor endothelium, respectively (6,11,12,51). The discrepancy between the avian and zebrafish models might be explained by the different methods used in the two model systems. Instead of the single-cell resolution photoconversion assay performed here, the transplantation assay employed in avian model is mainly based on anatomical recognition with lower resolution. Hence, the transplantation assay may not reflect the endogenous developmental trajectory of aortic roof and floor endothelium. Alternatively, the discrepancy is attributed to the species difference between avian and fish. Notably, several studies have suggested that, like zebrafish, the ontogeny of the DA in mice appears to arise from single mesodermal source (52,53), implying that mammals may employ a mechanism similar to that of zebrafish.
Another interesting observation in our study is that ectopic runx1 overexpression can promote HEC formation in the roof of the DA, and yet these roof-derived runx1 overexpressing cells undergo EHT only after they have repositioned to the floor of the DA. This observation indicates that in addition to the early fate determination during the time window of 14 to 18 hpf, extrinsic signals are required for the induction of EHT. One possible mechanism is that the floor niche produces extrinsic induction signals to trigger the initiation of EHT. Alternatively, the roof environment generates extrinsic inhibitory signals to Merge Eos Eos DA DA 3 dp 3 dp 3 dp 3 dp 3 dp 3 dp 3 dp 3 dp 3 dp 3 dp 3 dp 3 dp 3 dp 3 dp 3 dp 3 dp 3 dp 3 dp 3 dp 3 dp 3 dp d Tg Tg-Tg Tg-Tg Tg g g g g g g-g g-g g g-g g g g g prevent EHT in the roof of the DA. It is unclear which mechanism is correct, and perhaps both mechanisms are involved. Indeed, recent studies in zebrafish have broadened the understanding of the EHT process by showing that mechanical force generated by blood flow plays a critical role in controlling the ventral emergence and EHT budding direction of the HECs (36,37). Notably, a previous study, in which the authors utilized the whole-mount transparency method to quantify all hematopoietic clusters (c-Kit + ) temporally and spatially in E10.5 murine embryos, observed that some c-Kit + clusters were distributed along the dorsal wall of the DA in E10.5 murine embryos, while most of the c-Kit + clusters were ventrally localized (35). However, the wholemount transparency method could not provide direct evidence supporting that the c-Kit + clusters along the dorsal wall indeed arise in situ. It will be of interest to reinvestigate whether the dorsal wall of the murine DA consists of HECs. Finally, according to our data, we speculate that Runx1 expression does not seem to regulate the movements of the roof cells from the roof to the floor, since cell movements from the roof to the floor occur regardless of whether the roof cells express Runx1 or not, as indicated in runx1-overexpression or wild-type embryos, respectively. Instead, Runx1 is likely to play a role in HEC formation and EHT process (24,54,55).
Finally, we have proposed a model in which the Etv2 dosage plays a critical role in cEC and HEC fate specification (SI Appendix, Fig. S5). Previous studies identified that Fli1 (42,44), Scl (25,42), Notch1 (27), and Dll4 (56) are downstream of Etv2 in the regulation of EC development. Our study further reveals that low Etv2 dosage is sufficient to promote cEC commitment by inducing a low level of Fli1a, whereas high Etv2 dosage induces HEC lineage commitment through inducing a higher level of Fli1a and Sclβ expression and higher activity of Notch signaling. Interestingly, we also retrieved the scRNAseq gene expression matrix of four cell populations (Conflux AE, pre-HE, HE, and IAC) from a murine scRNA-seq dataset reported previously (17) and examined the expression of Etv2, Scl, Notch, Gata2, and Runx1. We observed an increase of both Etv2 expression level and the percentage of Etv2-positive cells between pre-HE and HE populations, accompanied with the increase of Runx1 expression (SI Appendix, Fig. S6 A-D). Likewise, the expression levels and positive cell proportions are also increased for Scl, Notch2, and Notch3 but not for Gata2 (SI Appendix, Fig. S6 D and E), which seems to be constitutively highly expressed in those four populations. The discrepancy of the Gata2 expression between zebrafish and murine datasets may be due to species difference since zebrafish contain duplicated gata2a and gata2b genes, and only gata2b shows HEC-specific expression (23). Nevertheless, the overall expression of the key transcriptional regulators seems to be highly conserved between these two datasets, suggesting a similar mechanism may be applied to the murine model. It will be of great interest to further investigate this issue in the murine model. Another interesting finding is that our epistatic analysis indicates Sclβ act upstream of gata2b and is required for gata2b expression in a Notch-independent manner. However, how Sclβ regulates gata2b expression remains unknown. Given the fact that the enhancers of Gata2, Fli1, and Scl are shown to be active during murine HSC ontogeny (57) and they have also been shown to form a recursively wired generegulatory circuit during early mouse hematopoietic development (58), we speculate that Sclβ may act collaboratively with Fli1a to regulate gata2b expression in zebrafish. Further studies will be required to clarify these issues.
Photoconversion Fate Mapping, Time-Lapse Imaging, Image Processing, and Unilocalization/Bilocalization Identification. For early fate mapping analysis, a single flk1-Eos + cell of Tg(flk1:eGFP;flk1:NLS-Eos) was converted by 405-nm UV laser at 15-and 18-hpf stages in each embryo as described (66). Single EC was randomly photoconverted within the 8-to 13-somite region at 15 hpf or 6-to 18-somite region at 18 hpf. With exception of three embryos, all photoconverted ECs could give rise to ECs of the DA in the AGM region (8to 17-somite region) at 26 hpf (pretested in this study). For time-lapse imaging to capture EHT, the ECs (5 to 10 cells) in the DA of Tg(flk1:eGFP;flk1:NLS-Eos) embryos were converted by 405-nm UV laser at around 26 hpf. The embryos were subsequently mounted in 1% low-melting-point agarose in E3 embryo medium with 0.01% tricaine for anesthetization as described (66). Time-lapse imaging was carried out on SP8 confocal microscope using a 20× dry objective (Leica). A heated microscope chamber at 28.5 • C was used for recording time-lapse videos. Stacks were taken every 15 min with a step size of 3.0 to 5.0 μm. Confocal stacks and time-lapse videos were analyzed using ImageJ (NIH). For unilocalization/bilocalization identification of descendant cells, Z projections of images were first performed. The roof/floor (cEC/HEC) ECs were identified according to the same criterion (0 to 45 • from the dorsal-ventral axis of the DA) relative to the roof/floor boundary indicated by the flk1-eGFP expression (37).
Heat Shock, Morpholino Injection, and mRNA Injection. Morpholinos were purchased from Gene-Tools (USA) and dissolved in distilled water. One ng/nL (high dose) or 0.1 ng/nL (low dose) of etv2 MO (42) was injected into one-cell-stage embryos. The high dose of etv2 MO completely blocked the etv2-d2eGFP expression of the Tg(etv2:NLS-d2eGFP) reporter (SI Appendix, Fig. S3C), while the low dose of etv2 MO (≤ 0.3 ng/nL) partially reduced etv2-d2eGFP expression (SI Appendix, Fig. S3C). sclβ MO knockdown was performed as previously described (26), and 8 ng of sclβ MO was injected into one-cell-stage embryos. The mRNA of fli1a and sclβ were synthesized by using mMessage mMachine Kit (Ambion Inc.). The primers for amplifying fli1a and sclβ fulllength coding sequences and injected dosage of fli1a (60 pg) and sclβ (200 pg) mRNA were referenced to a previous report (25). For heat shock experiments, embryos at the 10-somite stage (14 hpf) were collected in a 50-mL Falcon tube containing 5 mL of E3 embryo medium and submerged in a 40 • C water bath for 45 min.
WISH, Double Fluorescence In Situ Hybridization, and Immunofluorescent Antibody Staining. WISH was performed as described (66). AP-conjugated anti-DIG antibody (11093274910; Roche) and BM purple (11442074001; Roche) were used as the substrate for color reaction. Double fluorescence in situ hybridization (dFISH) assay was performed as described (67). , and eGFP (eGFP full-length coding sequence). For antibody staining, the pretreatment of embryos was similar to WISH as described above. Anti-cleaved caspase3 (559565, Biosciences, diluted at 1:500) and Alexa 555-anti-rabbit antibody (A31572, Invitrogen, diluted at 1:400) were used as primary and secondary antibody, respectively. WISH images were captured by Stereo Discovery V20 microscope (Zeiss) with digital camera (DS-Ri2, Nikon). dFISH and antibody staining images were captured by SP8 confocal microscope (Leica) and analyzed by ImageJ (NIH). FACS, cDNA Synthesis, and qPCR. The embryos at desired stages were pooled together (each biological replicate with 60 to 100 embryos) and homogenized by syringe with 27G needle. The cell suspension was prepared as previously described (68) and subjected to FACS analysis (FACS Aria IIIu). The cDNA was prepared from sorted cells ( 100 to 500 cells per well; two or four technical replicates per biological sample) and amplified by Smart-seq2 protocol (69).Realtime qPCR was performed to quantify the transcripts of gata2b, etv2, fli1a, scl, notch1a, notch1b, notch3, hey, hey2, and her12. Primers for qPCR used in this work are listed in SI Appendix, Table S1.
Single-Cell Isolation, cDNA Preparation, and RNA-seq. For single-cell isolation, 10 to 20 cells in the DA roof or floor of each Tg(flk1:eGFP;flk1:NLS-Eos) embryo were photoconverted at 21 or 28 hpf under a confocal microscope. Thirty embryos were anesthetized, and the trunk regions containing the converted cells were dissected with 27G needle in ice-cold phosphate-buffered saline (PBS). The dissected tissues were dissociated into single-cell suspension according to the protocol (70). After 20 min digestion, the single-cell suspension was pushed through a 40-μm Cell Strainer (352340, BD Falcon) and transferred to the 35-mm petri dish (150460, Invitrogen). A single red Eos + cell was collected manually with the micromanipulator system (NT-88-V3, Nikon) under a fluorescence microscope (Nikon), then washed by transferring the cell to a drop of 500 μL RNase-free PBS with 2% bovine serum albumin, and finally transferred into RNase-free 200 μL tube (Eppendorf) containing 4.4 μL lysis buffer (69). The yield and efficiency of manual single-cell collection were relatively low. In our hand, we were able to collect 20 to 30 single RNAqualified cells from 300 to 600 photoconverted cells within 2 h after digestion. The Smart-seq2 protocol (69) was used for whole-transcriptome amplification, and the cDNA library was prepared as previously described (68). Library quality was assessed using the Agilent Fragment Analyzer System, and the samples were sent to Novogene Company Limited for Illumina Hiseq X Ten 150-bp paired-end sequencing with an average depth of 6 × 10 6 raw reads per sample.
Bioinformatic Analysis of RNA-seq. Reads were aligned to the GRCz11.94 zebrafish reference genome using Spliced Transcripts Alignment to a Reference aligner (2.7.1a). Aligned reads were converted to counts for each gene using featureCounts (Rsubread 2.6.1). Dataset from 21-and 28-hpf roof cECs and floor HECs were analyzed by Seurat package (v.4.0.2). Cell quality control was performed before analysis based on minimal detected genes (default, n = 200). First, variable genes (default, n = 2,000) were identified based on variance stabilizing transformation using a normalized matrix. Principal component analysis was then performed (Dims = 1:20) followed by scaling scRNA-seq data. Finally, a plot graph was constructed by shared nearest neighbor-based clusters of the subtypes and transformed to UMAP for dimension reduction for intuitive visualization. The top 100 maker genes were identified based on spatial information (roof versus floor). GO analysis was performed with the clusterProfiler package (v.3.16.1) using the top 100 marker genes. For RNA-seq data integration analysis, GSE112294 (38) was reanalyzed by fetching counts data at 10, 14, and 18 hpf. Because etv2 expression decreases as endothelial precursors undergo differentiation, we set criteria as etv2 expression (counts > 0) at 10 and 14 hpf, flk1 expression (counts > 0) at 18 hpf, to extract potential endothelial/precursor cells, referred to as "Wagner's dataset." The Wagner's dataset was integrated with the 21-hpf dataset by the FindIntegrationAnchors function implanted in the Seurat package. Following the standard integration analysis tutorial (https://satijalab.org/seurat/articles/integration introduction.html), we set DefaultAssay as "integrated" and performed clustering analysis. All the analysis was performed in the R (version 4.0.2) environment.
Quantification and Statistical Analysis. Statistical parameters including the exact value of n are reported in Figs. 1, 2, 4, and 5 and the figure legends of Figs. 1-5. All statistical analyses were performed using GraphPad-Prism version 8. Paired or unpaired Student's t tests were used to calculate the twotailed P value for comparisons, clarified in the figure legends of Figs. 4 and 5. Data were presented as mean ± SD. Statistical significance is shown as follows: ns, P > 0.05; *P ≤ 0.05; **P ≤ 0.01; ***P ≤ 0.001; and ****P ≤ 0.0001.
Data and Software Availability. The raw fastq files for each sample reported in this paper and single-cell RNA-seq data are available in the NIH Gene Expression Omnibus database (accession no. GSE197757).