Shotgun metagenomics reveals signi cant gut microbiome features in different grades of acute pancreatitis

Shanshan Yu Peking Union Medical College Hospital&Peking Union Medical College&Chinese Academy of Medical Science and Peking Union Medical College Yangyang Xiong Peking Union Medical College Hospital&Peking Union Medical College&Chinese Academy of Medical Science and Peking Union Medical College Yangyang Fu Peking Union Medical College Hospital&Peking Union Medical College&Chinese Academy of Medical Science and Peking Union Medical College Guorong Chen Peking Union Medical College Hospital&Peking Union Medical College&Chinese Academy of Medical Science and Peking Union Medical College Huadong Zhu Peking Union Medical College Hospital&Peking Union Medical College&Chinese Academy of Medical Science and Peking Union Medical College Xun Mo The Second People's Hospital of Guiyang City Jun Xu Peking Union Medical College Hospital&Peking Union Medical College&Chinese Academy of Medical Science and Peking Union Medical College Dong Wu (  wudong@pumch.cn ) https://orcid.org/0000-0001-6439-3761

Shotgun metagenomic sequencing of the whole DNA provides valuable information about the functions of the microbial community [10]. The sequencing technology is used to obtain the entire genomic content of the microbiome and achieve accurate taxonomic classi cation and functional assignments. Also, it is able identify novel functional genes, antibiotic resistance genes, microbial pathways, and functional dysbiosis of the gut microbiome [11]. In this study we conducted metagenomic shotgun survey on microbial gut population of MAP, MSAP, SAP patients and three healthy controls, in order to characterize the AP-related composition and functional changes.

Subjects recruitment
The 2012 revised Atlanta classi cation strati ed the clinical severity of AP patients into three categories: mild AP (MAP), moderately severe AP (MSAP), and severe AP (SAP) [3]. All the patients with a diagnosis of AP based on the 2012 Revised Atlanta Classi cation and admitted to Peking Union Medical College Hospital, Beijing, China were eligible for inclusion if they were enrolled within 48 hours of the onset of symptoms. A total of three MAP, three MSAP, and three SAP patients were included. None of AP patients had been treated in other hospitals before admission.

Healthy control population
The healthy control population (HCR) were three healthy volunteers (community workers and students) identi ed from Peking Union Medical College Hospital. They were age, gender and BMI matched to the patient population and selected after the completion of the patient study. We used the same de nition of the HCR group with our previous work [9].
Exclusion criteria AP patients and HCR were excluded if they took antibiotics, probiotics, Chinese herbs, steroids, and other substances that may affect the structure of the ora in up to eight weeks prior to sampling. Other exclusion criteria referred to patients with a history of immunode ciency, allergy, asthma, celiac disease, colorectal cancer, diabetes, HIV, in ammatory bowel disease, irritable bowel syndrome, gastroenteritis, necrotizing enterocolitis, and arthritis [7].
Written informed consent was obtained from each participant. The institutional review boards of Peking Union Medical College Hospital approved this study (No. JS 1826).

