Molecular genetic typing techniques revealing intraspecies variability in pathogenic bacteria have been vigorously used in epidemiological research since the early 21st century [1]. Genotyping is especially important for the surveillance of multidrug-resistant infections, such as those caused by Neisseria gonorrhoeae. The World Health Organization (WHO) has included N. gonorrhoeae in the list of 12 pathogens that pose a global threat and require urgent development of novel antimicrobial drugs [2].
Surveillance of epidemiologically significant clones of N. gonorrhoeae, including multidrug-resistant isolates, and control of their global and regional spread rely on the use of Neisseria gonorrhoeae multi-antigen sequence typing (NG-MAST), Multilocus sequence typing (MLST) [3, 4] and whole-genome sequencing (WGS) [5]. WGS provides the most comprehensive information on the phylogeny of clinical isolates, but high costs and strict requirements for the quality of DNA samples still limit its application in routine practice. MLST was originally developed to type N. meningitidis and this method analyzes genes that are more conserved than porB and tbpB targeted by NG-MAST. Thus, NG-MAST currently remains the main and the most widely used technique to study the evolution of the pathogen and identify its transmission routes [4, 6]. Moreover, NG-MAST has demonstrated high resolution in the analysis of clinical isolates [7, 8].

NG-MAST consists in sequencing variable regions of two genes encoding the transmembrane porin protein (porB) and the transferrin-binding protein (tbpB). Based on sequencing data, each unique sequence is assigned a reference allele number, and the combination of allele numbers is assigned a sequence type (ST). The list of the identified alleles is constantly expanding. In January 2020, the database at www.ng-mast. net contained over 11,000 porB alleles and over 2,900 tbpB alleles that formed over 19,500 NG-MAST types.
The aim of this study was the NG-MAST typing of modern N. gonorrhoeae clinical isolates collected in Russia and the comparative phylogenetic analysis of bacteria causing gonococcal infection in Russia, European countries and Japan.

METHODS

Clinical isolates of N. gonorrhoeae (822 specimens) collected in 17 Russian regions under the Russian Gonococcal Antimicrobial Surveillance Programme (RU-GASP) were shipped to the State Research Center for Dermatovenerology and Cosmetology (Moscow). The isolates had been collected in 2013–2018 in 7 of 8 federal districts, except the Far Eastern district, specifically from Arkhangelsk, Astrakhan, Bryansk, Irkutsk, Kaluga, Novosibirsk, Omsk, Penza, Pskov, Ryazan, Tomsk, Chelyabinsk and Stavropol regions, the city of Moscow, the Republics of Tatarstan, Tyva and Chuvash. The specimens came from specialized healthcare facilities for dermatology and venerology (1 specimen per patient). The number of the obtained clinical isolates varied across the regions between 1 and 10% of the total reported cases of gonococcal infection in a given region (5 to 35 strains a year), depending on the population density and the incidence of the infection. Sample collection, shipment, culture, verification and storage conditions are described in detail in [9–11]. The number of clinical isolates analyzed during each year of our study are provided in tab. 1.
Molecular typing of N. gonorrhoeae consisted in sequencing the variable regions of the porB (490 bp) and tbpB (390 bp) genes according to the NG-MAST protocol [12]. After the first round of PCR amplification and purification of PCR products, the resulting DNA fragments were Sanger-sequenced using a 3730xl Genetic Analyzer (Applied Biosystems; USA). The obtained allele sequences were compared to the reference sequences from the NG-MAST database at www.ng-mast. net in order to identify a sequence type of each clinical isolate. Previously uncharacterized allele sequences detected in porB and tbpB or their combinations were submitted to the aforementioned NG-MAST database and assigned an NG-MAST sequence type.

Data on N. gonorrhoeae circulating in the European Union were retrieved from EUROGASP database at https://pathogen. watch/collection/eurogasp2013. We analyzed a total of 1,071 isolates from 21 EU countries, including Austria, Belgium, UK, Hungary, Germany, Greece, Denmark, Ireland, Iceland, Spain, Italy, Cyprus, Latvia, Malta, Netherlands, Norway, Portugal, Slovakia, Slovenia, France, and Sweden. NG-MAST data on 206 isolates of N. gonorrhoeae circulating in Japan (2015) were retrieved from the database at https://pubmlst.org/ (we selected samples of N. gonorrhoeae with known variants of porB and tbpB). Isolates with a mixed or unknown NG-MAST type were excluded from the analysis. A total of 29 isolates were excluded, which amounts to 1.4% of the initial sample size (EU + Japan).
Phylogenetic analysis was performed on the concatenated sequences of porB and tbpB. Model selection was done using Bayesian and Akaike information criteria. The best substitution model for our dataset was GTR with the calculation of the proportion of invariant sites and the G model of rate heterogeneity. Based on the obtained sequences and the chosen model, a maximum likelihood phylogenetic tree was constructed in RAxML ver. 7.4.8 [13]. The phylogenetic tree was created using the software at http://galaxy-dev.cnsi.ucsb.edu/ osiris/. The number of iterations applied was 999. The tree was constructed from the data on Russian and European clinical isolates of N. gonorrhoeae collected in 2013. For phylogenetic tree clades, bootstrap values of 90% or above were considered statistically significant.
The data on the population size and the proportion of isolates representing unique sequence types in each country were analyzed using Ward’s agglomerative hierarchical clustering procedure in the «cluster» package for R. The optimal number of clusters was computed in the «Nbclust» package for R [14]. Information on the population size in different countries and the net migration rate (number of migrants per 1,000 population) were taken from The World Factbook: https://www. cia.gov/library/publications/the-world-factbook/.

