PIK3AP1 and SPON2 Genes Are Differentially Methylated in Patients With Periodic Fever, Aphthous Stomatitis, Pharyngitis, and Adenitis (PFAPA) Syndrome

Periodic fever, aphthous stomatitis, pharyngitis, and adenitis (PFAPA) syndrome is the most common autoinflammatory disease in children and is often grouped together with hereditary periodic fever syndromes, although its cause and hereditary nature remain unexplained. We investigated whether differential DNA methylation was present in DNA from peripheral blood mononuclear cells (PBMC) in patients with PFAPA vs. healthy controls. A whole-epigenome analysis (MeDIP and MBD) was performed using pooled DNA libraries enriched for methylated genomic regions and identified candidate genes, two of which were further evaluated with methylation-specific restriction enzymes coupled with qPCR (MSRE-qPCR). The analysis showed that the PIK3AP1 and SPON2 gene regions are differentially methylated in patients with PFAPA. MSRE-qPCR proved to be a quick, reliable, and cost-effective method of confirming results from MeDIP and MBD. Our findings indicate that a B-cell adapter protein (PIK3AP1), as the PI3K binding inhibitor of inflammation, and spondin-2 (SPON2), as a pattern recognition molecule and integrin ligand, could play a role in the etiology of PFAPA. Their role and the impact of changed DNA methylation in PFAPA etiology and autoinflammation need further investigation.