Sample collection
Stool samples are generally considered the reference standard. However, they are not practical to obtain in patients with AP, since they are often a icted by reduced gastrointestinal motility [12]. Thus, rectal swabs are commonly used instead, as previously described [13]. In this study, microbial DNA was extracted from fecal samples obtained through the rectal swab immediately after admission. All individuals should have an empty stomach (did not eat for 12 hours) before sample collection. The sample of healthy volunteers was collected in the morning. The timing of the rectal swab sampling of AP patients was within 24 hours after disease onset in all individuals. Detailed processes of sample collection were described elsewhere [9]. DNA extraction, library construction and metagenomics sequencing DNA for metagenomics was extracted from fecal samples by using the E.Z.N.A.® DNA Kit (Omega Bio-tek, Norcross, GA, U.S.) according to the manufacturer's protocols. The DNA concentration and purity were quanti ed with TBS-380 and NanoDrop2000, respectively. DNA quality was examined with the 1% agarose gels electrophoresis system. DNA was fragmented to an average size of about 300 bp using Covaris M220 (Gene Company Limited, China) for paired-end library construction. The paired-end library was prepared by using TruSeqTM DNA Sample Prep Kit (Illumina, San Diego, CA, USA). Adapters containing the full complement of sequencing primer hybridization sites were ligated to the Blunt-end fragments. Paired-end sequencing was performed on Illumina HiSeq4000 platform (Illumina Inc., San Diego, CA, USA) at Majorbio Bio-Pharm Technology Co., Ltd. (Shanghai, China) using HiSeq 3000/4000 PE Cluster Kit and HiSeq 3000/4000 SBS Kits according to the manufacturer's instructions (www.illumina.com). All the raw metagenomics datasets have been deposited into NCBI Sequence Read Achieve database.
Sequence quality control and genome assembly 3' and 5' end was stripped using SeqPrep (https://github.com/jstjohn/SeqPrep). Low-quality reads (length<50 bp or with a quality value <20 or having N bases) were removed by Sickle (https://github.com/najoshi/sickle). Reads were aligned to the human genome (version 38) by BWA (http://bio-bwa.sourceforge.net) and any hit associated with the reads and their mated reads were removed. De bruijn-graph-based assembler SOAPdenovo (http://soap.genomics.org.cn) (version 1.06) was used to assemble short reads. K-mers, varying from 1/3~2/3 of reads length were tested for each sample. Scaffolds with a length over 500bp were retained for statistical tests; we evaluated the quality and quantity of scaffolds generated by each assembly and nally chose the best K-mer which yielded the minimum scaffold number and the maximum value of N50 and N90. Then, scaffolds with a length over 500 bp were extracted and broken into contigs without gaps. Contigs were used for further gene prediction and annotation. Open reading frames (ORFs) from each metagenomic sample were predicted using MetaGene (http://metagene.cb.k.u-tokyo.ac.jp/). The predicted ORFs with length being or over 100 bp were retrieved and translated to amino acid sequences using the NCBI translation table (http://www.ncbi.nlm.nih.gov/Taxonomy/taxonomyhome.html/index.cgi?chapter=tgencodes SG1). All sequences from gene sets with a 95% sequence identity (90% coverage) were clustered as the nonredundant gene catalog by the CD-HIT (http://www.bioinformatics.org/cd-hit/). Reads after quality control were mapped to the representative genes with 95% identities using SOAPaligner (http://soap.genomics.org.cn/), and gene abundance in each sample was evaluated.

Differences analysis in microbiota composition
The linear discriminant analysis (LDA) Effect Size (LEfSe) was used to support high-dimensional group comparisons to nd biological relevant characteristics. First, the non-parametric factorial kruskal-wallis sum-rank test was used to detect features with signi cant differential abundance between groups, with biological consistency investigated among subgroups using the Wilcoxon rank-sum test. Finally, the above-obtained data were used for LDA model analysis to rank the relative abundance difference with respect to groups. The signi cance for each test was set to 0.05.

Analysis of microbiome functional composition
The Venn diagram can be used to count the number of common and unique functions in multiple groups of samples. The relationship diagram of Circos samples and function is a visual circle diagram describing the corresponding relationship between samples and function, which not only re ects the composition ratio of dominant function of each sample, but also re ects the distribution ratio of each dominant function among different samples. In MAP, MSAP, and SAP, Venn diagram and Circos diagram were respectively used to analyze the functional composition of microbiota and the relationship between samples and microbiota functional composition by using R language.
Species annotation of microbiome in the amino acid sequence of non-redundant (NR) proteins database The NR gene set was compared with the NR database using DIAMOND (http://ab.inf.unituebingen.de/software/diamond/) software to obtain through taxonomic information in the NR database. The species abundance was calculated and counted in various taxonomic level including the Domain, Kingdom, Phylum, Class, Order, Family, Genus and Species to build corresponding taxonomic level abundance of spectrum (abundance pro le).
Functional annotation of microbiome in the clusters of orthologous groups of proteins (COG) database BLASTP (http://blast.ncbi.nlm.nih.gov/Blast.cgi) (Version 2.2.28+) was employed for taxonomic annotations by aligning non-redundant gene catalogs against NCBI NR database with an e-value cutoff of 1e-5. Cluster of orthologous groups of proteins (COG) for the ORFs annotation was performed using BLASTP against eggNOG database (v4.5) with an e-value cutoff of 1e-5.

Functional annotation of microbiome in the Antibiotic Resistance Genes Database (ARDB)
The ARDB (http://ardb.cbcb.umd.edu/) database includes bacterial resistance genes from different environmental sources (such as the intestinal tract), and their resistance spectra, mechanism of action, ontology and other annotation information, providing research basis for the study of drug action. Functional annotation was conducted using BLASTP search (Version 2.2.28+) against ARDB database (http://ardb.cbcb.umd.edu/) with an e-value cutoff of 1e-5.
Functional annotation of microbiome in the comprehensive antibiotic resistance database (CARD) CARD (http://arpcard. Mcmaster.ca) database contains a wide range of antibiotic resistance-related reference genes from a variety of organisms, genomes and plasmids, which can be used to guide the study of antibiotic resistance mechanisms in human and animal ora groups. Functional annotation was conducted using BLASTP search (Version 2.2.28+) against CARD database (http://ardb.cbcb.umd.edu/) with an e-value cutoff of 1e-5.

Analysis of microbiotafunctional composition difference
Based on the functional composition abundance data, relevant analysis methods were used to detect the functional composition diversity of microbial communities. In this study, the signi cance test was used to screen different functional compositions of microbiota in MAP, MSAP, and SAP. Detailed screening process was as follows: 1) Wilcoxon rank-sum test was used to test the average rank of the two groups of samples and determine whether there were differences in the distribution of the two samples; 2) multiple inspection and correction method for P value (FDR) was applied; 3) double tail test was used to specify the type of con dence interval desired; 4) Welch T test was used to calculate the con dence interval (con dence: 0.95); 5) The screening criteria is p<0.05.

