A total of 27 clinical isolates of K. pneumoniae were kindly provided by the Institute of Antibiotics, Huashan Hospital. The genomic DNA of these strains was extracted using a Bacterial Genome DNA kit (Tiangen, China), and the whole genome of each isolate was sequenced using the Illumina NovaSeq 6000 platform. Quality checking and de novo assembly were conducted using the Fastx-toolkit and Velvet (version 1.2.03), respectively (Zerbino and Birney, 2008). Draft genome annotation was performed by searching against the National Center for Biotechnology Information (NCBI) non-redundant (nr) database. The complete genomes of all sequenced strains of K. pneumoniae available in the NCBI database (as of April 2021) were downloaded, and background information on each strain was collected from the BioSample database. Strains were then further screened to select those with complete background information, including geographic location, isolation source, host disease, and collection date. In addition, K. pneumoniae strains HS11286, NTUH-K2044, Kp52.145, MGH78578, and ATCC 43816 were included in the analysis because they contained a clearly identified T6SS in the literature (Sarris et al., 2011; Liu et al., 2012; Rodrigues et al., 2020). From the initial screening, a total of 214 strains were selected for analysis in the study. Detailed information on all strains used in the study can be found in Supplementary Table 1.

The integrated database SecReT6 (Li et al., 2015) was used to identify T6SS loci. All genome sequences were analyzed by BLASTp 2.10.1 + (e-value 1e-4) against experimentally validated T6SS clusters in the SecReT6 database. BLASPn (e-value 1e-5, > 70% nucleotide identity and > 70% coverage) was used to identify T6SS-related genes in the plasmids of the 214 strains with complete genome sequences. As 27 strains were draft genomes and plasmid sequence could not be determined, the identification of T6SS-related genes in the plasmids of these strains was not analyzed.

The type i2 T6SS clusters in different species were compared using Easyfig (version 2.2.5) with default parameters (e-value 1e-3). The nucleotide sequences and annotations of Escherichia coli 042 (NC_017626), K. pneumoniae HS11286 (NC_016845), K. pneumoniae (NZ_FO834906), Enterobacter cloacae ATCC 13047 (NC_014121), and Burkholderia cenocepacia H111 (NZ_HG938371) were downloaded from the NCBI RefSeq database.

BLASTp (e-value 1e-5, > 50% amino acid identity) or RPS-BLAST was used to predict the potential effectors of K. pneumoniae T6SSs. Conserved protein domains were examined using the conserved domain database (CDD) at NCBI¹ and AlphaFold2 (Jumper et al., 2021). Putative effectors were then classified based on their amino acid sequence similarity and conserved domains. Sequences of known T6SS effector (Tse) proteins were downloaded from SecReT6 (as of September 2021). In addition, PLD (CDO15120.1) effector protein was included in the analysis (Hsieh et al., 2019). The designed method would identify putative effector proteins with >50% amino acid coverage in the analyzed strains. In addition, the functional domains of proteins in the strains were considered to avoid missing identification; thus, if a domain of strain protein was consistent with the known effectors, it would also be identified. Domain visualizations were conducted with Tbtools (Chen et al., 2020).

A homologous single-copy gene-based approach was employed to generate a phylogenetic tree of K. pneumoniae genome sequences. The STAG algorithm of OrthoFinder (version 2.5.4) was used to reconstruct a phylogenetic tree based on homologous genes (Emms and Kelly, 2019). Multiple sequence alignments of protein sequences, such as VgrG and PAAR, were performed using the Multiple Sequence Comparison by Log-Expectation (MUSCLE, version 3.8) algorithm (Edgar, 2004). The alignments were then refined using Gblock (Talavera and Castresana, 2007). FastTree was selected for the reconstruction of protein phylogenetic trees by the maximum-likelihood (ML) method. Tree visualization was performed using iTOL version 6² (Letunic and Bork, 2021).

Multilocus sequence typing (MLST) of the strains was performed with mlst version 2.19.0 (Seemann T, mlst Github)³ (Jolley and Maiden, 2010). The Comprehensive Antibiotic Resistance Database (CARD) (Alcock et al., 2020) was used for the identification and annotation of antibiotic resistance genes, while the Virulence Factor Database (VFDB) (Chen et al., 2005) and BLASTp (e-value 1e-5) were used to identify virulence factors.

Domain analysis of PAAR was conducted using the NCBI CDD⁴. The MEME online tool⁵ was applied for predicting the conserved motif structure of the PAAR family (Bailey et al., 2015).

