Relationships between Virulence Genes and Antibiotic Resistance Phenotypes/Genotypes in Campylobacter spp. Isolated from Layer Hens and Eggs in the North of Tunisia: Statistical and Computational Insights

Globally, Campylobacter is a significant contributor to gastroenteritis. Efficient pathogens are qualified by their virulence power, resistance to antibiotics and epidemic spread. However, the correlation between antimicrobial resistance (AR) and the pathogenicity power of pathogens is complex and poorly understood. In this study, we aimed to investigate genes encoding virulence and AR mechanisms in 177 Campylobacter isolates collected from layer hens and eggs in Tunisia and to assess associations between AR and virulence characteristics. Virulotyping was determined by searching 13 virulence genes and AR-encoding genes were investigated by PCR and MAMA-PCR. The following genes were detected in C. jejuni and C. coli isolates: tet(O) (100%/100%), blaOXA-61 (18.82%/6.25%), and cmeB (100%/100%). All quinolone-resistant isolates harbored the Thr-86-Ile substitution in GyrA. Both the A2074C and A2075G mutations in 23S rRNA were found in all erythromycin-resistant isolates; however, the erm(B) gene was detected in 48.38% and 64.15% of the C. jejuni and C. coli isolates, respectively. The machine learning algorithm Random Forest was used to determine the association of virulence genes with AR phenotypes. This analysis showed that C. jejuni virulotypes with gene clusters encompassing the racR, ceuE, virB11, and pldA genes were strongly associated with the majority of phenotypic resistance. Our findings showed high rates of AR and virulence genes among poultry Campylobacter, which is a cause of concern to human health. In addition, the correlations of specific virulence genes with AR phenotypes were established by statistical analysis.


Introduction
The most prevalent cause of bacterial gastroenteritis, Campylobacter spp., accounts for 5-14% of all diarrheal diseases worldwide [1]. Among the human-associated Campylobacter species, 95% of campylobacteriosis is caused by the C. jejuni and C. coli species [2], causing 96 million cases of diarrhea each year globally. In addition, Campylobacter is the second most prevalent agent of diarrhea in Europe (after Salmonella) [3]. Campylobacter was the second aetiological agent of outbreaks related to food and water poisoning in 2018 [4]. Contrary to European and developing countries, there are few reports of human campylobacteriosis in North African countries, including Tunisia, presumably owing to the low prevalence of the disease or the sporadic cases of infections. Globally, high rates of antimicrobial resistance (AR) have been increasingly noted, which dramatically reduced treatment options for campylobacteriosis. The National Antimicrobial Resistance This study sought to determine whether specific virulence genes, resistance genes, and AR characteristics were associated with one another in Campylobacter isolates. To achieve this, we determined the antimicrobial susceptibility and investigated by PCR-selected genes of virulence and AR in a collection of Campylobacter isolates collected from laying hens and eggs. The relationship between the different aforementioned traits was then investigated using a variety of statistical and computational methodologies.

Ethics Statement
The Biomedical Ethics Committee of the Institut Pasteur de Tunis gave its approval to this study, and the sampling protocol was performed according to internationally recognized guidelines ISO 10272-1:2006 (Annex E) for the detection of Campylobacter spp. [33].

Bacterial Strains
One hundred seventy-seven Campylobacter isolates have been reported previously [34]. These isolates include 124 C. jejuni and 53 C. coli recovered from five laying hen farms located in the northeast of Tunisia between October 2017 and May 2018.

PCR Detection of Genes Encoding AR
Fluoroquinolone resistance is commonly encoded by single point mutation (Thr-86-Ile) in the quinolone resistance-determining region (QRDR) of the GyrA subunit of the DNA gyrase enzyme [37]. For C. jejuni isolates, MAMA-PCR was performed as previously reported [38], while for C. coli, the used protocol was as cited by Zirnstein et al. [37]. MAMA-PCR was also used to detect point mutations at positions 2074 and 2075 in domain V of the 23S rRNA gene, which are related to erythromycin resistance, as described previously [39]. For all isolates, the following genes were detected by the classical PCR method: erm(B) (erythromycin resistance) Qin et al. (2014), tet(O) (tetracycline resistance), aph-3-1 (aminoglycosides resistance), cmeB (multidrug efflux pumps), and bla OXA-61 (beta-lactam resistance) (Table A2). Positive control strains from our collection were used in every PCR analysis [36].