Results
Clinical characteristics of patients and healthy subjects A total of 9 AP patients (3 MAP, 3 MSAP, and 3 SAP) and 3 HCR were enrolled in this study. The HCR group were all women and they were 42.0±14.2 years old with a BMI of 20.9±2.3. The clinical characteristics of AP groups were shown in Table 1. Age (p=0.473) and gender (p=0.99) were comparable between the AP groups and the HCR group. Among the AP patients statistically signi cant differences were observed in local complications, duration of organ failure, ICU admission, length of ICU stay, and length of hospital stay ( Table 1). None of the AP patients died.

Differences analysis in microbiota composition
In MAP, the taxa of Thermoprotei and Crenarchaeota was signi cantly more enriched compared to HCR (supplementary Figure 1A). Sulfolobus was remarkably more abundant in MSAP compared to HCR (supplementary Figure 1B). In SAP, enrichment in the taxa of Sulfolobales was signi cantly observed compared to HCR (supplementary Figure 1C). Methanobrevibacter_ruminantium and Methanosarcina_thermophila were uniquely abundant in MSAP compared to MAP (supplementary Figure  1D). Methanomicrobiales_archaeon_53_19 was remarkably abundant in SAP compared to MSAP (supplementary Figure 1E).

Analysis of microbiome functional composition
The Venn diagram showed that KEGG was the most enriched functional composition in MAP (supplementary Figure 2A). CAZy and VFDB were the most enriched functional compositions in MSAP (supplementary Figure 2B). From the Circos diagram, we found that macrolide was dominant functional composition of ARDB in MAP (supplementary Figure 3A); Escherichia-coli was the dominant microbiota specie of NR in MSAP (supplementary Figure 3B); In SAP, GT4 (glucosyltransferase) was a dominant functional composition of CAZy (supplementary Figure 3C).