Phyre2⁶ was utilized to predict protein structures similar to PLD superfamily protein in K. pneumoniae (Kelley et al., 2015). PLD three-dimensional structure prediction was conducted using AlphaFold2, and the predicted local distance difference test (plDDT) was used to evaluate the quality of the structure model (Jumper et al., 2021). The root-mean-square deviation (RMSD) calculated by Pymol was used to assess the similarity between protein structures (Delano, 2002). Alignments of the PLD protein sequences from Arabidopsis thaliana (NP_188194.1), Homo sapiens (AAH15033.1), K. pneumoniae (CDO15120.1), and P. aeruginosa PAO1 (AAG08474.1 and AAG06875.1) were performed using ClustalW multiple sequence alignment⁷ and colored with ESPript⁸ (Robert and Gouet, 2014). Domain analysis and visualization were performed using CDD and Tbtools, respectively.

The prevalence of T6SS loci was evaluated in 241 strains of K. pneumoniae, including strains available from the NCBI databases. Strains with background information (Supplementary Table 1) were included in the analysis. Genes encoding a T6SS cluster were present in all analyzed strains of K. pneumoniae. The majority of these isolates (229 of 241) harbor intact T6SS loci encoding conserved genes, and the vgrG gene is followed by variable regions that tend to have predictive immune-effector protein pairs (Figure 2). The T6SS gene cluster was incomplete in the remaining 12 strains. Further analysis showed that there were 1–4 T6SS loci in K. pneumoniae, and two clusters were the most common among the strains (178 of 241 strains), followed by three T6SS clusters (52 of 241 strains) (Supplementary Table 2). Numerous “orphan” genes encoding Hcp, VgrG, and PAAR were identified as being located outside the major T6SS gene clusters. This is congruent with those found in other species. Putative effectors are commonly encoded in the vicinity of these “orphan” genes (Huang et al., 2017; Ma et al., 2017a). There is no TagA protein encoded in the T6SS gene cluster of K. pneumoniae, but another protein, TagL, constitutes the T6SS protein complex and has been confirmed to be involved in the formation of the T6SS transmembrane complex in E. coli (Nguyen et al., 2021). The intact T6SS cluster in K. pneumoniae can be divided into two main types. One was represented by strain HS11286 (Figure 2), and the effector-immune protein genes appeared in clusters in variable regions, which was different from the effector-immune pairs in other bacteria, such as Acinetobacter baumannii. Another T6SS gene distribution pattern in K. pneumoniae is typical of strain NTUH-K2044, with no effector protein genes identified in the variable region and the location of PAAR differed from that of the first distribution type. In some strains, such as Kp52.145, PAAR proteins were not even identified. In addition, the T6SS cluster was not predicted in plasmids, but components of T6SS-related genes were identified (Supplementary Table 3).

Employing the schemes proposed by SecReT6 and previous reports (Barret et al., 2011, 2013; Russell et al., 2014b), T6SS clusters are classified into i, ii, or iii types according to the structure of gene clusters and the evolutionary relationship of core components. The type i class contains six subtypes (i.e., i1, i2, i3, i4a, i4b, and i5), while type ii class is endemic to the Francisella pathogenicity island-encoded system. The type iii class was identified in Bacteroidetes. The complete T6SS cluster in K. pneumoniae belongs to the type i2 class according to the data in the SecReT6 database. In addition to the T6SS cluster in K. pneumoniae, the type i2 class of T6SS has also been identified in E. coli, E. cloacae, and B. cenocepacia. The T6SS gene cluster has previously been compared in K. pneumoniae (Sarris et al., 2011; Barbosa and Lery, 2019), and in this study, the type i2 T6SS clusters in strains of other bacteria were compared with those of K. pneumoniae. K. pneumoniae HS11286 and K. pneumoniae Kp52.145 were selected for comparison because they were strains with two and three T6SS gene clusters, respectively, as reported previously (Liu et al., 2017; Storey et al., 2020). Conservation of the type i2 class of T6SS gene cluster was analyzed in E. coli 042, E. cloacae ATCC 13047, B. cenocepacia H111, and the two previously mentioned strains of K. pneumoniae. There was high conservation of the T6SS-related genes within K. pneumoniae (Figure 3), and the type i2 T6SS gene cluster in K. pneumoniae was highly similar to that of E. cloacae ATCC 13047. However, the sequence and homology of the type i2 T6SS in B. cenocepacia H111 were quite distinct from those of K. pneumoniae. The middle variable regions in each strain were different and contained genes that encode unknown functional proteins. These results indicate the diversity of effector proteins. The key genes hcp, vgrG, and paar in the T6SS were depicted in different colors in all strains. The type i2 T6SS in the five strains all included vgrG and paar genes, but that of E. cloacae ATCC 13047 did not include hcp genes. The direction and sequence of some genes in the type i2 T6SS gene cluster of B. cenocepacia H111 were inconsistent with those of the other four strains. In addition, although some strains were identified as containing a type i2 T6SS cluster, they lacked some conservative T6SS genes, e.g., strain Kp52.145 did not have a paar gene. There were also strains with only one unknown type of T6SS identified, and the function of these T6SS clusters has yet to be elucidated.

