Clinical details for APECED patients were collected from medical records, through a questionnaire and patient interviews, which included questions on medical history, prior infections, medications, and other relevant parameters. All patients were clinically examined by a medical doctor (S.L.). The grading of the severity of GI symptoms experienced by the patients was done based on discussion and a unanimous decision by two physicians (I.H. and S.L.). The severity of diarrhea was classified according to the degree of symptoms experienced: 1) no symptoms, 2) intermittent diarrhea, episodes of diarrhea separated by obstipation or symptom-free periods, 3) chronic difficult diarrhea.

Fecal samples for the patients and controls were collected at home and immediately stored at -20°C. They were transported to a study center within 1 week. An uninterrupted frozen cold chain was ensured in the provision and handling of the fecal samples. Bacterial DNA was extracted from ca. 250 mg of fecal matter using the Repeated Bead Beating (RBB) method (13) with the following modifications for automated DNA purification: 340 μl and 145 μl of lysis buffer was added to first and second round of bead beating, respectively. 200 µl of the clarified supernatant collected from the two bead beating rounds was used for DNA extraction with the Ambion Magmax™ -96 DNA Multi-Sample Kit (Thermo Fisher Scientific, USA) using the KingFisherTM Flex automated purification system (Thermo Fisher Scientific, USA). DNA was quantified using Quanti-iT™ Pico Green dsDNA Assay (Invitrogen, San Diego, CA, USA).

Blood samples were collected in the morning between 7 and 10 am after an 8- to 12-hour fast before taking morning medications. Sera were isolated with a standard protocol and stored at -80°C until analyses.

The ASCA were quantified from patient sera with the Anti-Saccharomyces cerevisiae Antibodies IgG test kit (Bio-Rad, USA) as instructed by the manufacturer. Values over 15 U/mL were considered positive.

Immunoprecipitation of radiolabeled antigens was used to screen sera for autoantibodies against tryptophan hydroxylase-1. Human complementary DNA from the antigen in expression vector was used to perform the assay. In vitro transcription and translation were performed in the presence of ³⁵S-methionine, according to the manufacturer’s protocol (Promega TNT Systems). Immunoprecipitation was performed in 96-well plates overnight at 4°C at 300 rpm with serum samples (2.5 µl) and 30,000 cpm of radiolabeled protein. A positive control (positive patient serum to each antigen) and a negative control, 4% bovine serum albumin, were included in each plate. All samples were analyzed in duplicates. The immune reaction was transferred to filter plates (Millipore) and immune complexes were captured to Protein A sepharose (nProtein A Sepharose 4 fast flow, GE Healthcare) during a 45 min incubation at 4°C at 300 rpm. After ten washing steps with wash buffer (150mM NaCl, 20 mM TrisHCl pH 8, 0.15% Tween 20, and 0.1% BSA), plates were dried and scintillation fluid (Optiphase, HiSafe 3, PerkinElmer) was added. Radioactivity was then measured in a beta counter (1450 Microbeta Trilux, Wallac). Autoantibody index values were calculated according to the following: (sample value – negative control value)/(positive control value – negative control value) x 100.

For qPCR targeting fungal ITS1-region, 5 ng of DNA was mixed with 2 μM of ITS1F and ITS2 primers (14) and 2X Power SYBR green PCR Master Mix (Thermo Fisher, USA) in a total volume of 20 μl. The PCR reaction was performed in the C1000 Touch Thermal Cycler (Bio-Rad, USA) as follows: 95°C for 15 min, followed by 40 cycles of 95°C 15 sec, 55°C 30 sec and 72°C 30 sec and ending in 95°C 1 min and 60°C 1 min. Results were analyzed using the comparative Ct method (15).

Library preparation and Illumina MiSeq sequencing of the hypervariable V3-V4 regions of the 16S rRNA gene and sequencing data preprocessing for the patient and control samples as well as negative controls were performed as previously described (16).

The pre-processing of the sequencing reads, their taxonomic annotation, and statistical analysis were performed using R package mare (Microbiota Analysis in R Easily) (17) as previously described (16, 18). Samples from cases and controls were sequenced in two different MiSeq runs with non-template negative controls. The lowest number of reads per sample was 28175 and the mean number of reads was 72610 for patients and 73661 for controls. The three negative control samples in the run with APECED samples contained 139, 259, and 728 reads assigned to two uncultured lactobacilli. As the reads assigned to Lactobacillus in the APECED samples were also mostly assigned to the two uncultured lactobacilli in negative controls, we excluded Lactobacillus from all data analyses. To account for the varying sequencing depth, the number of reads per sample was used as an offset in all statistical models. The β-diversity was estimated using Bray-Curtis dissimilarity as the distance measure and the contribution of different variables to microbiota variation was calculated using permutational ANOVA (PERMANOVA) adonis function in vegan (19) on genus, family and phylum levels. The BMI or gender were not associated with microbiota variation (permutational ANOVA 4% and 3%, respectively, FDR-P > 0.3 for both) and were not used as confounding factors. Also, age did not have a statistically significant effect on permutational ANOVA, but it correlated with the relative abundance of some of bacteria of interest in patients (Faecalibacterium, Megasphaera and Prevotella) so it was used as a confounder in all models. The differences in the microbiota between patients and controls or between patients with different clinical characteristics were analyzed using generalized linear models with negative binomial distribution implemented in the mare package for prevalent bacteria i.e. bacterial genera detected in > 60% of the samples. P-values were adjusted by the Benjamini-Hochberg method for multiple testing and reported as FDR-P. All significant findings from differential abundance testing underwent visual inspection when applicable to eliminate statistically significant findings driven by few extreme values. Only statistically significant and visually validated results are reported. In the univariate data, a statistical difference was evaluated using the Independent-Samples Mann-Whitney U test in SPSS 25 (IBM). For all tests, P-values < 0.05 were considered as statistically significant.