Difference analysis of microbiotafunctional composition of NR
In the analysis of NR functional composition, there were 364, 181, and 133 signi cantly different microbiota species in MAP, MSAP, and SAP, respectively. Streptococcus and Anaerostipes hadrus was respectively the most signi cantly increased (p=0.01211) and decreased (p=0.0009861) microbiota species in MAP (supplementary Figure 4A). Unclassi ed Bacteria and Anaerostipes hadrus was respectively the most signi cantly increased (p=0.01912) and decreased (p=0.001623) microbiota species in MSAP ( Figure 4B). In SAP, Enterococcus and Blautia was respectively the most signi cantly increased (p=0.03043) and decreased (p=0.001954) microbiota species ( Figure 4C).

Difference analysis of microbiotafunctional composition of COG
In the analysis of COG functional composition, there were 0, 5, and 3 signi cantly different functional compositions in MAP, MSAP, and SAP, respectively. U (intracellular tra cking, secretion, and vesicular transport) and Z (cytoskeleton) was respectively the most signi cantly increased (p=0.002019) and decreased (p=0.02135) COG functional composition in MSAP ( Figure 1A). In SAP, C (energy production and conversion) and T (signal transduction mechanisms) was respectively the most signi cantly increased (p=0.03573) and decreased (p=0.01441) COG functional composition ( Figure 1B).

Difference analysis of microbiotafunctional composition of KEGG
In the analysis of KEGG functional composition, there were 21, 37, and 27 signi cantly different functional compositions in MAP, MSAP, and SAP, respectively. Glutathione metabolism and Neomycin, kanamycin and gentamicin biosynthesis was respectively the most signi cantly increased (p=0.0005001) and decreased (p=0.004317) KEGG functional composition in MAP (Figure 2A). Lipopolysaccharide biosynthesis and starch and sucrose metabolism was respectively the most signi cantly increased (p=0.006658) and decreased (p=0.0007967) KEGG functional composition in MSAP ( Figure 2B). In SAP, valine, leucine and isoleucine degradation and fatty acid metabolism was respectively the most signi cantly increased (p=0.00221) and decreased (p=0.0002432) KEGG functional composition ( Figure 2C).

Difference analysis of microbiotafunctional composition of ARDB
In the analysis of ARDB functional composition, there were 2, 22, and 2 signi cantly different functional compositions in MAP, MSAP, and SAP, respectively. Gentamicin was the most signi cantly increased (p=0.005638) ARDB functional composition in MAP ( Figure 4A). Penicillin and bacitracin was respectively the most signi cantly increased (p=0.0006484) and decreased (p=0.02834) ARDB functional composition in MSAP ( Figure 4B). In SAP, tetracycline was the most signi cantly decreased (p=0.009731) ARDB functional composition ( Figure 4C).  Figure 6A). Capsular polysaccharide and Two-component system was respectively the most signi cantly increased (p=0.0002623) and decreased (p=0.004764) VFDB functional composition in MSAP ( Figure 6B). In SAP, alginate biosynthesis and Twocomponent system were respectively the most signi cantly increased (p=0.003317) and decreased (p=0.0008046) VFDB functional composition ( Figure 6C).

