WGS of Commensal Neisseria Reveals Acquisition of a New Ribosomal Protection Protein (MsrD) as a Possible Explanation for High Level Azithromycin Resistance in Belgium

In this study, we characterized all oropharyngeal and anorectal isolates of Neisseria spp. in a cohort of men who have sex with men. This resulted in a panel of pathogenic Neisseria (N. gonorrhoeae [n = 5] and N. meningitidis [n = 5]) and nonpathogenic Neisseria (N. subflava [n = 11], N. mucosa [n = 3] and N. oralis [n = 2]). A high proportion of strains in this panel were resistant to azithromycin (18/26) and ceftriaxone (3/26). Whole genome sequencing (WGS) of these strains identified numerous mutations that are known to confer reduced susceptibility to azithromycin and ceftriaxone in N. gonorrhoeae. The presence or absence of these known mutations did not explain the high level resistance to azithromycin (>256 mg/L) in the nonpathogenic isolates (8/16). After screening for antimicrobial resistance (AMR) genes, we found a ribosomal protection protein, Msr(D), in these highly azithromycin resistant nonpathogenic strains. The complete integration site originated from Streptococcus pneumoniae and is associated with high level resistance to azithromycin in many other bacterial species. This novel AMR resistance mechanism to azithromycin in nonpathogenic Neisseria could be a public health concern if it were to be transmitted to pathogenic Neisseria. This study demonstrates the utility of WGS-based surveillance of nonpathogenic Neisseria.