Bacterial metagenome content was predicted from the 16S rRNA gene-based microbial compositions. Functional inferences were made from the Kyoto Encyclopedia of Gene and Genomes (KEGG) catalog (20) using PICRUSt2 (21). The “DESeq” function in DESeq2 (22) was used to test for differentially abundant KEGG pathways between the groups.

Random forest model was used to evaluate the predictive performance of all identified genera (N = 65) or predicted functions (N = 222) based on leave-one-out cross-validation. The random forest analysis was performed using the default settings of the “randomForest” function implemented in the randomForest R package. We then used the receiver operating characteristic (ROC) curve and calculated the area under the ROC curve (AUC), implemented in the pROC package, as the measure of performance. The importance of input features was evaluated by mean decrease Gini values.

The microbiome samples of altogether 15 patients (7 females) with APECED were investigated. The patients were examined at the median age of 46.4 years (range, 19.3 – 70.1). Characteristics of the patients are presented in Table 1 in comparison to 22 patients with APECED from whom we were unable to acquire stool samples. The patient characteristics differ little between the two groups suggesting that the sample of patients participating in the study represents well the overall APECED population.

All patients harbored at least one copy of the Finnish founder mutation c.769C>T, p.(Arg257Ter) in the AIRE gene (NM_000383.4). Altogether 13/15 of patients had primary adrenal insufficiency, and 12/15 had hypoparathyroidism. The patients had median of 7 (range 4-11) disease components of APECED. Additional symptoms, each present in one individual, included high blood pressure, asthma, depression, migraine, APC resistance, hypercholesterolemia and mild hypertrophic cardiomyopathy. All 15 patients had experienced mucocutaneus candidiasis. One patient was using anti-fungal medication with treatment dose and three patients with prophylaxis dose. Either periodic or chronic diarrhea was evident in 60% of the patients. None of the patients had used antibiotics within 1 month prior to the sampling. Four patients were using probiotics at the time of sampling, but these patients did not significantly differ from other patients in their microbial composition (data not shown). None of the controls had used antibiotics within 3 months prior to the sampling and none of them used probiotics.

We performed 16S rRNA gene amplicon sequencing to investigate the microbial composition of the stool samples. After processing we had on average 73135 (28175-164233) high quality reads representing overall of 395 operational taxonomic units (OTUs) distributed to 78 bacterial genera, 33 families, 19 orders, 13 classes, and 5 phyla. The microbiota richness ranged from 94 to 220 (mean, 168 vs. 152 in patients and controls, respectively; not significant; Mann-Whitney U test). Alpha diversity was significantly higher in patients compared to healthy controls [mean 13.0 (range, 4.2-18.4) vs. 6.8 (4.3-10.5); P < 0.001; Mann Whitney U test].

We compared the difference between patients and controls in the overall gut microbiota using PERMANOVA based on Bray-Curtis dissimilarity. Patients differed from controls already on the phylum level (14% of variation explained, P = 0.006, Figure 1A ). Bacteroidetes and Proteobacteria were expanded in patients’ gut microbiota, while a lower proportion of Firmicutes was found ( Figure 1B ) compared to the controls. On the family level, PERMANOVA indicated a 10% variation in the gut microbiota between the patients and controls (P = 0.003). Two patients distinctly differed from others with markedly expanded relative abundance of Veillonellaceae and diminished Ruminococcaceae ( Figure 1C ). Interestingly, these subjects were two of the three patients reporting the most severe diarrhea.