Discussion
In this study, we performed the Shotgun metagenomic sequencing approach to a cohort of 12 individuals, including 9 AP patients (3 MAP, 3 MSAP, and 3 SAP) and 3 HCR. We found remarkably different outcomes (table 1) and signi cant dysbiosis of microbiome composition and function in the AP patients. Moreover, the composition and function of microbiota in SAP were different from MAP and MSAP, implying the association of the dysbiosis of microbial composition and function to AP severity.
Firstly, in differences analysis in microbiota composition, we found that Thermoprotei and Crenarchaeota were signi cantly enriched in MAP; and Sulfolobus was highly abundant in MSAP and SAP. We also found that Methanobrevibacter_ruminantium and Methanosarcina_thermophila were uniquely abundant in MSAP compared to MAP, and Methanomicrobiales_archaeon_53_19 was highly abundant in SAP compared to MSAP. Methanobrevibacter is a major intestinal methanogenic archaeon [14]. Gaci N et al. reported that the mammalian gut is occupied by methanogenic archeaea such as methanomicrobiales [15]. The methanomicrobiales was also detected in fecal samples in foals [16]. This was the rst time for microbiotas to be identi ed in AP patients. Further mechanism research is needed.
In the differential analysis of microbiota composition, we found that Streptococcus and Anaerostipes hadrus were respectively the most signi cantly increased and decreased microbiota species in MAP; Anaerostipes hadrus was the most signi cantly decreased microbiota species in MSAP. In addition, from the Circos diagram, we found that Escherichia-coli was the dominant microbiota specie in MSAP. Streptococcus, a type of Gram-positive bacteria, is a common pathogen in infected pancreatitis necrosis [17]. It has also been reported that Streptococcus is related to chronic alcoholic pancreatitis and pancreatic cancer [18]. Anaerostipes hadrus converts indigestible carbohydrates into fermentation products, including short-chain fatty acids (SCFAs, such as butyrate) [19]. Since SCFA is essential in maintaining gut barrier, host metabolism and immunity, the reduction of SCFAs-producing bacteria can be detrimental in patients with AP. It is not surprising that Anaerostipes hadrus is associated with the pathogenesis of chronic alcoholic pancreatitis [20]. Escherichia-coli is the predominant microbiotas in pancreatic cyst uids and consists a common pathogen of infected pancreatic necrosis in non-mild AP [21]. Our ndings suggested that the alteration of microbiotas may play roles in the development of MAP and MSAP.
In SAP, we found that Enterococcus and Blautia were respectively the most signi cantly increased and decreased microbiota species. Previous studies have demonstrated that Enterococcus could adhere and invade the surface of host cells, cross host epithelial barriers, and get access to other organs and the systemic circulation, leading to sepsis and organ failure [7]. The relative abundance of Enterococcus (potential harmful pathogens) is related to infection and systemic in ammation in patients with AP [22].
The inoculation of Il10-/mice with Enterococcus species causes characteristic chronic in ammatory diseases, supporting its pathogenic role [23]. It is also found that the population of intestinal Enterococcus is higher in the gut ora of SAP patients [24]. Blautia is the predominant taxa (bene cent microbiota) in the intestinal tract of the healthy population [25]. Signi cant reduction of Blautia levels in patients with SAP, as revealed in this study, promotes the overgrowth of intestinal bacteria, increases intestinal permeability, and ultimately results in higher concentrations of endotoxins and activation of in ammatory cascades [9] .
These ndings altogether indicated that the disbalance of Enterococcus and Blautia might play a key role in the in ammatory process of SAP.
Besides the dysbiosis of microbiome composition, we also found microbial dysfunction in AP patients. In differential analysis, we revealed that intracellular tra cking, secretion, and vesicular transport and cytoskeleton was the most signi cantly altered COG functional composition in MSAP. It is reported that disruption of polarized intracellular tra cking may leads to interstitial leakage of pancreatic enzymes and other macromolecules, and nally pancreatitis [26]. In the acinar cell, unde ned intracellular tra cking defects and secretory blockade are believed to culminate in pathologic proenzyme activation during AP [27]. AP induced by supramaximal secretagogue stimulation is characterized by the early and rapid disruption of the apical actin cytoskeleton [28]. Our result is consistent with that of Jungermann J and Ueda T et al., and they found that the tubulin cytoskeleton was disrupted in an experimental model of AP [29,30]. In SAP, we found that energy production and conversion was the most signi cantly increased COG functional composition. The gut microbiota plays a crucial role in energy metabolism to maintain the host homeostasis [31]. On account of acute systemic in ammatory response syndrome, AP (especially SAP) always triggers a hypercatabolic state, resulting in increased energy requirements [32]. Our result further suggested the important roles of intracellular tra cking, secretion, cytoskeleton, and energy metabolism in the process of MSAP and SAP.
In differential analysis, we found that glutathione metabolism was the most signi cantly increased KEGG functional composition in MAP. Glutathione system assists certain pathogens' survival in in ammatory tissues [33]. Previous studies showed that changes in glutathione metabolism occurred at the early stage of acute pancreatitis [34]. In addition, we found that lipopolysaccharide biosynthesis and starch and sucrose metabolism was the most signi cantly increased and decreased KEGG functional composition in MSAP. Lipopolysaccharide damages the intestinal barrier, which in turn increases mucosal permeability and promotes bacterial translocation [35]. Lipopolysaccharide also induces early damage and remarkable inflammation in the pancreas [36]. Raetz CR and Exley AR et al. detected the existence of lipopolysaccharides in the plasma of patients with SAP at an early stage of the disease [37,38]. In healthy individuals, the starch that entered the large bowel is metabolized by various saccharolytic bacteria (inhabiting that region of the intestine), which in turn promote the production of SCFA [39]. In patients with SAP, we found that valine, leucine, and isoleucine degradation and fatty acid metabolism was the most signi cantly increased and decreased KEGG functional composition. Valine can be used to accurately discriminate hyperlipidemic AP from healthy controls [40]. The level of valine is signi cantly decreased in the serum of chronic pancreatitis patients [41]. Some authors believed that valine could be a potential biomarker for early-stage AP [42]. Leucine and isoleucine has been found in the serum of patients with chronic pancreatitis [41]. In pancreatitis, excess free fatty acids cause oxidative stress, free radical accumulation, and acinar necrosis [43]. Saharia P and Halangk W et al. reported that perfusion of mice pancreas with fatty acid induced pancreatic edema could activate trypsinogen, which ultimately initiate AP [44,45]. Of note is that high amounts of unsaturated fatty acids have been detected in the serum of SAP patients [46]. Our study indicated that the above microbiota functional compositions of KEGG might be involved in the pathogenesis and progression of AP.
Furthermore, we found that GT41 (glycosidyltransferase) was the most signi cantly increased CAZy functional composition in MAP. In addition, from the analysis of microbiome functional composition (Circos diagram), we found that GT4 (glucosyltransferase) was a dominant functional composition of CAZy in SAP. It is reported that the glucosyltransferase activity of the toxin impacts secretion of in ammatory mediators [47]. This suggested that glycosidyltransferase may play a key role in the in ammatory response of MAP and SAP. We also showed that gentamicin was the most signi cantly increased ARDB functional composition in MAP. From the analysis of microbiome functional composition (Circos diagram), we found that macrolide was a dominant functional composition of ARDB in MAP.
Gentamicin is a macrolide prescribed for infected necrosis in AP patients [48]. Macrolides have a number of biological effects potentially bene cial for AP, including antibacterial and anti-in ammatory activity [49,50]. In MSAP, we found that penicillin was the most signi cantly increased ARDB functional composition.
Penicillin is considered a Class IV drug that could induce AP [51]. It has been demonstrated that penicillin is ineffective in treating patients with AP [52]. In addition, we found that tetracycline was the most signi cantly decreased ARDB functional composition in SAP. In fact, tetracycline has been claimed to increase the risk of AP and is ineffective in treating AP [52,53]. In this study, we found the microbial resistance in MAP, MSAP, and SAP, which suggested that gentamicin, macrolide and tetracycline would not be considered as optimal therapeutic drugs for infectious complications of AP.
In the differential analysis of microbiota functional composition of CARD, we found that aminocoumarin antibiotic and glycopeptide antibiotic was the most signi cantly increased and decreased CARD functional composition in MAP. 3-Aminocoumarin derivatives possess a number of biological activities, such as antibacterial. The N-glycopeptides is identi ed in the pancreas tissue [54]. Some N-glycopeptides are elevated in chronic pancreatitis [55]. In MSAP, we found that uoroquinolone antibiotic was the most signi cantly decreased CARD functional composition. Fluoroquinolone has adequate tissue penetration and bactericidal properties in infected pancreatic necrosis [56]. Clinically, uoroquinolone has been recommended for antimicrobial therapy in AP [57]. Our ndings indicated that microbial resistance to aminocoumarin, glycopeptides and uoroquinolone antibiotics in MAP and MSAP may provide the research basis for the study of drug action in AP.
In the differential analysis of microbiota functional composition of VFDB, we found that FeoAB was the most signi cantly decreased VFDB functional composition in MAP. The feoAB operon of Escherichia coli is the most-characterized ferrous iron transport system, which promotes intracellular infection [58]. In MSAP, we found that capsular polysaccharide was the most signi cantly increased VFDB functional composition.
Several pathogenic bacteria, such as Escherichia-coli could produce capsular polysaccharide known as glycosaminoglycans for molecular camou age [59]. Maekawa T et al. found that patients with chronic pancreatitis were characterized by increased serum level of antibody against the capsular polysaccharide of Enterococcus faecalis [60]. In addition, we found that two-component system was the most signi cantly decreased VFDB functional composition in both MSAP and SAP. It is reported that peptidoglycan recognition proteins bind to the outer membrane of Escherichia coli and activates functionally homologous CpxA-CpxR two-component system, which modulates microbiome and in ammation [61]. Our study suggested that above virulent factors, including FeoAB, capsular polysaccharide, and the twocomponent system, may promote microbial infection and cause AP.