Statistical Analysis
Statistical analysis was performed to investigate a possible association between virulence genes and AR in all isolates. We studied the antimicrobial susceptibility phenotypes (resistance/susceptibility) against the eight tested antibiotics (Amp, Amc, Cip, Nal, Ery, Tet, Chl, and Gen), and associated the latter with the presence/absence of the investigated virulence genes (cadF, ciaB, racR, flaA, dnaJ, cdtA, cdtB, cdtC, virB11, pldA, wlaN, ceuE, and cgtB). This was investigated first for all Campylobacter isolates, and then for the isolates of each species. The association test of each virulence gene with the AR phenotype was computed by Pearson's chi-square or Fisher's exact test using R software via RStudio (version 1.4.1103). Fisher's exact test was used when the expected cell counts for the contingency table held less than five isolates. If the p-value < 0.05, the association was deemed statistically significant.

Network Generation
Two groups of networks were built connecting phenotypical AR with virulence genes, as well as AR genes with virulence genes. The networks were displayed via Cytoscape (https://cytoscape.org/) (20 February 2022) (version 3.8.1) (https://pubmed.ncbi.nlm.nih. gov/14597658/) (20 February 2022). These networks were built with the aim of revealing co-occurrence patterns and identifying interactions that could reveal information on the patterns of the incidence of virulence genes and AR across all Campylobacter isolates (only virulence genes that showed a statistically significant association were used to build the network).

Predictive Analysis Using Machine Learning Random Forest Algorithm
Following the statistical association test, a predictive analysis was performed using the machine learning Random Forest algorithm, via the randomForest R package (https: //link.springer.com/article/10.1023%2FA%3A1010933404324) (20 February 2022), in order to determine the most important virulence genes that could be related with a specific AR phenotype. Classification trees are used in the analysis to establish, for each variable, its importance in classifying the data and determining the outcome through the production of an importance score [40]. Only virulence genes that showed a statistically significant association with AR through Pearson's chi-square/Fisher's exact test for all Campylobacter isolates (both C. coli and C. jejuni species together) were considered for this classification. The Random Forest measures the contribution of each virulence gene to a particular resistance phenotype. The algorithm produces a MeanDecreaseGini score that gives a valuable estimation of the significance of the variable in the model and thus, in our case, valuable information to determine which gene is more likely to be linked to an increased probability of a specific AR [41].

Virulotypes and Phenotypic Profiling of AR
One-hundred-and-seventy-seven Campylobacter isolates (124 C. jejuni and 53 C. coli) were analyzed to determine the virulotype (content of genes encoding virulence factors) and phenotypic AR profiles. All isolates (n = 177, 100%) harbored the flaA, cadF, ciaB, and cdt genes, closely followed by the racR gene (n = 161, 90.96%) ( Figure 1A). A close result was obtained when analyzing the 124 C. jejuni isolates. Indeed, the flaA, cadF, ciaB, and cdt genes were present in all isolates (100%), followed by the dnaJ (n = 119, 95.97%) and ceuE (n = 115, 92.74%) genes ( Figure 1B). There were no discernible differences found in the C. coli species for the most common virulence genes. Indeed, all isolates contained the flaA, cadF, racR, ciaB, and cdt genes, whereas, the pldA gene was detected in 51 (96.22%) isolates. Interestingly, a major difference was observed concerning the ceuE gene, which was absent in all C. coli isolates but highly present in the C. jejuni ones (92.74%) ( Figure 1C).

Molecular Detection of AR Genes
All isolates (n = 177) carried the tet(O) and cmeB genes, according to the PCR data. In the β-lactam-resistant C jejuni and C coli isolates, the bla OXA-61 gene was found in 18.82% and 6.25%, respectively. For the quinolone-resistant isolates, the Thr-86-Ile mutation in GyrA was found in all C. jejuni and C. coli isolates. Similarly, all erythromycin-resistant isolates harbored the A2075G and A2074C mutations, while the erm(B) gene was detected in 53.71% (94/175) of the erythromycin-resistant Campylobacter isolates, being in 60 (48.38%) and 34 (60.15%) of the erythromycin-resistant C. jejuni and C. coli isolates, respectively. There was no isolate harboring the aphA-3 gene.

Statistical Analysis of Phenotypic AR with Virulence Genes 214
Pearson's chi-square and Fisher's exact tests were executed to study the association 215 between the set of virulence genes and AR in all isolates showing resistance to 4-6 anti-216 biotics, as well as those resistant to more than six antibiotics (  Table 2). 223 However, no significant relationships were observed concerning the C. coli isolates (Table   224 3).

Network Analysis of Resistance, Virulence Genes, and Phenotypic AR
In order to examine the co-occurrence patterns, we generated networks describing the connections between (i) phenotypic AR with virulence genes and (ii) AR genes with virulence genes to provide information on the patterns and incidence of virulence genes and AR across all Campylobacter isolates. Figure 3 reveals three distinct networks that describe links between AR and the presence/absence of certain virulence genes for each isolate.
Focusing only on the virulence genes that showed a statistically significant association, we noticed the coexistence of some connections between phenotypic AR and specific virulence genes among some isolates more frequently than other ones. Approximately, for a third of the isolates (n = 50), there was a high frequency of connections linking nalidixic acid (Nal), tetracycline (Tet), erythromycin (Ery), ciprofloxacin (Cip), ampicillin (Amp), and chloramphenicol (Chl) resistance with the virulence genes pldA and racR ( Figure 3A). Similarly, when looking into the networks generated for 100 antimicrobialresistant isolates and the total number of Campylobacter isolates (n = 177), the same highfrequency connections were created between phenotypic AR and the virulence genes pldA and racR as shown in Figure 3B,C, respectively.
there was a high frequency of connections linking the following resistance genes: cmeB, 258 tet(O), Cj-gyrA, and 23S rRNA (mutated) with the virulence genes ceuE, pldA, and racR 259 and the ermB with pldA and racR ( Figure 4A). However, for 100 and 177 isolates, the latter  Figure 3. Visualization of the co-occurrence pattern of phenotypic AR and virulence genes across Campylobacter isolates (only virulence genes that showed a statistically significant association were used). Red lines indicate a high incidence of links between AR and virulence among isolates. The line thickness between the nodes reveals the frequency of isolates with identical coincident connections. Nodes in orange and blue are AR and virulence genes, respectively. (A) Connections between phenotypic AR and virulence genes across 50 Campylobacter isolates out of 177 isolates. (B) Connections across 100 Campylobacter isolates out of 177 isolates. (C) Connections across all Campylobacter isolates (n = 177). Amp, ampicillin; Amc, amoxicillin/clavulanic acid; Cip, ciprofloxacin; Nal, nalidixic acid; Ery, erythromycin; Tet, tetracycline; Gen, gentamicin, and Chl, chloramphenicol.
In Figure 4, we displayed three networks that show the connections between resistance genes and the virulence genes for each Campylobacter isolate. For 50 isolates, there was a high frequency of connections linking the following resistance genes: cmeB, tet(O), Cj-gyrA, and 23S rRNA (mutated) with the virulence genes ceuE, pldA, and racR and the ermB with pldA and racR ( Figure 4A). However, for 100 and 177 isolates, the latter connections were conserved by adding new connections linking ermB with the virulence gene ceuE. New added links have shown a high frequency of connection between bla OXA-61 with racR, pldA, and cgtB and Cc-gyrA with racR and pldA ( Figure 4B,C). specific AR, Random Forest produces a MeanDecreaseGini value, and the higher this 280 value is, the higher the significance of the variable in the model. 281 This investigation showed that one virulence gene, racR, displayed the most important value with 282 two antibiotics, nalidixic acid, and ciprofloxacin ( Figure 5C,D). On the other hand, another gene,

Predictive Analysis of AR/virulence Genes Links Using the Machine Learning Random Forest Algorithm
The Random Forest algorithm was used to further explore the possible association of the virulence genes that showed a significant association upon the statistical analysis for all Campylobacter isolates. In order to predict which one could be the best indicator of a specific AR, Random Forest produces a MeanDecreaseGini value, and the higher this value is, the higher the significance of the variable in the model.
This investigation showed that one virulence gene, racR, displayed the most important value with two antibiotics, nalidixic acid, and ciprofloxacin ( Figure 5C,D). On the other hand, another gene, ceuE, has shown the most important value with five other antibiotics, Amoxicillin, Erythromycin, Tetracycline, Chloramphenicol, and Gentamicin ( Figure 5A,E-H). Finally, the pldA gene showed an important value for Ampicillin only ( Figure 5B).

Antimicrobial Resistance and Corresponding Genotypes
The treatment of Campylobacter infections is currently jeopardized by the emergence of AR, which has become a complex challenge and a major issue for global public health. The Tunisian government lacks an integrated program for monitoring AR in primary human and production animal pathogens such as C. jejuni, C. coli, and C. fetus, making it difficult to implement new antimicrobial control and restriction measures. Furthermore, unlike other European countries, Tunisia has no specific legislation mandating campylobacteriosis testing. AR studies are thus critical for characterizing the circulating Campylobacter strains in Tunisian poultry. Mobile genetic elements, including plasmids and transposons, which can also carry virulence determinants, are highly associated with the global spread of AR. In Tunisia, research on AR in Campylobacter isolates from laying hens and eggs is scarce. Thus, herein, we analyzed 177 Campylobacter isolates (124 C. jejuni and 53 C. coli) from the layer hens and eggs collected in the north of Tunisia. The isolates were investigated to determine their virulotypes and AR phenotypes.
This research revealed no discernible differences in the status of certain virulence genes between Campylobacter isolates that are resistant to 4-5 antibiotics or to more than 6 antibiotics. Our findings revealed that resistance to erythromycin, tetracycline, quinolone, and ciprofloxacin is common, which can considerably restrict the number of the available treatment options of infections caused by such strains. High rates of resistance are anticipated because these antibiotics have been on the market for a long time and have been used widely in both legal and illegal situations.
The interaction(s) between AR and virulence is still poorly understood. However, there is strong scientific proof that the development of AR by the overexpression of genes encoding AR or multidrug-resistant efflux pumps causes a fitness cost to bacteria, such as lower growth rates and pathogenicity [39,42]. However, many other studies have found that pathogens' acquisition of AR improves their fitness and virulence [30,43].
The majority of our isolates were resistant to tetracycline, ciprofloxacin, and nalidixic acid. Several other studies, especially recent ones, have revealed similar high rates of resistance [35,44]. Indeed, the selection and development of antimicrobial-resistant Campylobacter are enhanced by the widespread use of these antibiotics in the treatment, management, and disease prevention in livestock.
In the majority of isolates, AR phenotypes corroborate well with the presence of genes and genetic mutations encoding AR. Tetracycline resistance has been linked to the tet(O) gene encoding the ribosomal protection protein TetO, which is commonly detected in a variety of Gram-positive and Gram-negative bacteria [45,46]. Furthermore, tetracycline is overused in the avian industry because of its low cost and simplicity of administration through drinking water [47]. It is worth noting that the chicken's body temperature (42 • C) promotes conjugation and thus contributes to the sharing of plasmids carrying various AR genes [48].
Our isolates showed a high resistance rate to fluoroquinolones. The cmeABC operon, encoding multidrug efflux, is the major molecular cause of this resistance in Campylobacter [44]. This operon was detected in all our isolates independently of their resistance or susceptibility to quinolones/fluoroquinolones. The second resistance pathway involves one or more point mutations in the QRDR of the GyrA protein, namely the Thr-86-Ile substitution, which is frequently observed in quinolones/fluoroquinolones-resistant isolates [49]. The widespread use of a specific fluoroquinolone (enrofloxacin) in avian industries has caused the selection and wide dissemination of resistant Campylobacter strains, which explains the rising resistance trend globally [50]. The world health organization (WHO) classified fluoroquinolones-resistant Campylobacter strains as high-priority pathogens resistant to antibiotics, requiring the development of new antibiotics [51,52].
Similarly, numerous studies have shown that the misuse of macrolides in poultry production has resulted in high rates of macrolide resistance in avian Campylobacter strains. All erythromycin-resistant Campylobacter isolates had the two-point mutations A2075G and/or A2074C in the gene encoding 23S rRNA [49]. The erm (B) gene, which can be carried by a variety of multi-drug resistance gene islands (MDRGI), was found in 53.1% [94/177: 48.38% (60/124) C. jejuni and 64.15% (34/53) C. coli] of our erythromycin-resistant Campylobacter isolates. Since the discovery of the erm(B) gene in Campylobacter in China, it has also been detected in turkey isolates in Spain in 2016 [53] and in the United States in 2016 in a human who previously visited Malaysia [54]. Being found on MDRGIs alongside resistance genes to other antimicrobials including ampicillin, ciprofloxacin, and tetracycline makes noteworthy the presence of this gene in our Campylobacter isolates [6].
Since macrolides, in particular erythromycin and azithromycin, are the preferred antibiotics for treating human Campylobacter infections, these findings are worrisome.
Campylobacter spp. is intrinsically resistant to beta-lactam antibiotics, including ampicillin [55]. However, acquired resistance has been reported. Indeed, enzymatic inactivation by the beta-lactamase encoding gene bla OXA-61 , detected in 18.82% and 6.25% of β-lactam resistant C jejuni and C coli isolates, respectively, is the main mechanism of acquired ampicillin resistance; in addition, other molecular mechanisms such as porins and the reduced affinity of penicillin-binding protein (PBP) have also been reported [55,56]. The majority of our isolates were gentamicin-susceptible, which is in agreement with previous reports [57][58][59]. This might be linked to its limited use for systemic infections [60], and it is not used in poultry production [58].

Virulence Power of Campylobacter Isolates
The virulome of the Campylobacter species contributes to their pathogenicity [61], hence the virulence factors of avian Campylobacter need to be investigated for consumer safety. All our isolates had the flaA, cadF, and ciaB genes, which are related to adhesion, colonization, and invasion, as well as the cdtA, cdtB, and cdtC genes, which are critical for CDT expression. The detected frequencies of these genes were analogous to those reported previously from Korea [62], Poland [63], and Italy [64], but higher than those reported from South Africa and Chile [65,66]. The presence of the cadF and ciaB genes promotes Campylobacter adhesion and internalization in cell models [47,67]. The pldA gene encoding the outer membrane phospholipase A was detected at a higher rate in C. jejuni than in C. coli, which is consistent with findings from South Africa [59], Japan [68], and Iran [69]. In addition, regardless of species, all Campylobacter isolates contained the virB, racR, and dnaJ genes.

Relationship between Virulence Genes and Phenotypic and Genotypic Antimicrobial Resistance
We identified a possible link between virulence genes and antibiotic resistance by analyzing the antibiotics to which Campylobacter isolates are more resistant or susceptible. Interestingly, using Pearson's chi-square and Fisher's exact tests, the virulence genes racR, pldA, CeuE, and cgtB were found to be closely associated with MDR Campylobacter isolates. The same analysis was also performed for each species. C. jejuni isolates showed a significant relationship between AR and the different virulence genes, specifically for racR, virB11, pldA, and cgtB. The racR, pldA, and virB11 genes facilitate bacterial adhesion and intracellular invasion [63]. In addition to the above-mentioned virulence genes linked to Campylobacter adhesion and invasion, the ceuE gene is one of the four most significant predictor genes in resistant Campylobacter isolates. Interestingly, the cgtB gene is also a significant gene that has demonstrated a substantial correlation between overall antibiotic resistance status and the prevalence of virulence genes. It is thought to play an important role in the manifestation of Guillain-Barré syndrome, the most severe side effect of human Campylobacter infection [70,71]. Since the cgtB gene enables bacteria to survive certain stressors, it can also be predicted to be associated with increased AR. However, no significant relationship was observed for C. coli, correlating with previous findings [36].
The co-occurrence network demonstrated three distinct networks that illustrate the links between phenotypic AR and the presence or absence of certain virulence genes in each isolate. We observed the coexistence of certain connections between AR and specific virulence genes among the isolates more frequently than others when we focused only on the virulence genes that showed a statistically significant association. There was a high frequency of connections linking nalidixic acid, tetracycline, erythromycin, ciprofloxacin, ampicillin, and chloramphenicol resistance with the virulence genes pldA and racR in nearly one-third of isolates (n = 50). Similarly, when the networks for resistant isolates (n = 100) and the total isolates (n = 177) were examined, the same high-frequency connections were observed between phenotypic AR and the virulence genes (pldA and racR).
When we looked at the relationship between virulence genes and AR using various approaches, we noticed that our network visualization matches the Random Forest analysis. We used the Random Forest approach to forecast the value of each virulence gene in order to figure out which gene is more significant for increasing the likelihood of phenotypic resistance in Campylobacter isolates. The virulence genes racR and ceuE were revealed to be the most important predictors of phenotypic resistance in Campylobacter. Finally, only one antibiotic, ampicillin, has proven significant value for the pldA gene.

Conclusions
Using statistical and computational tools, we demonstrated the relationship between the distribution of bacterial virulence genes and their phenotypic AR pattern and AR genes among Campylobacter isolates from layer hens and eggs. Furthermore, we have shown that the virulence genes racR, pldA, virB11, ceuE, and cgtB and the AR genes tet(O), cmeB, and bla OXA-61 , as well as mutations in rRNA 23S and gyrA, warrant further investigation using a wide range of antimicrobials to prove links that may increase virulence in bacteria. The findings of this study will be valuable in determining the association between phenotypic traits and genetic characteristics such as the status of virulence and AR genes. Our findings thus open up the possibility for further research into the pathophysiology and the underlying causes of antibiotic resistance.  Institutional Review Board Statement: The study was approved by the Biomedical Ethics Committee of the Pasteur Institute of Tunis, with reference number: 2018/12/I/LR16IPT03.

Informed Consent Statement: Not applicable.
Data Availability Statement: The datasets generated and/or analyzed during the current study are available from the corresponding author on reasonable request.