Beta-diversity between the patients and controls was significantly different also on genus level (10% variation in PERMANOVA, P = 0.001). In comparison to healthy controls, the proportion of Faecalibacterium was significantly reduced in patients [fold change (fc) = 0.45, FDR-P = 0.05, Table 2 , Figure 2 ]. Multiple genera were significantly overrepresented in the patients. The largest 15-20 fold increase was found in the fraction of Butyrivibrio, Atopobium, and Haemophilus (fc = 15.1, FDR-P = 7x10^-5; fc = 18.8, FDR-P = 7x10^-5; fc = 17.3, FDR-P = 0.0003 respectively, Table 2 and Figure 2 ). Atopobium was almost undetectable in the majority of the controls while it was observed in all patients accounting up to 20.3% of the total microbiota in one patient. The relative abundance of other gram-negative bacteria, Veillonella and Prevotella, was significantly increased in patients (fc = 5.4, FDR-P = 0.02; fc = 2.0, FDR-P = 0.04; respectively, Table 2 and Figure 2 ). Of interest, the reads derived from Megasphaera spp. within the Veillonnellaceae family accounted for up to 37.1% and 0.8% of reads in two patients, respectively, but was otherwise virtually undetectable in the majority of patients and controls ( Figure 2 ). All genera with relative abundances that are significantly different between the patients and controls are summarized in Table 2 .

To understand functional implications of the observed taxonomic difference between patients with APECED and controls, we inferred metagenomes using the PICRUSt2 algorithm. This is a computational approach that reconstructs functional composition of a metagenome connecting the sequenced genes to reference genomes. The imputed functions of patients’ microbiota differed significantly from those of healthy controls, explaining 9% of the variation in the gut microbiota (P = 0.02, Figure 3A ), which was attributable to the functional modules associated with LPS synthesis that were enriched in patients ( Figure 3B ).

Random forest (RF) classification models were subsequently trained by 5-fold cross-validation using the relative abundances of all genera (N = 65) or imputed functional modules (N = 211). Both RF models including bacterial taxa or predicted functions classified patients with APECED with a comparably robust accuracy, achieving an area under the curve (AUC) of 0.83 and 0.87, respectively ( Figure 3C ). The most important genera selected by the model included Faecalibacterium and Atopobium ( Supplementary Figure 1A ). The functional modules related to LPS synthesis were among the most important features selected by the random forest model ( Supplementary Figure 1B ).

Having established the differences in the gut microbiota between patients and controls at the level of the entire bacterial community and individual taxa, we next examined whether patient characteristics were linked to gut microbiota deviations. We used the ASCA quantification that was available for all but one of the patients with APECED, and 8/14 (57.1%) were positive, indicating a systemic immune response to fungal products or fungi. Five of the ASCA positive patients had diarrhea including the three patients with most severe symptoms. The ASCA positivity explained 10% of variation in the gut microbiota on the phylum level (P = 0.05, Figure 4A ). As mentioned previously, a lower proportion of Firmicutes was characteristic to patients compared with healthy controls, which was more pronounced in ASCA positive patients (FDR-P = 1.1x10^-6 between ASCA+ and ASCA- patients, Figure 4B ).

Interestingly, the antigen for ASCA antibodies can also be produced by C. albicans and infection by this pathogenic yeast has been shown to be able to induce ASCA (23). Patients with APECED are susceptible to chronic candidiasis (1, 2) and all of our patients had suffered from chronic candidiasis at some site of the body with varying severity. We used ITS1 primers to quantify fungi in patient stool microbiota by qPCR. Fungal DNA was amplified in all patient samples and based on the amplification results, the patients were grouped into three categories using a comparative Ct method (median Cq of patients +/- 1 cycles, respectively). Three of the patients had a markedly higher level of fungi in their stool, while two patients had a significantly lower level than other patients. Bacterial composition of these groups differed significantly (21% variation explained in PERMANOVA, P = 0.04, Figure 4C ). Four patients were receiving anti-fungal medication at the time of sampling (one in the low group, the rest in the median group).

Autoimmunity against neuroendocrine cells as measured by TPH antibodies was detected in 60% of the patients (9/15) and was associated with a higher abundance of Haemophilus (FDR-P = 0.001, fc = 65 to patients without anti-TPH antibodies).

We examined the associations between the intestinal microbiota and GI symptoms of the patients. Four out of the fifteen patients (26.7%) suffered from intermitting obstipation episodes, but their stool microbiota did not differ from the other patients. The opposite was, however, true for patients reporting symptoms of diarrhea. We divided patients into three groups according to the severity of their symptoms of diarrhea (no symptoms, N = 6; intermittent diarrhea, N = 6; chronic difficult diarrhea, N = 3) and found the fecal microbiota composition to significantly differ among the groups (24% of variation explained in PERMANOVA on the family level, P = 0.02; 15% variation explained in PERMANOVA on the genus level, P = 0.02, Figure 4C ). Atopobium spp. were overrepresented in patients with intermitted diarrhea and chronic difficult diarrhea (fc = 6.3, FDR-P = 2x10^-35; fc = 117, FDR-P = 6x10^-305, respectively, compared to patients without symptoms) while Faecalibacterium spp. were less common if a patient had chronic difficult diarrhea (fc = 0.3, FDR-P = 0.05 compared to patients without symptoms).

The microbial composition of patients without any GI symptoms (N = 5), also differed from the healthy controls (11% of variation explained in PERMANOVA on the genus level, P = 0.04).