Conclusions
We conducted a shotgun metagenomics survey on the gut microbiome of AP patients with three different severity grades for the rst time. We identi ed several new AP-related bacteria and functional gene pathways, which extend the current knowledge about the role of gut microbiota in AP. Our ndings may be useful for future studies investigating the worsening mechanism of AP and developing strategies for the prevention, diagnosis, and treatment of AP. Needless to say that there are certain limitations in our study. Firstly, the sample size is small, and further longitudinal studies with a larger number of subjects are needed. Secondly, causal links between these microbiota and functional pathways and AP susceptibility need a deeper investigation.

Consent for publication
All the authors have consented for publication of this paper in the current form.
Availability of data and material All data generated or analyzed during this study are included in this published article.

Competing interests
The authors declare that they have no competing interests.  Difference analysis of COG functional composition in MSAP (A) and SAP (B) group. On the left, X and Yaxis represents the average relative abundance of the COG functional composition in various groups and the COG functional composition names at a certain classi cation level, respectively. On the right, X and Y axis represents different COG functional compositions between various groups and p-value of signi cance. *0.01<p<0.05, **0.01<p<0.001. U: intracellular tra cking, secretion, and vesicular transport; Z: cytoskeleton; C: energy production and conversion; T: signal transduction mechanisms; O: posttranslational modi cation, protein turnover, chaperones; W: extracellular structures.