Introduction
Neisseria gonorrhoeae has developed antimicrobial resistance (AMR) to every class of antimicrobials used to treat it [1]. Of particular concern are the rapid increases in azithromycin resistance and the emergence of combined ceftriaxone and high-level azithromycin resistance [1][2][3][4]. Much of these resistance mechanisms are acquired via horizontal gene transfer (HGT) from commensal Neisseria [5][6][7][8]. The importance of this pathway for the emergence of AMR has been well established for extended spectrum cephalosporins, such as ceftriaxone. Phylogenetic analyses confirmed that the HGT of a section of the penA gene from commensal Neisseria played a crucial role in the genesis of ceftriaxone resistance in N. gonorrhoeae [7,9]. As a result, one of the crucial first steps in the emergence of ceftriaxone resistance in N. gonorrhoeae was the selection of ceftriaxone resistance in commensal Neisseria [7,[9][10][11]. Similarly, HGT played an important role in the genesis of AMR to macrolides such as azithromycin [8].
An important reason why AMR would be expected to emerge in commensal Neisseria prior to N. gonorrhoeae is the considerably higher prevalence of commensals in human populations. Whilst the prevalence of N. gonorrhoeae is typically a fraction of a percent in general populations and only reaches 10% in core groups [1,12], the prevalence of commensal Neisseria is close to 100% in the oropharynx, where they form an important part of a healthy microbiome [6,13]. This makes them more likely to develop AMR in response to high antimicrobial consumption compared to N. gonorrhoeae [10,14]. Because this AMR can be transferred to N. gonorrhoeae and monitoring AMR in commensals is not complicated, it has been proposed that the surveillance of AMR in commensals could be used as an early warning system to detect the risk of AMR in N. gonorrhoeae and other bacteria [10,14,15]. AMR in N. gonorrhoeae, Mycoplasma genitalium and Treponema pallidum has frequently emerged in core-groups with high rates of partner turnover and high antimicrobial consumption [14,[16][17][18]. This provides the rationale for monitoring the antimicrobial susceptibility of commensal Neisseria in these core-groups. In our setting, men who have sex with men (MSM) attending our STI clinic represent one such core-group [17,19]. In a pilot study, we assessed ceftriaxone and azithromycin susceptibilities of all Neisseria spp. isolated (n = 26) pre-and post-treatment from 10 MSM attending our STI clinic with a diagnosis of anorectal gonorrhea. The minimum inhibitory concentrations (MICs) were found to be alarmingly high [15]. The most prevalent commensal, N. subflava, for example, was found to have a median azithromycin MIC of 176 mg/L (IQR 0.047-256) and a median ceftriaxone MIC of 0.38 mg/L (IQR 0.023-2) [15]. These values were considerably higher compared to Belgian historical samples and moreover, were the highest MIC values for N. subflava ever published [5,14,15,20].
In this paper, we reported the antimicrobial resistance-associated mutations detected in these Neisseria species with a focus on those implicated in macrolide and extendedspectrum cephalosporin resistance.

Characterization of Strain Collection
Twenty-six clinical Neisseria species were isolated and subjected to whole genome sequencing. De novo assemblies resulted in genome sizes ranging from 2.09 Mbp to 2.82 Mbp with the number of putative coding DNA sequence (CDS) ranging from 1983 to 2535 (Table S1).
Prior species identification by matrix-assisted laser desorption/ionization-time of flight mass spectrometry (MALDI-TOF-MS) was verified by analyzing the whole genome assemblies using BIGSdb and the 50S ribosomal gene, rplF [15,21,22]. There was complete agreement between these three methods except that MALDI-TOF misidentified N. mucosa as N. macacae, and rplF gene annotation could not discriminate between N. mucosa and N. oralis. In the present study, screening by BIGSdb was therefore used in species identification. This resulted in a Neisseria panel with a composition of pathogenic Neisseria: N. gonorrhoeae (n = 5) and N. meningitidis (n = 5) and nonpathogenic Neisseria: N. subflava (n = 11), N. mucosa (n = 3) and N. oralis (n = 2).
The antimicrobial susceptibility testing (AST) of this panel showed higher azithromycin MICs in commensal Neisseria isolates than in pathogenic Neisseria (commensal MIC median 140 mg/L (IQR 6-256), pathogens median 0.75 mg/L (IQR 0.31-1.00; p-value < 0.0001; Table 1; Figure S1a). The same was true for ceftriaxone (commensal MIC median 0.082 mg/L (IQR 0.043-0.125), pathogens median 0.016 mg/L (IQR 0.016-0.016; p -value < 0.0001; Table 1; Figure S1b). For illustrative purposes, we applied The European Committee on Antimicrobial Susceptibility Testing (EUCAST) definitions of gonococcal resistance to azithromycin and ceftriaxone to define AMR for all the Neisseria spp. [23] (Figure 1). Of note, all the commensal Neisseria were classified as resistant to azithromycin (MIC > 1 mg/L). Resistance was particularly evident in the N. subflava strains, of which 8/11 had azithromycin MICs greater than 256 mg/L. Three of these eight isolates were classified as resistant to ceftriaxone.   3 Minimum inhibitory concentrations of azithromycin (AZM) and ceftriaxone (CRO). 4    The combined MICs of all post-treatment isolates (visit 2) where higher than those of the pre-treatment isolates (visit 1): the azithromycin MIC increased from a median of 1.5 mg/L (IQR 0.016-0.047) to 256 mg/L (IQR 24-256; p-value 0.0012) and ceftriaxone MIC from 0.016 mg/L (IQR 0.016-0.047) to 0.094 mg/L (IQR 0.047-0.125; p-value 0.0059). This effect was mainly driven by the elimination of the more susceptible pathogenic Neisseria between visits 1 (n = 8) and 2 (n = 0; Figure 2). There was no statistically significant increase in azithromycin or ceftriaxone MIC obtained from visit 2 compared to visit 1 in the N. subflava isolates alone. The sample sizes were too small to evaluate MIC changes within other species.
Pathogens 2021, 10, x FOR PEER REVIEW 4 of 15 The combined MICs of all post-treatment isolates (visit 2) where higher than those of the pre-treatment isolates (visit 1): the azithromycin MIC increased from a median of 1.5 mg/L (IQR 0.016-0.047) to 256 mg/L (IQR 24-256; p-value 0.0012) and ceftriaxone MIC from 0.016 mg/L (IQR 0.016-0.047) to 0.094 mg/L (IQR 0.047-0.125; p-value 0.0059). This effect was mainly driven by the elimination of the more susceptible pathogenic Neisseria between visits 1 (n = 8) and 2 (n = 0; Figure 2). There was no statistically significant increase in azithromycin or ceftriaxone MIC obtained from visit 2 compared to visit 1 in the N. subflava isolates alone. The sample sizes were too small to evaluate MIC changes within other species.
No RAMs in PilQ or 23S rRNA were found in either the commensal or pathogenic Neisseria.
porB: In the commensals, up to four copies of full length porB were detected per isolate. In N. subflava, four paralogs were detected. These varied in size between 997 to 1198 bp. In N. mucosa, three paralogs were found (length of 984 bp, 1074 bp and 1122 bp), whereas in N. oralis, two paralogs of 1098 bp and 1125 bp were present. Mutations in the loop III region in porB, known as penB, are associated with reduced susceptibility to antimicrobials such as ceftriaxone [24,25]. Whilst the transmembrane domains of porB genes in commensals are relatively conserved, the loop III structure is hypervariable among Neisseria species and therefore, penB mutations are difficult to interpret in commensals ( Figure S2). mtrCDE: Alignment of the mtr promoter showed variable and conserved regions between the Neisseria species (Figure 4). More specifically, no differences were found in the -10 mtrCDE and -10 mtrR promoter sites, but there were variations in the mtrR binding site and in the inverted repeat. An A to C conversion in the inverted repeat of the -35 mtrCDE promoter was detected in N. gonorrhoeae (1/5). This A to C conversion was also present in N. meningiditis (5/5), N.oralis (2/2) and N. mucosa (3/3). All strains of N. subflava (10/10) contained a G at this position. One strain of N. meningitidis (MO00021/1) contained a single T insertion in the inverted repeat region, which has previously been described in a clinical N. gonorrhoeae strain [26]. One N. gonorrhoeae strain (MO00014/1, MIC of 2 mg/L for azithromycin) contained the A to C conversion in the inverted repeat region as well as a G to A conversion in the -35 mtrR promoter site. This G to A conversion was also detected in 2 N. subflava strains (MO00036/1 and MO00036/2). A commonly described disruption of the -35 mtrCDE promoter, deletion of one A in the inverted repeat was not found in this panel [27]. porB: In the commensals, up to four copies of full length porB were detected per isolate. In N. subflava, four paralogs were detected. These varied in size between 997 to 1198 bp. In N. mucosa, three paralogs were found (length of 984 bp, 1074 bp and 1122 bp), whereas in N. oralis, two paralogs of 1098 bp and 1125 bp were present. Mutations in the loop III region in porB, known as penB, are associated with reduced susceptibility to antimicrobials such as ceftriaxone [24,25]. Whilst the transmembrane domains of porB genes in commensals are relatively conserved, the loop III structure is hypervariable among Neisseria species and therefore, penB mutations are difficult to interpret in commensals ( Figure S2). mtrCDE: Alignment of the mtr promoter showed variable and conserved regions between the Neisseria species ( Figure 4). More specifically, no differences were found in the -10 mtrCDE and -10 mtrR promoter sites, but there were variations in the mtrR binding site and in the inverted repeat. An A to C conversion in the inverted repeat of the -35 mtrCDE promoter was detected in N. gonorrhoeae (1/5). This A to C conversion was also present in N. meningiditis (5/5), N.oralis (2/2) and N. mucosa (3/3). All strains of N. subflava (10/10) contained a G at this position. One strain of N. meningitidis (MO00021/1) contained a single T insertion in the inverted repeat region, which has previously been described in a clinical N. gonorrhoeae strain [26]. One N. gonorrhoeae strain (MO00014/1, MIC of 2 mg/L for azithromycin) contained the A to C conversion in the inverted repeat region as well as a G to A conversion in the -35 mtrR promoter site. This G to A conversion was also detected in 2 N. subflava strains (MO00036/1 and MO00036/2). A commonly described disruption of the -35 mtrCDE promoter, deletion of one A in the inverted repeat was not found in this panel [27].

AMR Gene Screening
In 13 commensals and in one pathogenic Neisseria isolate, gene acquisition was detected ( Figure 5). In total, six different genes related to resistance to five different classes of antimicrobials were found. Only one of these is implicated in resistance to macrolides -msr(D).

AMR Gene Screening
In 13 commensals and in one pathogenic Neisseria isolate, gene acquisition was detected ( Figure 5). In total, six different genes related to resistance to five different classes of antimicrobials were found. Only one of these is implicated in resistance to macrolides -msr(D).

Msr(D):
The msr(D) gene has been shown to confer macrolide resistance in numerous species by acting as a ribosomal protection protein and displacing the macrolides from the ribosome [28][29][30]. It has not been reported in Neisseria before. However, it was present in 9 out of 11 isolates of N.subflava in our panel. Its genomic location and organization were highly conserved in all nine N. subflava isolates ( Figure S3). A 1944 bp sequence was present at 30 bp downstream of a DNA uptake sequence (DUS) and continued 457 bp downstream and 23 bp upstream of msr(D) (Figure 6). This entire 1944 bp sequence was 100% identical to a portion (position 2006 to 3950) of the macrolide efflux genetic assembly (MEGA) element in Streptococcus pneumoniae [31][32][33]. This element encodes five open reading frames (ORFs); mef(A), msr(D), ORF6, ORF7 and umuC-mucB. The integrated portion in N. subflava includes the complete msr(D) gene and a 338 bp section from the C terminus of mef(A). The msr(D) gene is part of the antibiotic resistance ATP Binding Cassette (ARE ABC-F) protein family, which all have a similar basic structure. A unique domain in this family is the P site tRNA Interaction Motif (PtlM) which is thought to play a crucial role in displacing bound macrolides from the 50S ribosome [34]. This motif as well as two other ABC-F unique motifs were found in the msr(D) gene ( Figure S4). To assess whether msrD may be playing a role in azithromycin resistance, genotypic azithromycin MIC prediction models with and without msrD were compared. These models revealed that the known azithromycin RAMs were only able to explain 21% of the variation in azithromycin MICs in nonpathogenic Neisseria (adjusted R-squared: 0.21, p-value: 0.18), compared to 96% in pathogenic Neisseria (adjusted R-squared 0.96, p-value: 0.001). The inclusion of a binary variable indicating the presence or absence of the msr(D) gene improved the predictive ability of the model from 21 to 84% in nonpathogenic Neisseria (adjusted R-squared: 0.84, p-value: <0.001).

Msr(D):
The msr(D) gene has been shown to confer macrolide resistance in numerous species by acting as a ribosomal protection protein and displacing the macrolides from the ribosome [28][29][30]. It has not been reported in Neisseria before. However, it was present in 9 out of 11 isolates of N.subflava in our panel. Its genomic location and organization were highly conserved in all nine N. subflava isolates ( Figure S3). A 1944 bp sequence was present at 30 bp downstream of a DNA uptake sequence (DUS) and continued 457 bp downstream and 23 bp upstream of msr(D) (Figure 6). This entire 1944 bp sequence was 100% identical to a portion (position 2006 to 3950) of the macrolide efflux genetic assembly (MEGA) element in Streptococcus pneumoniae [31][32][33]. This element encodes five open reading frames (ORFs); mef(A), msr(D), ORF6, ORF7 and umuC-mucB. The integrated portion in N. subflava includes the complete msr(D) gene and a 338 bp section from the C terminus of mef(A). The msr(D) gene is part of the antibiotic resistance ATP Binding Cassette (ARE ABC-F) protein family, which all have a similar basic structure. A unique domain in this family is the P site tRNA Interaction Motif (PtlM) which is thought to play a crucial role in displacing bound macrolides from the 50S ribosome [34]. This motif as well as two other ABC-F unique motifs were found in the msr(D) gene ( Figure S4). To assess whether msrD may be playing a role in azithromycin resistance, genotypic azithromycin MIC prediction models with and without msrD were compared. These models revealed that the known azithromycin RAMs were only able to explain 21% of the variation in azithromycin MICs in nonpathogenic Neisseria (adjusted R-squared: 0.21, p-value: 0.18), compared to 96% in pathogenic Neisseria (adjusted R-squared 0.96, p-value: 0.001). The inclusion of a binary variable indicating the presence or absence of the msr(D) gene improved the predictive  tetM: The tetM gene was detected in one N. gonorrhoeae and four N.subflava strains ( Figure 5). Whilst the tetM gene was carried on a plasmid in N. gonorrhoeae, in N. subflava it was integrated into the chromosome ( Figure S5).

tetM:
The tetM gene was detected in one N. gonorrhoeae and four N.subflava strains ( Figure 5). Whilst the tetM gene was carried on a plasmid in N. gonorrhoeae, in N. subflava it was integrated into the chromosome ( Figure S5).

Core Genome (cg) MLST
A gene-by-gene approach was used to assess the cgMLST/wgMLST allelic loci differences. In total, 3799 allelic loci were identified using the study strains (n = 26) and the reference strains (n = 7). For constructing the cgMLST tree, 1038 allelic loci were used (Figure 7). Distinct clusters were observed and the strains clustered closely according to their respective species. Maximum allelic differences between N. gonorrhoeae reference strains and the N. subflava core genome were >3000. The maximum allelic difference of the two strains which lacked the msr(D) had allelic loci differences of >4000. Interestingly, both strains originated from different clones.

Discussion
In this small study, we explored the genetic determinants of the high azithromycin and ceftriaxone MICs identified in Neisseria isolates from an exploratory study among Belgian MSM. The ceftriaxone MICs could be fairly accurately predicted by known RAMs in penA. The high prevalence of these RAMs in N. subflava is a cause of concern as the HGT of these RAMs from commensal Neisseria has been shown to be responsible for the emergence of cephalosporin resistance in the past [9]. Although known macrolide RAMs were

Discussion
In this small study, we explored the genetic determinants of the high azithromycin and ceftriaxone MICs identified in Neisseria isolates from an exploratory study among Belgian MSM. The ceftriaxone MICs could be fairly accurately predicted by known RAMs in penA. The high prevalence of these RAMs in N. subflava is a cause of concern as the HGT of these RAMs from commensal Neisseria has been shown to be responsible for the emergence of cephalosporin resistance in the past [9]. Although known macrolide RAMs were identified in commensals, these did not fully explain the high azithromycin MICs values in these commensals. In particular, none of the rRNA 23S RAMs known to cause high level azithromycin resistance in N. gonorrhoeae were found [26]. The addition of msr(D) to the regression models improved the model prediction from 21 to 84%.
An important limitation of this study is that we did not experimentally validate the effect of the msr(D) gene on macrolide MICs. Furthermore, we did not evaluate whether this gene could be transformed into N. gonorrhoeae or N. meningitidis. The ability of the ABC-F proteins, including msr(D), to induce macrolide resistance in a range of Gram-positive and -negative bacterial species is however well established [34][35][36][37][38]. We demonstrated that the complete msr(D) gene has been incorporated into the chromosomal DNA of nine of 11 circulating N. subflava strains. These nine isolates exhibit 100% sequence identity in this gene and the surrounding DNA to the MEGA element in S. pneumoniae [31]. These findings suggest that the msr(D) gene may have been acquired from S. pneumoniae (or other bacterial species [35,36]) via transformation or less likely a plasmid. The presence of a DUS just upstream from the integration site increases the probability that transformation was responsible [39]. Core genome MLST analysis suggests that these acquisitions in N. subflava either took place on more than one occasion or that the msr(D) has been taken up and lost in sub-lineages.
In S. pneumoniae, the MEGA element includes the full length mef(A) efflux pump, which is thought to act synergistically with msr(D). Msr(D) displaces macrolides from ribosomes and the mef(A) then expels the macrolides before they can reattach [28,36]. The truncated version of mef(A) in N. subflava is unlikely to be able to perform this efflux function. It is possible that other efflux pumps perform this function. Of note, all 11 of the N. subflava's in this study contained the K823E mutation in mtrD which has been shown to enhance the ability of the mtrCDE efflux pump to export macrolides [8]. The retained portion of the MEGA element contains the putative upstream regulatory MYLIFM sequence [34]. Various lines of evidence suggest that this sequence could enable the msr(D) gene to be induced by macrolides [29,35]. Direct experimental evidence for this effect is however only available for msr(A) and vmlA in the ABC-F gene family [28,34].
We also found evidence that another gene that confers antimicrobial resistance via interactions with the ribosome, tetM, has been taken up by a commensal Neisseria and integrated into its chromosome. This uptake and integration was originally described in 1987, where once again a Streptococcal species was thought to be the donor [40][41][42]. More recently, Fiore et al. have found tetM to be present in N. subflava. They did not, however, specify whether it was carried on a plasmid or integrated into the chromosome [5].
In addition to the acquisition of the msr(D) and tetM genes, several SNPs in commensals were identified which could be a source of RAMs in N. gonorrhoeae and N. meningitidis. In keeping with previous studies, we found that the mtr promoter region, penA, mtrD in commensals are a potential rich source of RAMs [7][8][9]. Whilst the G70S and G70A mutations in rplD have been shown to lead to reduced susceptibility to macrolides in N. gonorrhoeae, their presence in commensal Neisseria has not, to the best of our knowledge, been established before [43].
The high ceftriaxone and azithromycin MICs in our circulating Neisseria are a cause of concern. They are in all likelihood a result of the high antimicrobial consumption of the surveyed population [10,15]. The consumption of macrolides and cephalosporins is high in the general Belgian population compared to other European countries [44,45]. This has been linked to a combination of cultural factors (such as a high uncertainty avoidance index), and various structural factors that retard antibiotic stewardship campaigns [46][47][48][49]. Consumption is considerably higher in MSM attending our STI clinics than the general population [50]. One of the drivers of this is intensive screening and treatment of asymptomatic MSM for N. gonorrhoeae and C. trachomatis which results in macrolide exposures of around 12 defined daily doses per 1000 inhabitants [50]. Considering that this exposure is around 6-fold higher than resistance inducing thresholds in a range of bacteria, it is perhaps not too surprising that commensal Neisseria in our study population have acquired a range of mechanisms that enable them to withstand the macrolide selection pressure [51]. Taken in conjunction with the lack of evidence of a benefit of screening for N. gonorrhoeae and C. trachomatis in MSM, our results provide further evidence to support the reconsideration of this practice [12]. Our results also provide further motivation for the surveillance of AMR in commensal Neisseria (and other species) in populations at high risk for the emergence of AMR in N. gonorrhoeae and other bacteria [10,15].

Strain Collection and MIC Determination
Between January and May 2019, 10 MSM attending the Institute of Tropical Medicine (ITM) STI clinic with a diagnosis of anogenital gonorrhea were enrolled into this study [15]. After informed consent was obtained, oropharyngeal and anorectal swabs were taken. They were then treated with 500 mg ceftriaxone intramuscularly and 2 g azithromycin orally. The same swabs were taken 14 days later. All swabs were inoculated onto blood and modified Thayer-Martin agar and incubated in 5% carbon dioxide at 36.5 • C for 24 h.
All colonies with a morphology compatible with Neisseria were subcultured. Gram staining and oxidase test were performed, and Neisseria species were identified by matrixassisted laser desorption/ionization-time of flight mass spectrometry (MALDI-TOF-MS). Minimum inhibitory concentrations (MICs) of ceftriaxone and azithromycin were determined by E-test. EUCAST breakpoints for N. gonorrhoeae were used to define the reduced susceptibilities of all the Neisseria species in the panel (azithromycin (AZM) resistant (R) >1 mg/L (ECOFF); ceftriaxone (CRO) resistant (R) ≥0.125 mg/L) [23].
To evaluate the difference in MIC distribution in (i) pathogenic and commensal Neisseria and (ii) pre-and post-treatment isolates, a Wilcoxon signed-rank test was performed in R (version 3.6.3) (R Foundation for Statistical Computing, Vienna, Austria).
Genes were aligned by using MEGA with the MUltiple Sequence Comparison by Log-Expectation (MUSCLE) algorithm (v10.1.7) [65], except for the hyper variable sequences (porB and mtr promoter) which were aligned with MAFFT-einsi (v7.471) [66]. SNPs were manually checked based on known RAM positions in N. gonorrhoeae.
The protocol from David Eyre et al. was used to determine the number of RAMs of the 23S rRNA alleles [67]. In brief: a bwa index was made of 23S gene sequence (FA1090) to map reads with BWA-MEM followed by sorting and indexing with SAMTOOLS (v1.0) [68]. The ratio of base counts at both RAM positions (A2045 and C2611) were determined with a Python script available at https://github.com/davideyre/gc_mic_prediction_chapter (accessed on 4 March 2020). Abricate (v1.0.1) was used to screen for AMR genes on contigs with default settings (database used: NCBI AMRFinderPlus) [52,69].

Multiple Linear Regression Models of Azithromycin MIC Prediction in Neisseria
Multiple linear regression in R (version 3.6.3) (R Foundation for Statistical Computing, Vienna, Austria) was used to predict the MIC values of azithromycin by using the lm function. Log2 transformed values of azithromycin MIC were used as the outcome variable. Known RAMs for azithromycin in N. gonorrhoeae/N. meningitidis in mtrD, mtrRP, ponA rplD were used as predictors in the models. Two models were constructed-one excluding and one including a binary variable representing the absence or presence of msrD in the isolate, respectively. The adjusted R squared of both models were compared.

Data Availability Statement:
The data presented in this study are openly available in Genbank at https://www.ncbi.nlm.nih.gov/sra/PRJNA703317, BioProject accession number PRJNA703317.

Conflicts of Interest:
The authors declare no conflict of interest.