RESULTS

Based on the results of molecular typing, 822 clinical isolates of N. gonorrhoeae collected in Russia between 2013 and 2018 were attributed to 301 different NG-MAST types. The most prevalent sequence types were 807, 228, 1993, 5714, and 9476 (8.3%, 3.3%, 3.2%, 3.2%, and 2.7% of the total sample size, respectively). In 2013, isolates of NG-MAST type 807 made up over 13% of the Russian N. gonorrhoeae population; in the next few years, this number gradually declined, reaching 4.6% in 2018 (tab. 1). Representatives of sequence types 807 and 1993 were also detected among European specimens, but occurred there sporadically: type 807 was found in the Spanish and Slovakian populations, whereas type 1993, in Denmark. In 2018, sequence type 228, which is unique to Russia, was the most common type in Russia, amounting to 14.6% of the gonococcal population. The proportion of Russian isolates classified as unique Russian sequence types (marked by an asterisk in tab. 1) remained constant relative to the total number of isolates analyzed between 2013 and 2018 (80% on average).
Notably, the proportion of sequence types represented by only one specimen was high and stable over the entire studied period, amounting to 25–31% of the total number of isolates. In the European population of N. gonorrhoeae, it was as high as 24% in 2013 [6].
Results of N. gonorrhoeae molecular typing in European countries [5, 6] and the most common sequence types are shown in tab. 2. In 2013, 389 sequence types were circulating in the EU, the most prevalent being sequence type 1407, which occurred in 13 countries and constituted 7.3% of the total analyzed clinical isolates.
The Japanese population of N. gonorrhoeae (tab. 2) was very different from both Russian and European populations. Japanese isolates represented 65 molecular types, of which only 2 sequence types were present in Russia and European countries: NG-MAST 1407 (found in EU and Russia) and NG-MAST 4186 (one specimen from Sweden). The rest 63 molecular types were unique to Japan.
In order to establish a phylogenetic relationship between Russian and European clinical isolates of N. gonorrhoeae, a maximum likelihood phylogenetic tree was constructed. The phylogenetic tree colored according to the presence of NG-MAST type in each country is shown in fig. 1.
The analysis of sequence types distribution on the tree for the isolates from Russia and EU identified 16 clades with bootstrap values higher than 90% (fig. 1). Some clades fully or partially corresponded to the European genogroups established for N. gonorrhoeae samples collected in Europe in 2013 in the report of the European center for disease prevention and control [6].

Clade α (bootstrap value = 100) comprises 40 NG-MAST-types (82 isolates) from 17 European countries and Russia; of them 6 sequence types are represented by 7 Russian isolates (8.5%). The clade also includes sequence type 5624 (represented by 15 specimens from 7 countries) and corresponds to the previously described European genogroup G5624 [6].

Clade β (bootstrap value = 94) is formed by 33 sequence types (91 isolates), including type 225 commonly found in Europe (24 specimens from 12 European countries) and type 292 (15 specimens from 6 countries; of those specimens, one comes from Russia). Russian contribution to this clade is represented by 19 isolates (21.3%) belonging to 12 different sequence types. Clade β corresponds to the European genogroup G225 [6].

Clade Υ (bootstrap value = 96) is constituted by 14 sequence types (98 isolates) collected in 17 European countries (1 isolate comes from Russia), of them 78 isolates represent sequence type 2992, which is the second most common NG-MAST type in Europe. Clade Υ corresponds to the European genogroup G2992 [6].

Clade δ (bootstrap value = 98) comprises 3 sequence types (Portugal, Norway, UK). Clade ε (bootstrap value = 100) consists of 4 sequence types (4 isolates from Norway, Greece and UK). Clade ζ (bootstrap value = 100) contains only sequence types unique to Russia (7 specimens representing 5 NG-MAST types).