Figure 2
Difference analysis of KEGG functional composition in MAP (A), MSAP (B) and SAP (C) group. On the left, X and Y-axis represents the average relative abundance of the KEGG functional composition in various groups and KEGG functional composition names at a certain classi cation level, respectively. On the right, X and Y axis represents different KEGG functional compositions between various groups and p-value of signi cance. *0.01<p<0.05, **0.01<p<0.001, ***p<0.001.  Difference analysis of ARDB functional composition in MAP (A), MSAP (B) and SAP (C) group. On the left, X and Y-axis represents the average relative abundance of the ARDB functional composition in various groups and the ARDB functional composition names at a certain classi cation level, respectively. On the right, X and Y axis represents different ARDB functional compositions between various groups and p-value of signi cance. *0.01<p<0.05, **0.01<p<0.001, ***p<0.001.

Figure 5
Difference analysis of CARD functional composition in MAP (A), MSAP (B) and SAP (C) group. On the left, X and Y-axis represents the average relative abundance of the CARD functional composition in various groups and the CARD functional composition names at a certain classi cation level, respectively. On the right, X and Y axis represents different CARD functional compositions between various groups and p-value of signi cance. *0.01<p<0.05, **0.01<p<0.001. Detailed annotation of each CARD functional composition name is available in the website of https://card.mcmaster.ca/ontology/39421.

Figure 6
Difference analysis of VFDB functional composition in MAP (A), MSAP (B) and SAP (C) group. On the left, X and Y-axis represents the average relative abundance of the VFDB functional composition in various groups and the VFDB functional composition names at a certain classi cation level, respectively. On the right, X and Y axis represents different VFDB functional compositions between various groups and p-value of signi cance. *0.01<p<0.05, **0.01<p<0.001, ***p<0.001.