The genomes of 241 strains of K. pneumoniae, which represented both historical and present isolates from a range of geographical locations, were analyzed, with the majority (75 of 241) belonging to ST11, followed by ST147 (18), ST15 (14), ST23 (12), and ST231 (12). To identify the putative effectors in K. pneumoniae, BLASTp and RPS-BLAST were performed. Given that Hcp and VgrG proteins are both components of T6SS cluster structure and effector proteins, and exist in each analyzed strain, they are not classified in Table 1. According to the classification methods in previous studies (Russell et al., 2014a; Li et al., 2015), effector proteins predicted in this study were classified into four types, namely T6SS DNase/RNase effector (Tde), T6SS lipase effector (Tle), T6SS effector (Tse), and T6SS peptidase effector (Tpe). The effector proteins were then divided into 11 groups based on their domains (the list of effectors for each strain is supplied in Supplementary Table 4), and a number was assigned to each group for identification. Table 1 and Figure 4 show the domains and functions of each group of putative effectors. Putative functions of effectors in K. pneumoniae included nuclease (group 1, 3, 10), lipase (group 6), Mn-containing catalase (group 7), and polysaccharides hydrolytic enzymes (group 11). Among them, effector proteins in group 1 were identified as S-type pyocin (pfam06958), which are polypeptide toxins produced by bacteria. S-type pyocins cause cell death through the breakdown of DNA due to endonuclease activity. Groups 4/5/6/10 were located in the downstream variable region of vgrG in the T6SS cluster. Group 4 contained a DUF2235 domain (pfam09994) whose specific function is unknown and usually found in clusters with multiple immune proteins. Group 5 effectors contained an alpha/beta hydrolase domain of unknown function. The effectors of group 6 were a class of lipase with a PLD superfamily domain (cl15239) that contains an H-x-K-x(4)-D signature motif. These PLD effectors were identified in no more than 16 strains and may be associated with the virulence of K. pneumoniae (Lery et al., 2014). The effector protein in group 10 was PAAR, a structural component of the T6SS, that contained the PAAR_CT_2 domain (cd14744) and the C-terminal toxin domain (pfam15607) with RNase activity. The group 10 effector proteins listed in this study refer to PAAR with a C-terminal toxic domain. Both groups 4 and 10 effector proteins were located in the intact T6SS cluster. Generally, if group 4 effector proteins were identified, PAAR was located in front of the tssA gene and contained only the PAAR_CT_2 domain. However, if the C-terminal toxin domain of PAAR was predicted, it was usually located downstream of tssM. The results predicted by Alphafold2 are shown in Supplementary Table 5, but the specific functions of most proteins in this study have yet to be elucidated.

A comparison of the background information of 241 strains of K. pneumoniae (Figure 5) revealed that all ST11 strains had two T6SS clusters and one of the clusters was intact. ST14 strains had three to four T6SS clusters, while all ST15 strains had three T6SS loci. However, due to the small number of ST14 and ST15 strains, the characteristics of the T6SS gene clusters require verification. For the distribution of effectors in isolates of different ST types, groups 1, 2, 4, 7, 8, and 9 effector proteins were identified in ST11 K. pneumoniae, and most ST11 isolates (60 of 75) had a group 4 effector protein with the DUF2235 domain. However, these effectors were not identified in ST14, ST15, and ST231 isolates. Group 5 effector proteins containing an alpha/beta hydrolases domain were predominantly distributed in ST14 and ST15 strains, while group 10 effector proteins (PAAR with a C-terminal toxin domain) were predominantly distributed in ST23 strains. In addition, all ST323 strains analyzed in this study contained group 6 proteins, which have predicted PLD activity. The corresponding relationship between ST type and effector is shown in Supplementary Table 2. Although the identified repertoire of potential effector proteins is diverse, almost all examined ST lineage strains shared identical predicted effectors regardless of the country of isolation or site of infection. To determine whether related strains of K. pneumoniae expressed similar groups of effectors, an ML phylogenetic tree of all strains included in these analyses was generated from a homologous single-copy gene (Figure 6). Functionally similar effectors were largely found in closely related strains, such as alpha/beta hydrolases (group 4). This result was consistent with the above findings.