Clades η and θ (bootstrap values 98 and 90, respectively) contain only European sequence types not found in Russia. Clade η is formed by 9 NG-MAST types (16 samples), clade θ consists of 11 sequence types (34 samples). NG-MAST-types of these 2 clades represent the European genogroup G5333 [6].

Clade ι (bootstrap value = 98) consists of 4 sequence types (4 samples collected in Portugal, Italy and the Netherlands).

Clade Κ (bootstrap value = 94) comprises 4 European sequence types (36 samples), including sequence type 4995 (the fourth most common type in Europe found in 10 European countries) represented by 31 isolates. Clade Κ corresponds to the European genogroup G4995 [6].

Clade λ (bootstrap value = 95) includes 3 sequence types from European countries represented by 10 samples, of which 8 represent NG-MAST type 5441 detected in 5 countries.

Clade μ (bootstrap value = 94) comprises 4 sequence types (10 samples from Greece, Denmark, Slovenia, and Portugal). Clade ν (bootstrap value = 98) is formed by 3 sequence types (5 samples from 5 European countries).

Clade ξ (bootstrap value = 100) consists of 3 sequence types (4 samples from Ireland and Denmark). Clade ο (bootstrap value = 98) is formed by 5 NG-MAST types (17 samples), including types 1993 and 5714 that are frequently found in Russia.

Clade π (bootstrap value = 97) includes 7 sequence types (12 isolates) found both in Russia and European countries. Those include sequence types 5792, 9485, 9490, and 9491, all of which are unique to Russia.

In each analyzed country there were isolates representing unique sequence types not found in any other country. For European countries, the proportion of isolates representing unique sequence types (relative to the total number of isolates in a given country) varied between 25 and 56%. The lowest proportion was observed in the UK and Belgium, the highest, in Austria, Slovenia and Sweden (> 50%). In Russia and Japan, the proportion of unique sequence type isolates exceeded 80% (tab. 3).

To identify groups with similar distribution patterns for unique sequence types and their relationship with the population size in a country of interest, cluster analysis was carried out. Twenty-one countries were included in the analysis (tab. 3). Cyprus and Iceland were excluded due to the small number of isolates available (6 and 5, respectively). After applying Ward’s agglomerative hierarchical clustering procedure, the optimal number of clusters was determined in NbClust for R using the majority rule. Nine out of 30 methods (Silhouette, Duda, PseudoT2, Beale, Ratkowsky, PtBiserial, McClain, Dunn, SDindex) showed that the optimal number of clusters was 2 (bootstrap value = 0.01); 6 methods (Hartigan, Scott, Marriot, TrCovW, TraceW, Ball) showed that the optimal number of clusters was 3. Following the majority rule, it was decided that two clusters were the optimal choice; at the same time there was a statistically based rationale for breaking down cluster 2 into 2 sub-clusters (fig. 2А).

Cluster 1 included Russia and Japan (countries with a population size of over 125 million people). Cluster 2 comprised all European countries with a population size from 0.4 to 81 million. Cluster 2 was divided into 2 subclusters: 2a and 2b. Subcluster 2a was composed of the UK, Germany, Spain, Italy, and France (countries with a population size of 47–81 million); subcluster 2b included European countries with populations below 17 million (Ireland, Norway, Denmark, Slovakia, Malta, Latvia, Slovenia, Netherlands, Belgium, Portugal, Greece, Austria, Hungary, Sweden).
In general, the greater was the population size, the more isolates of unique sequence types were present in the country, as demonstrated by the countries with a population size over 47 million included in cluster 1 and subcluster 2a (fig. 2B). This pattern was not observed for the countries with smaller population sizes constituting cluster 2b. Transition from cluster 2 to cluster 1 occurred for the countries with the population size over 125 million (Russia, Japan).
It was also interesting to analyze the contribution of migration to the distribution of unique sequence types across different countries. If only net migration rates and proportion of isolates of unique sequence types were analyzed, we were unable to find statistically significant clusters. When 3 parameters were included in the analysis (migration, population size and the proportion of unique sequence types), the agglomerative coefficient was 0.83, which is lower than 0.93, the value yielded by the analysis of 2 parameters (population size and the proportion of unique sequence types).

DISCUSSION