INTRODUCTION
Periodic fever syndrome with adenitis, pharyngitis, and aphthous stomatitis (PFAPA) belongs to the group of autoinflammatory diseases (AID). This group of disorders is marked by increased inflammation associated with the innate immune system, and most of the disorders are inherited in a Mendelian pattern (1,2). Of the various autoinflammatory diseases, many now have a confirmed known genetic cause. However, PFAPA syndrome still has an unknown genetic background and pathogenesis (3).
PFAPA was first described by Marshall et al. (4). Its most common feature is periodic fever, while its other features are more variable: pharyngitis, aphthous stomatitis, and cervical adenopathies (5,6). Episodes have an early onset, usually in patients around 2-3 years of age, mainly before the age of five. Episodes can reoccur for several years, followed by disease remission, and leave no long-term health consequences (5)(6)(7). The syndrome is considered sporadic and has been described as a non-inherited syndrome; however, familial cases of PFAPA have been reported, suggesting a potential genetic origin (3,(8)(9)(10)(11). PFAPA could be caused by cytokine dysregulation linked to genetic variants of autoinflammatory disease-associated genes (12). Studies also point to altered complement activation and IL-1 production in PFAPA patients (13), as well as to IL-1β dysregulation (14). Since PFAPA shares clinical similarities with monogenic fever syndromes, several studies have investigated the possible involvement of the genes responsible for Familial Mediterranean Fever (MEFV gene), TNF-Receptor Associated Periodic Syndrome (TNFRSF1A gene), Mevalonate kinase deficiency (MVK gene), and Cryopyrin-Associated Periodic Syndrome (NLRP3 gene) in PFAPA cohorts (15). A variant in CARD8 has been associated with more severe cases of PFAPA. CARD8 acts as a negative regulator of the inflammasome; its decreased ability to bind NLRP3 could partly contribute to exacerbated inflammasome responses, although probably not on its own (16). Genes known to be involved in inflammation or in autoinflammatory disorders seem to contribute to a predisposition to PFAPA syndrome, suggesting complex genetic inheritance and interaction with non-genetic factors (12,14,(17)(18)(19)(20).
Several inflammasome-related genes have been found to have increased expression when demethylated during the differentiation of monocytes to macrophages, specifically AIM2, NLRC5, PYCARD, CASP1, and PSTPIP2, and their targets IL1A, IL1B, and IL1RN. Furthermore, untreated patients with Cryopyrin-Associated Periodic Syndrome (CAPS) have exacerbated DNA methylation-dependent regulation of the inflammasome product genes IL1B, ILR1N, NLRC5, and PYCARD, while patients with CAPS undergoing anti-IL-1β treatment have displayed demethylation levels in stimulated monocytes similar to those seen in healthy subjects (21). Different methylation levels of the MEFV gene have been observed in Familial Mediterranean Fever patients compared to healthy controls (22). Additionally, DNA methylation plays an important role during hematopoietic differentiation to a myeloid vs. a lymphoid lineage (23).
Since studies suggest that DNA methylation plays a role in other autoinflammatory diseases, we hypothesized that specific methylation patterns may be aberrant in at least a small portion (population) of peripheral blood mononuclear cells (PBMC) in PFAPA patients. To identify the potentially relatively small change in the methylation patterns due to the specifics of a PBMC-derived DNA sample and to confirm whether the PFAPA cohort differs in DNA methylation patterns from healthy controls, a whole-epigenome analysis was performed using pooled DNA libraries enriched for methylated genomic regions using Methylated DNA Immunoprecipitation (MeDIP) and Methyl-CpG-binding domain (MBD). We identified several candidate genes with differential methylation, two of which (PIK3AP1 and SPON2) were chosen based on their involvement in the inflammation pathways and were further evaluated with MSRE-qPCR. To our knowledge, no research had previously been performed regarding differentially methylated DNA in PFAPA patients.

Participants
Clinical data and samples of 75 patients (44 boys and 31 girls) with PFAPA syndrome at the University Children's Hospital Ljubljana were collected from 2008 to 2016. The median age of the patients when blood was taken for the DNA analysis was 4.1 [Interquartile range, IQR 3-5.8] years. Patients were not in an active state of the disease when samples for DNA isolation were taken. Prior to donating blood for DNA analysis, 15 (20%) had received methylprednisolone in the past but not in the month prior to blood donation (except for one individual, who had received treatment 2 weeks prior). Samples from 65 apparently healthy children (35 boys and 30 girls) of Slovenian ethnicity whose median age was 5.3  were included in the study as healthy controls. The parents of each child included in the study were informed about the aim of the study and signed a written informed consent form for inclusion in the study. The study was approved by the Ethics Committee of the Republic of Slovenia and was conducted according to the principles of the Helsinki Declaration.

Samples
Five ml of peripheral blood was taken for DNA isolation. Blood was taken from patients during routine venipuncture at follow-up visits. Blood of the healthy controls was taken during routine health examinations of the children. DNA isolation was performed using the FlexiGene isolation kit (Qiagen, Germany) according to the recommended protocol. The DNA was stored at 4 • C.

MeDIP and MBD Data Analysis
For each sample pool (PFAPA vs. healthy), three separate NGS libraries were generated for each enrichment method (MBD or MeDIP). All 36 libraries were combined and simultaneously sequenced on the MiSeq Illumina sequencer [MiSeq Reagent Kit v3 (150-cycle)] by a 2 × 75 paired-end run. Sequencing was repeated until each separate library had acquired at least 10 million pair-end reads. The acquired dataset was filtered for low-quality reads and duplicates aligned by the BWA-MEM aligner (24) and followed by MACS2 (25) analysis of narrow and broad peak regions. The callable regions identified by the MACS2 algorithm were used to count reads in each peak region for each library by the featureCounts algorithm (26). Using DeSeq2 (27), the algorithm peaks were compared to identify differentially methylated regions via differences in normalized read counts. The Benjamini-Hochberg procedure, which controls false discovery rate (FDR), was used to identify true differentially methylated regions (DMRs). DMRs were annotated for the overlapping or proximal gene regions. DMRs associated with the immune system were further verified using MSRE-qPCR. After successful evaluation for specificity and annealing temperature, evaluation of the qPCR assay was performed. We used Luna Universal qPCR Master Mix (New England BioLabs) with a fast cycling profile. Standard curves were prepared for each primer pair to determine the efficiency (E) of the designed primers. Efficiency calculations were calculated online with the NEBioCalculator for qPCR Quantification (https:// nebiocalculator.neb.com/#!/qPCRGen) ( Table 2).

Enzyme Restriction
The regions selected for MSRE-qPCR had multiple CpGs between forward and reverse primer annealing sites. The selected enzymes (MspJI and McrBC; Tables 3, 4) cover almost all CpGs; however, for successful digestion and subsequent evaluation, at least one of the targeted CpGs must be methylated. Furthermore, the MspJI enzyme can also recognize methylated C in front of any base (Table 3). However, since, in humans, mostly CpGs are  methylated, only CpGs were counted for the number of candidate restriction sites ( Table 4). Beta globin amplicon, which does not contain any CpGs, was used as a reference. Equal amounts of DNA were used in control and restriction reaction. First, 200 ng of DNA was digested in accordance with the manufacturer's protocol by 2 units of MspJI and 2 units of McrBC (both New England BioLabs) in a 25-µl reaction with CutSmart buffer for 1 h at 37 • C, and this was followed by heat inactivation for 20 min at 65 • C. Samples from the restriction reaction were then purified using Sample Purification Beads from a TruSight R One Sequencing Panel Kit according to the manufacturer's protocol.

Quantitative PCR
The purified digested and undigested DNA samples were used for relative qPCR quantification. qPCR assays were run in triplicates in 96-well plates using Luna Universal qPCR Master Mix (SYBR). In addition, melt curves were performed. Five pairs of primers [two for PIK3AP1, two for SPON2, and one for normalization (HBB, beta globin)] were analyzed on digested and undigested DNA samples for relative quantification using the C T method (28). This method analyzes relative changes in methylation of the target gene to the reference gene, and it assumes uniform PCR amplification efficiency across all reactions. Due to varying efficiencies in qPCR assays, we used efficiency corrected calculation (Equations 1-3) (28) for each obtained CT to produce more accurate estimates in relative quantification.
Calculations for the C T method were done in Excel. Statistical data analysis was performed in the GraphPad Prism 8 program. A Mann-Whitney test was performed; the data represent relative expression ratios obtained with Equation (1). P < 0.05 was considered statistically significant, and p < 0.01 was considered statistically highly significant.

Clinical Characteristics of PFAPA Patients
All included participants fulfilled the clinical criteria for PFAPA syndrome (5). Included in our study were 44 (58.7%) boys and 31 (41.3%) girls. All patients were asymptomatic during the afebrile period. Pharyngitis (94.7%) and adenitis (89.3%) were the most present symptoms during febrile episodes, followed by abdominal pain (61.3%) and aphthous stomatitis (58.7%) ( Table 5). Included in our cohort were three pairs of siblings (twin brothers, twin sisters, and two brothers). Positive family history, meaning that at least one first-degree relative had recurrent fevers or a tonsillectomy, was found in 73.4% of patients. Positive family history, either in the first or second degree, was present in 87.5% of patients. To reduce the number of false-positive results in DMR analysis using DeSeq2, the FDR cutoff value was modified in such a way that the potential DMR set contained up to 1 (one) false-positive result. Consequently, the FDR cutoff value was set to q < 0.06 for the MeDIP dataset and q < 0.15 for the MBD dataset. After count tables

Identification of DMRs With MSRE-qPCR
Of the differentially methylated sequences previously identified by using MDB and MeDIP, two regions of interest (ROIs) were chosen to perform our analysis based on the location and function of the gene. Of the top significant candidate genes, one was chosen from MeDIP (PIK3AP1, region 10:98425909-98426208, GRCh37) and one from the MBD method (SPON2, region 4:1163349-1163641, GRCh37). Relative quantification of digested and undigested DNA was performed with MSRE-qPCR in order to compare patients with healthy controls. The second intron region of PIK3AP1 (10:98425909-98426208) was found to  Table 6).
Methylation results in the fifth intron region of the SPON2 gene are conflicting. MBD defined the region of the fifth intron of SPON2 (4:1163349-1163641) as more methylated in PFAPA patients, while MSRE-qPCR later showed that this is not the case for the whole region. The first half (4:1163329-1163502) of the fifth intron was, in fact, less (P = 0.001) methylated and the second half (4:1163480-1163727) was more (P = 0.0191) methylated in PFAPA patients compared to the healthy controls.
Results from both MeDIP and MBD were confirmed with MSRE-qPCR, with some differences. Firstly, MSRE-qPCR identifies smaller regions where differential methylation occurs, because qPCR reaction efficiency is limited with amplicon length. Secondly, because MSRE-qPCR has an amplicon size limit, it was able to reveal that the whole region of the fifth intron of the SPON2 gene does not have higher methylation in PFAPA patients. Moreover, MSRE-qPCR requires only one methylated CpG at the ROI per single DNA molecule to detect a difference, as one cut site is enough to prevent PCR amplification of a particular DNA molecule.
An unpaired Mann-Whitney test with a 95% confidence level was performed. Patients vs. controls were compared for each ROI. Exact two-tailed P-values are listed in Table 6. The data are visualized with boxplots with added scatter plots and the Pvalue of the Mann-Whitney test in Figure 1. Each dot represents the relative expression ratio of each control or patient obtained with Equation (1). The number of controls and patients used for MSRE-qPCR analysis was lower than the initial 65 controls and 76 patients due to the limited amount of DNA for specific samples and difficulties with qPCR.

DISCUSSION
PFAPA syndrome is the most common pediatric fever syndrome, and its etiology is unknown. So far, no clear genetic causes have been found. In this study, we analyzed differential methylation patterns in a cohort of 75 PFAPA patients. All patients included in the study fulfilled the clinical criteria for PFAPA syndrome (5). The majority of patients analyzed were boys (58.7%). During the fever episode, pharyngitis and adenitis were the most common symptoms (94.7 and 89.3%, respectively). All three major symptoms (adenitis, pharyngitis, and aphthous stomatitis) were present in 52% of the patients. The included patients had previously been analyzed in an independent study, and genetic analysis was performed for four genes: AIM2, NLRP3, MEFV, and MVK. No clinically significant variants were found (18).
Relative quantification with the MSRE-qPCR method proved to be a quick and efficient way to estimate differentially methylated DNA regions, which were identified beforehand by MeDIP and MBD. Our results showed that there are two regions, one inside the second intron of the PIK3AP1 gene and a second one inside the fifth intron of SPON2, which have a statistically significant difference in methylation in PFAPA patients when compared to healthy controls.
Both MeDIP and MBD have the advantage of identifying regions with differential methylation over the whole genome, while MSRE-qPCR is more suitable to analyze the few previously identified candidate regions. MeDIP favors regions with low CpG density and MBD favors regions with higher CpG density (29), which was also observed in our case. All three methods are a suitable way to identify differentially methylated genomic regions, though they cannot determine how many and which CpGs are methylated specifically.
The PIK3AP1 gene was identified as more methylated in PFAPA patients by MeDIP, which the MSRE-qPCR method confirmed, but only for a part of the said region. The MBD results identified a region of SPON2 as more methylated in PFAPA patients, while MSRE-qPCR later revealed that this is not the case for the whole identified region. The first half had lower methylation levels in PFAPA patients, while the second half had higher methylation levels. A possible reason for such divergence with the results of the first half of SPON2 could be the specifics of the process of broad peak identification by the MACS2 algorithm, where a wider area is incorporated into the final signal peak identification, partially influenced by the final NGS library insert size as well. Nevertheless, there is a measurable difference in the methylation pattern in the fifth intron region of the SPON2 gene.
The B-cell adapter protein (BCAP, PIK3AP1 gene) is a phosphoinositide 3-kinase (PI3K) binding protein (30) and FIGURE 1 | Comparison of MSRE-qPCR ratios between Control and Patient group for each ROI. Boxplots include scattered dots that represent relative expression ratio of Control (blue) or Patient (pink) sample. Also included at the top are exact two-tailed P-values of unpaired Mann-Whitney tests. All scatter plots are plotted with all data, except the plot for PIK_2 ROI, which was plotted with data without outliers (ROUT method, Q = 0.1%) for visualization purposes. Mann-Whitney tests for all ROIs were performed with all data (outliers included).
an important inhibitor of proliferation and myeloid cell differentiation that works in a cell-intrinsic manner (31). The PI3K signaling cascade influences cell proliferation and survival, metabolic reprogramming, and cellular migration. As a PI3K adaptor, BCAP is a key regulator of PI3K signaling and T-cell development into effector and memory cells (32). It is expressed in lymphoid as well as in myeloid cell populations (31)(32)(33)(34). BCAP regulates inflammatory response. BCAP deficiency results in exaggerated innate immune response, leading to higher CD4 + activation (35) and to more proliferative cells (31). As a macrophage signaling adaptor protein, it can dampen NLRP3 and NLRC4 inflammasome activation through interaction with the caspase-1 inhibitor (36). Altogether, BCAP could play a part in influencing systemic inflammation, which PFAPA syndrome is known for.
Spondin-2 (also called mindin, SPON2 gene) is an extracellular matrix protein that functions as a pattern recognition molecule (PRM) for initiating innate immune responses, as well as an integrin ligand for inflammatory cell recruitment and T-cell priming (37)(38)(39). It has been shown that, in vivo, it is crucial for the efficient clearance of bacterial (37) and viral (40) infections. Spondin-2/mindin-induced signaling could be as important as other, better-defined signaling pathways, such as TLR signaling. He et al. have proposed that mindinmediated carbohydrate recognition of microbial pathogens represents a secondary stimulation essential for activation of innate immune cells (37). Additionally, it has been suggested that it has a role in the immune response against tumor cell growth and migration (41). Little else is known of spondin-2 regarding immune regulation; however, as an extracellular PRM, it could have a role in inflammation, since extracellular PRMs are able to complement activation, opsonization, agglutination, neutralization, and regulation of inflammation (42).
Epigenetic modifications, including DNA methylation, can affect gene expression and, consequently, can be used as a disease biomarker (43). Methylation of CpG dinucleotides is one of the principal epigenetic mechanisms (44). Regions rich in CpGs are called CpG islands (CGI). Unmethylated promotor CGIs are generally associated with transcriptionally active genes, whereas hypermethylated promotor CGIs result in gene transcription repression. Furthermore, methylation levels inside the gene correlate with gene expression, especially the first intron methylation (45). The majority of research on disease-related DNA methylation has been done in the field of cancer (43); however, there is increasing evidence that epigenetic dysregulation plays a role in autoinflammatory diseases as well (46).
The role of DNA methylation of intragenic regions is less clear, though evidence exists that methylated intragenic regions influence gene transcription (47). Partial methylation of the coding region can inhibit gene expression (48,49), and even a few methylated cytosines can inhibit a flanking promoter, but a threshold of modified sites is required (50). Elevated DNA methylation in intragenic regions usually correlates with silencing of the associated gene (51). Blattler et al. showed that changes in DNA methylation within bodies of genes played a much larger role than changes in promotor regions (52). The methylation level of the first intron is inversely associated with gene expression, and this association is conserved regardless of the species or tissue. Methylation tends to increase with distance from the first exon-first intron boundary, and its effect on gene transcription decreases with downstream distance from the first exon (45,53). The study of DNA methylation during monocyte differentiation revealed that a high proportion of the changes occurred in non-promoter regulatory regions, mainly enhancers in gene bodies and intergenic regions (54). Methylation patterns in the MEFV exon in pediatric patients with Familial Mediterranean Fever correlated with expression of the same gene; the observed slight increase in DNA methylation of the second exon in patients correlated with decreased expression (22).
The observed differential methylation of the second intron in PIK3AP1 (10:98425953-98426190) could be of great significance regarding changed expression. BCAP functions as a checkpoint to restrict TLR signaling and production of inflammatory cytokines (55). We hypothesize that the higher methylation found in the second intron of PIK3AP1 in PFAPA patients could lead to lower expression of the BCAP protein and cause a disrupted inhibition of inflammation, leading to exaggerated inflammation or response to environmental stimuli (for example an infection). The cause of this particular differential methylation is unknown, and DNA methylation can be influenced by a range of external factors, such as diet, drugs, and infections (56,57). There is evidence of infection-induced hypermethylation of PIK3AP1 promoter and downregulation of its expression (58) and evidence of hypomethylation and increased expression of PIK3AP1 triggered by low levels of folic acid (59). Infections are proposed contributors to the pathogenesis of PFAPA (15,60).
The second half of the spondin-2 region (4:1163480-1163727, GRCh37) that was investigated and showed higher methylation in PFAPA patients is also part of a CTCF binding region (4: 1163601-1163800, GRCh37). CTCF is a transcription factor, chromatin organizer, and insulator protein that was initially discovered as a transcriptional repressor; however, its exact mechanisms remain unknown (61). We hypothesize that the observed increased methylation could potentially inhibit CTCF binding. As an insulator protein, CTCF can influence gene transcription by restricting the binding of transcription enhancers. Its binding sites are located far from the transcriptional start sites, and their distribution is strongly correlated with genes (62). However, without gene expression analysis, we cannot confirm whether the changed methylation in PIK3AP1 and SPON2 influences their expression and in what way.
DNA was isolated from whole blood, meaning we analyzed methylation in DNA from the whole blood cell population with unknown ratios. Samples were collected during an afebrile period, patients were asymptomatic, and we assumed that blood cell populations were normal and comparable based on the facts that changes in white blood cell counts are associated with febrile episodes (13,63) and that the concentrations of white blood cells in afebrile PFAPA patients are comparable to those in healthy children (63). Our results indicate that there is a portion of blood cells with differentially methylated DNA in these two genes and that the number of cells or the specific clonal subpopulation of cells was high enough to generate a measurable signal. Currently, there is no clear evidence which cell subpopulation is involved. High monocytes and high neutrophils seem to be the most consistently reported changes in blood cell populations during PFAPA febrile attacks (13,63,64) and could be the possible source of the differential DNA methylation we observed.
A limitation of our study, in addition to the mixed cell populations as our source of DNA, is the fact that 15 patients (20%) had previously received methylprednisolone, which could affect their DNA methylation patterns. We do not know if it affected the regions in question or in what way and for how long. Steroids do have an effect on the disease, but since their effect is not permanent, they only stop the episode and not the disease itself (65)(66)(67). In order to overcome this limitation, we reanalyzed data excluding patients that received steroid treatment and obtained similar results. Significance did not change, except for the second half of the SPON2 region, which crossed the significance line but still showed a tendency of hypermethylation. Loss of significance could also be partly attributed to reduced sample size. Additional studies investigating the effect of methylprednisolone should be conducted. Another limitation would be that the test and control groups were not of the same age (4.1 and 5.3 years, respectively). The median age of our control group was intentionally above five, which is when the fevers generally begin to appear (5-7). Altogether, the results from our methods do not inform us of the identity and number of CpGs involved. We also cannot tell in which cells this difference in methylation occurs, since the source of our analyzed DNA was whole blood. For a more detailed look into the possible role of DNA methylation in the pathogenesis of PFAPA, DNA isolated from different cell lineages should be analyzed.
Both the B-cell adapter protein through the PI3K activation pathway and spondin-2 as an extracellular pattern recognition molecule are involved in the primary stages of immune responses. BCAP acts as an inhibitor of receptor signaling, mainly TLR and spondin-2, possibly as a recruiter or activator of T cells, or by working as an opsonin that activates complement. The role of changed DNA methylation in autoinflammation certainly needs further investigation, as does the role of BCAP and spondin-2 in the etiology of PFAPA. Differential methylation of genes PIK3AP1 and SPON2 is not reported in association with other similar autoinflammatory diseases; however, PIK3AP1 was associated with autoimmunity (68), and its increased expression was shown to promote TLR7-driven lupus-like disease (69). Differential expression and methylation of both genes have been linked with cancer. Changed methylation of PIK3AP1 is associated with neuroblastoma (70), and its upregulation is associated with Waldenström macroglobulinemia (71). Hypomethylation of the promoter of SPON2 and its increased expression is associated with prostate cancer (72,73); its upregulation is also associated with colorectal cancer (74).
As far as the clinical impact of these data is concerned, this kind of change in methylation would not be detected by Sanger sequencing, because the proportion of differentially methylated DNA molecules is too small. Therefore, a more advanced quantification method is needed, such as MSRE-qPCR. Nevertheless, the cutoff values for the positive MSRE-qPCR results should be carefully examined, set, and validated.

CONCLUSIONS
Whole-genome methylation screening analysis methods, such as MeDIP and MBD are a useful tool for identifying differentially methylated genomic regions. MSRE-qPCR proved to be a reliable, quick, and cost-effective method of confirming results and showed potential applicability in translation of the research into clinical practice. The changed methylation patterns in PIK3AP1 and SPON2 that we observed in PFAPA patients point to novel and still unknown roles of BCAP and spondin-2 in the etiology of PFAPA. Furthermore, likely transient changes of DNA methylation patterns are potentially a novel direction in research into the molecular mechanisms leading to PFAPA development.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Ethics Committee of the Republic of Slovenia. Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin.

AUTHOR CONTRIBUTIONS
EL, JK, TT, and TR contributed the conception and design of the study and performed the actual experiments and analysis. EL wrote the first draft of the manuscript. JK wrote the sections of the manuscript. NT and DP selected the included patients and collected data and samples. MD and TA contributed equally and were responsible for the final approval of the submitted version. All authors contributed to the article and approved the submitted version.