Effector proteins are usually predicted in variable regions downstream of vgrG. Analysis of the T6SS cluster of K. pneumoniae revealed that four Tses were located downstream of vgrG, and the domains of these effectors included DUF2235 (pfam09994), bacterial toxin 44 (pfam15607), PLD superfamily (such as cl15239), and alpha/beta hydrolases (pfam06958). One of the effectors was PAAR, which is also a component of the T6SS. The VgrG corresponding to the first two effectors mentioned above belongs to the complete T6SS gene cluster. To delineate the relationship(s) among the VgrG proteins, a total of 201 VgrG amino acid sequences corresponding to the above four effector proteins were analyzed by ML phylogenetic tree (Figure 7). This phylogenetic analysis revealed that closely related VgrG proteins were always found upstream of similar effector-encoding genes. The VgrG of 47 strains could predict both DUF2235 and PAAR with the C-terminal toxic domain. Most VgrG proteins contained the same VgrG domains, including the vgr_GE superfamily (cl36942) and T6SS_Vgr (pfam13296) domains, while the C-terminal region of each VgrG protein contained an amino acid region with unknown function (DUF2345). However, most of the VgrG proteins (28 of 29) corresponding to bacterial toxin 44 do not have T6SS_Vgr domains.

The analysis of 241 K. pneumoniae genomes predicted that three strains did not encode any PAAR proteins in intact T6SS clusters. Another 12 strains with incomplete T6SS clusters were not included in the analysis; thus, a total of 226 strains were included in the analysis. The paar genes of most strains (158 of 226) were located in the variable region after vgrG. However, in other strains, the paar genes were behind tssA genes, and some (30 of 67) of them included a bacterial toxin 44 domain at the C-terminal. One strain, designated strain 47, had both kinds of PAAR-encoding genes in the intact T6SS cluster. Figure 8 shows the ML phylogenetic tree based on the amino acid sequences of the PAAR proteins. The PAAR proteins were separated into 10 groups, designated 1–10. The conservative domain of PAAR could be divided into two main types, namely, one type is PAAR with Ntox44 (pfam15607) and PAAR_CT_2 (cd14744) domains, while the other type only contains a PAAR_CT_2 domain. PAAR with bacterial toxin 44 at the C-terminus belongs to PAAR8 and PAAR9.

All the PAAR proteins in this analysis contained a PAAR-CT_2 domain and were members of the PAAR_like domain superfamily. This superfamily is subclassified into eight subgroups (cd14737–cd14744) based on the different predictive functions of the additional N- and C-terminal domains in the family members (Shneider et al., 2013). To further understand the characteristics of PAAR in K. pneumoniae, the MEME online program was utilized to analyze the motifs of different subgroups (as shown in Figure 9). Analysis of motif distributions showed that the length of the conserved motif varied from 7 to 17 amino acid residues, and all subgroups shared motif 2 (blue box). Compared with the other PAAR_like subgroups, the PAAR_CT_2 subgroup lacked motif 1 and motif 3, which might indicate that these subgroups have experienced divergence of gene functions.

A novel virulence factor PLD encoded in a T6SS locus was previously reported in K. pneumoniae Kp52.145 and may be involved in the pathogenicity of this bacterial species. PLD in K. pneumoniae belongs to the cardiolipin synthase subfamily and plays a role in balancing phosphatidylglycerol (PG) and cardiolipin (Lery et al., 2014). However, whether it is secreted by T6SS has not been verified. BLASTp revealed that proteins of the PLD superfamily were identified in no more than 16 strains in this study. As stated above (Table 1), this is a class of Tle effector proteins with lipase function. PLD belongs to the PLD superfamily and is widely distributed in prokaryotes and eukaryotes. Members of the PLD superfamily usually contain a conserved region of the catalytic HKD motif, H-x-K-(x)4-D, where x is any amino acid (Metrick et al., 2020). The sequence alignment results of the PLD family protein in K. pneumoniae (by Phyre2 fold recognition server) demonstrated 26% identity with known PLD alpha (PLDα; PDB:6KZ9) of Arabidopsis thaliana and 25% identity with phospholipase D2 (PLD2) catalytic domain (PDB:6OHM) of Homo sapiens. The structure and quality assessment parameters for the K. pneumoniae PLD protein are shown in Supplementary File 3. Structural deviations were assessed by RMSD values. The RMSD between PLD in K. pneumoniae and PLDα or PLD2 was 2.648 and 3.451, respectively. The predicted structure of PLD in K. pneumoniae showed a highly similar overlap with PLDα and PLD2 (Figure 10A), indicating a conserved organization of these enzymes. The PLDA (AAG06875.1) and PLDB (AAG08474.1) effectors of the T6SS in P. aeruginosa PAO1 also belong to the PLD family of proteins. Therefore, a multiple-amino-acid-alignment outcome profile for the PLD proteins of different species was generated (Supplementary File 2). The alignment outcome reveals significant differences among the five protein sequences, and the catalytic HKD motif was highly conserved (Figure 10B). While these PLD proteins have not been well conserved in the primary sequences, their domain structures show similar topologies. In addition, the conserved nature of the HKD catalytic motifs implies that the members of the PLD superfamily have similar reaction mechanisms against the substrate.