The comparative analysis of N. gonorrhoeae molecular typing data revealed significant differences between Russian, European and Japanese gonococcal populations. The majority of Russian isolates belonged to the sequence types that were not detected in European countries or Japan in the analyzed period. In Russia, only 3 (0.4%) isolates represented the epidemiologically significant sequence type 1407, which prevailed in Belgium, Hungary, Denmark, Spain, Norway Portugal, Slovakia, and Japan. It is worth mentioning that in 2010 this sequence type amounted to over 10% of all isolates in many European countries, including Austria, Belgium, UK, Netherlands, Spain, Italy, Portugal, Romania, Slovenia, and Poland [15–17], as well as Japan and USA [18, 19]. This sequence type poses a great threat because it carries multiple determinants of antimicrobial resistance, including the mosaic penicillin-binding protein PBP2 and a mutation of the Ala501 residue in PBP2, both of which cause resistance to cephalosporins [16, 20, 21].
In the Russian dataset, there were only 3 (0.4%) N. gonorrhoeae isolates representing sequence type 2992 and 7 (0.9%) isolates representing type 2400. Those sequence types were second and third most common types in Europe in 2013 (7.0 and 4.7%, respectively). By contrast, types 228, 5714, and 1751 were not detected in the European population, although they were highly prevalent in the Russian population.
For Poland, NG-MAST data were retrieved from [17], as they were missing in the EUROGASP database. The distribution of NG-MAST types in Poland was closer to that in European Union than in Russia; the Polish gonococcal population was dominated by sequence type 1407 (tab. 2).

The analysis of the phylogenetic tree constructed for NG-MAST types of clinical isolates collected in Europe and Russia reveals that N. gonorrhoeae populations are heterogenous both in Russia and Europe. There are sequence types that are present in many countries; at the same time, isolates coming from one and the same country can belong to phylogenetically distant clades. The most common European sequence types (NG-MAST 1407, 2992, 2400, 4995, 21, 225) that amount to the total of 25.9% of all European strains are phylogenetically distant from each other and distributed throughout the entire tree. Notably, there is no close relationship between the most common European sequence types and Russian sequence types, including NG-MAST 807, 1152 and 5941. Sixteen clades of the tree are characterized by high bootstrap values and partially or fully correspond to the European genogroups [6].
Further phylogenetic analysis of NG-MAST types drove us to the conclusion that geographically distant populations of N. gonorrhoeae evolve more or less locally. This could be one of the reasons underlying the difference in the incidence of gonococcal infection in Russian and EU [22]: in Russia, the incidence of the infection is declining [9, 23] and there are no isolates resistant to ceftriaxone (the primary drug for gonorrhea treatment) [9, 10].
We have analyzed the associations between the proportion of isolates of unique sequence types and the population size of the country of interest. In Russia and Japan, the majority of isolates (> 80%) belonged to NG-MAST types unique to this country. The cluster analysis of the population size and the proportion of isolates belonging to unique sequence types detected in a country allowed us to identify two clusters: the first cluster included Russia and Japan (population over 125 million people), the second cluster included European countries with population between 0.4 and 81 million people. On the whole, there was a certain distribution pattern for unique sequence types: the greater the country’s population, the higher the proportion of samples with unique sequence types.

CONCLUSIONS

NG-MAST typing of N. gonorrhoeae clinical isolates collected in Russia has demonstrated their significant difference from the gonococcal populations circulating in the EU and Japan. The Russian N. gonorrhoeae population was dominated by sequence types 807, 1152 and 5941, which occurred only sporadically in other countries; sequence types 228, 5714 and 1751 were not detected in Russia at all. Between 2013 and 2018, in Russia there were only 3 isolates of epidemically significant sequence type 1407 (0.4% of the total samples under study). This type poses a serious threat because it carries multiple determinants of antimicrobial resistance and prevails in many European countries and Japan. The phylogenetic analysis of NG-MAST types for Russian and European isolates has demonstrated their high heterogeneity and genetic distance between common European and common Russian NG-MAST types, suggesting that the Russian population of N. gonorrhoeae has been forming and evolving locally.
The majority (> 80%) of Russian and Japanese isolates are unique to the two countries. This number is higher than the proportion of isolates with unique sequence types in European countries. The cluster analysis of the proportion of isolates representing unique sequence types and the population size allowed us to identify two clusters: one cluster for Russia and Japan and another cluster for European countries. The following trend was observed for the countries with a total population size over 47 million (Spain, UK, Germany, Italy, France, Russia, Japan): the greater the population size, the higher the proportion of isolates representing unique sequence types. We were unable to establish an association between the proportion of isolates representing unique sequence types and the net migration rate.
Thus, the phylogenetic analysis has revealed a relative isolation of the currently existing Russian population of N. gonorrhoeae, which follows its own evolutionary patterns. Nevertheless, there is a need for continuous surveillance of the spread of gonococcal infection and antimicrobial resistance of N. gonorrhoeae in Russia. The objective of this surveillance is timely detection and effective elimination of globally spreading sequence types with multiple determinants of antimicrobial resistance.