**Author details**

Fumiaki Uchiumi\* and Masashi Asai Department of Gene Regulation, Faculty of Pharmaceutical Sciences, Tokyo University of Science, Noda, Chiba-ken, Japan

\*Address all correspondence to<EMAIL_ADDRESS>
© 2018 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

**11**

*Introductory Chapter: Gene Expression Controlling System and Its Application to Medical Sciences*

[9] Gyoja F. Basic helix-loop-helix transcription factors in evolution: Roles in development of mesoderm and neural tissues. Genesis. 2017;**55**(9):e23051

1996;**262**:294-313

2003;**303**:11-34

[10] Wintjens R, Rooman M. Structural classification of HTH DNA-binding domains and protein-DNA interaction modes. Journal of Molecular Biology.

[11] Oikawa T, Yamada T. Molecular biology of the ETS family of transcription factors. Gene.

[12] Hsu T, Trojanowska M, Watson DK. ETS proteins in biological control

and cancer. Journal of Cellular Biochemistry. 2004;**91**(5):896-903

[13] Wei GH, Badis G, Berger MF, Kivioja T, Palin K, Enge M, et al. Genome-wide analysis of ETS-family DNA-binding in vitro and in vivo. The EMBO Journal. 2010;**29**:2147-2160

[14] Karim FD, Urness LD, Thummel CS, Klemsz MJ, McKercher SR, Celada A, et al. The ETS-domain: A new DNA-binding motif that recognizes a purine-rich core DNA sequence. Genes & Development. 1990;**4**:1451-1453

M. Homeodomain proteins: An update. Chromosoma. 2016;**125**(3):497-521

[17] Kamachi Y, Kondoh H. Sox proteins: Regulators of cell fate specification and differentiation. Development.

transcriptional activation of CCAAT promoters. Critical Reviews in

[15] Bürglin TR, Affolter

[16] Sturm RA, Herr W. The POU domain is a bipartite DNAbinding structure. Nature. 1988;**336**(6199):601-604

2013;**140**(20):4129-4144

[18] Dolfini D, Gatta R, Mantovani R. NF-Y and the

*DOI: http://dx.doi.org/10.5772/intechopen.80676*

[1] Carey MF, Peterson CL, Smale ST. A primer on transcriptional regulation in mammalian cells. In: Carey MF, Peterson CL, Smale ST, editors. Transcriptional Regulation in Eukaryotes: Concepts, Strategies, and Techniques. 2nd ed. Cold Spring Harbor, NY: Cold Spring Harbor

Laboratory Press; 2009. pp. 1-45

[2] Yang C, Bolotin E, Jiang T, Sladek FM, Martinez E. Prevalence of the initiator over the TATA box in human and yeast genes and identification of DNA motifs enriched in human TATA-less core promoters. Gene. 2007;**389**(1):52-65

[3] Hahn S, Donczew R. Mechanistic differences in transcription initiation at TATA-less and TATA-containing promoters. Molecular and Cellular Biology. 2018;**38**(1):e00448-e00417

[4] Yella VR, Bansal M. DNA structural features of eukaryotic TATA-containing and TATA-less promoters. FEBS Open

[5] Marbach-Bar N, Bahat A, Ashkenazi S, Golan-Mashiach M, Haimov O, Wu SY, et al. DTIE, a novel core promoter element that directs start site selection in TATA-less genes. Nucleic Acids Research. 2016;**44**(3):1080-1094

Bio. 2017;**7**:324-334

[6] Uchiumi F, Miyazaki S,

[7] Soutourina J. Transcription regulation by the mediator complex. Nature Reviews. Molecular Cell Biology.

[8] Carlberg C, Molnár F. Transcription factors. In: Carsten C, Molnar F, editors. Mechanisms of Gene

Regulation. 2nd ed. Dordrecht: Springer Science+Business Media; 2016. pp. 57-73

2018;**19**(4):262-274

Tanuma S. The possible functions of duplicated ETS (GGAA) motifs located near transcription start sites of various human genes. Cellular and Molecular Life Sciences. 2011;**68**(12):2039-2051

**References**

*Introductory Chapter: Gene Expression Controlling System and Its Application to Medical Sciences DOI: http://dx.doi.org/10.5772/intechopen.80676*


Increasing the driving increases coherence generation. On the open branches, this coherence results in rather significant frictional losses that quickly bring us to the turnover point where the friction dominates the cycle performance and the engine no longer produces work. This turnover point occurs in the Carnot cycle for much smaller values of the coherence than in the Otto cycle, presumably because, in the Carnot cycle, the open branches generate additional coherence. It is also the reason our graphs of the efficiency and the power for the Carnot cycle cannot reach lower cycle times, cf. Figures 10, 11, 15 and 16. The smooth behavior of both the efficiency and power as a function of the cycle time for local cycles gives way to oscillations at short times for global cycles (Figures 15 and 16). Coherence by nature oscillates, and these oscillations result in effectively constructive and destructive interference with the oscillation during the following stroke in our global cycles. Note that this feature also shows up for the global Otto, but at much faster cycle times. In general, the global Otto cycle is less sensitive to coherence (the coherence related operators and energy are on the same scale). While the sensitivity to coherence depends on the temperature gap between the system and bath during the open strokes; this lower sensitivity to coherence for the Otto cycle holds for any comparable gaps.

#### *9.3. Comparing to the Harmonic Working Fluid*

Engine models with the qubit and harmonic oscillator working medium have been the most popular quantum systems in the study of quantum heat devices [63,95,107,143–149]. These models share many common features, including the tradeoff between power and efficiency, and obtain the Carnot bound in the limit of large cycle time. Moreover, in the limit of low temperatures, the harmonic oscillator converges to the qubit model, and the thermodynamic performance should be equivalent. Despite the similarities, there are qualitative differences in the thermodynamic performance. The major differences between the two models can be traced to the dynamical algebra of the two, *SU*(2) and the Heisenberg-Weyl group *H*3. The former algebra is compact, while the latter is non-compact. A direct consequence is that the heat capacity of the harmonic oscillator increases with the temperature,

saturating for high temperatures. In contrast, the capacity of the qubit reaches a maximum value and then asymptotically decreases as *T*−2 in the high temperature regime.

The different algebra influences the dynamics as well. For example, for a constant adiabatic parameter protocol (non-adiabatic driving), the effective frequency of the qubit increases, while in the harmonic case the effective frequency decreases. In turn, the effective frequency determines the relaxation rate towards the instantaneous attractor. In both models, this rate increases monotonically with the effective frequency, and the relaxation rate will be influenced in an opposite manner. In addition, the detailed balance condition is also modified, which means the internal temperature of the qubit is reduced in the presence of non-adiabatic driving. When comparing the present global cycle to an analogous harmonic cycle [75], we find that the qubit cycle has a greater sensitivity to the presence of coherence (short cycle times). A possible explanation of this result is the shift to a lower internal temperature and higher relaxation rate, which destroys the coherence and nulls the extracted work.

In the operation of the engines, both working mediums allow for performing shortcut protocols [106] on the unitary and open strokes. For frictionless and shortcut cycles, the harmonic Otto cycle exhibits a maximum efficiency when optimizing the compression ratio, which corresponds to the classical endoreversible result *ηCA*, Equation (22) [95,107]. This result is independent of the power of the engine. In contrast, the qubit model reaches the Curzon–Ahlborn efficiency for an endoreversible cycle in the high temperature limit, Section 3.2 [64].

## *9.4. High Temperature Limit*

With the motto of learning from example, we can employ the qubit model to elucidate the path from the quantum first principle derivation to the classical FTT results. The key is the high temperature limit. This means that the polarization |*S*¯| is small and it can be used to expand the thermodynamical expressions to first order. In the elementary Carnot-type cycle this expansion leads to the Curzon–Ahlborn efficiency at maximum power *ηCA*, Equation (50), without relying on the linear Newtonian heat transfer law or the low dissipation limit [67,68,78]. If we consider the cost of driving (Section 5.4), we find that the entropy production rate Σ*u* at the high temperature limit can be cast into the template of the Onsager relations, Equation (71).

The qubit engine model operates in the low dissipation limit when the cycle period is very large. We observe the expected limit (insert of Figure 16), as the dissipated power scales inversely with the square of the cycle period P*diss* ∝ 1/*τ*2 *cyc*. This is generically true, as discussed in the next subsection.

## *9.5. Dissipation*

The problem of operating a heat engine while trying to minimize dissipation has a simple general answer: turn off the engine so nothing happens. This, of course, has a dissipation cost of zero; you cannot do better. In order to get an interesting answer to the minimum entropy production question, we have to require something to happen. For a heat engine, one natural choice is to carry entropy ΔS from the hot bath at temperature *Th* to the cold bath at temperature *Tc*. Once such a constraint is specified, the interesting optimizations for a finite-time heat engine range from minimum entropy production to maximum power, with maximum efficiency as merely an intermediate point [8].

The low dissipation limit has been recently employed to bound the dissipation in a heat engine while using the notion of thermodynamic distance [60,78,150]. This distance, defined on thermodynamic states using the second derivative of the entropy, bounds the finite-time cost of driving a system along a given trajectory in the linear response regime, i.e., the slow process limit. In fact the minimum cost of driving a system with relaxation time *τR* along a path of length L in time *τ* is *τR*L2/*τ*. Applied to our engine, when the qubit traces the cycle of length L, the dissipated power for slow processes must be at least *τR*L2/*τ*2 as matching our observations, cf. Figure 12. How to geometrically bound the dissipation in the context of our non-adiabatic master equation formalism for faster driving is not clear. A major difference is the coupling of energy and coherence. This issue is left for future efforts.

## *9.6. Experimental Connections*

Realization of engines with a working medium composed of an ensemble of spins is an expected development, such as in an NMR experiment [151]. The surprise is the ability to operate an engine with a single spin.

The experimental realization of single qubit engines and refrigerators is in the process of rapid development. In part, this progress is part of a larger effort in developing quantum technology. This breakthrough is due to the ability to cool the ambient environment to temperatures in the range, or colder than, the qubit energy gap [152]. Moreover, the rapid progress in manipulations, designed for quantum information processing, can be employed for the unitary strokes of quantum heat engines. Recently, experimental realizations of four-stroke cycles were demonstrated [42–44] as well as a two-stroke engine [153]. Another application for quantum engines is in quantum sensing and in particular, thermometry. Suggestions that are based on the transition point between a quantum engine cycle and a refrigerator have recently been proposed [154–157].

A more immediate goal is quantum refrigeration [152,158]. Because any quantum device operates in ultracold temperatures, the drive for miniaturization will require an on-the-chip quantum refrigerator replacing the cumbersome dilution refrigerators of today.

## **10. Conclusions**

The qubit thermal engine has been a source of insight concerning finite-time thermodynamics for 30 years, with its origins dating back to a time when qubits were still two-level-systems. Among the lessons from the model was the role of coherence in friction, like phenomena, in the unitary strokes. Further analysis revealed that the generation of coherence occurs on the unitary strokes and is separated from its dissipation, which occurs when the qubit is in contact with the thermal bath. This insight led first to the notion of quantum friction and later to the exploration of shortcuts to adiabaticity (STA). The analysis of the qubit engine generated a unifying overview of these finite-time thermodynamic phenomena. We have tried to present such an overview alongside our new findings.

Our study explores a new chapter in the behavior of the qubit: the study of a driven isothermal process. Recent progress in open system dynamics allowed for the treatment of thermalization processes, driven processes with time dependent Hamiltonians in contact with a heat bath [34,91,99]. This allowed for the Carnot cycle to be analyzed for shorter times than the previous linear response treatments.

This breakthrough was achieved by using a basis of eigenvectors of the instantaneous propagator of the system dynamics in the interaction representation while exploiting the dynamical *SU*(2) algebra of the qubit. The key to the breakthrough was the non-adiabatic master equation (NAME) [34] that could correctly describe the dynamics of thermalization in a thermodynamically consistent way and eliminate issues of time reordering [99]. In the qubit model, the instantaneous attractor is rotated from the energy direction by an amount that depends on the speed of the driving and it reveals the details of the coupling between energy and coherence.

In the present manuscript, we have analyzed local cycles where the coherence is required to vanish at the switching points between strokes as well as global cycles where coherence is set only by the driving protocols and it is carried from one stroke to the next. For local cycles we designed and implemented an STE protocol that mimics an isotherm. For the global cycles we used only constant adiabatic parameter *μ* trajectories. Both these types of cycles are analyzed and compared. As a result, for the qubit engine we are now able to assess the role of coherence along the isothermal strokes. The global cycles exhibit oscillations in their efficiency and power once we reach small cycle times (fast driving). These oscillations are due to the oscillations in coherence and how the timing of the switching between strokes happens to catch the coherence oscillation of the previous stroke. The Carnot-like cycles show enhanced sensitivity to coherence. This coherence can be reduced by weak measurement of the energy causing pure dephasing. Dephasing is damaging for short cycles periods, but it can be beneficial for intermediate cycle times. Another coherence related finding concerns the global Otto cycle in the sudden limit near which the coherence of the engine acts as a useful flywheel.

Our analysis calculates the entropy production rate from first principles. Using inertial coordinates for the qubit, we generally find that this entropy production naturally takes a flux-times-force form. In the high temperature limit, we find a linear relation between these fluxes and forces and an associated Onsager relation.

The present manuscript managed to compare the known behavior of local and global Otto and Carnot-like cycles. That is a lot to compare. We tried to focus on the new emergent phenomena in a coherent Carnot-type cycle. We found that, to discuss these phenomena, we needed a backdrop of related results to compare to. The above attempted synthesis is the outcome.

**Author Contributions:** This study is a chapter in a 25 year collaboration between P.S. and R.K., initiated by the study of an Otto spin engine with Tova Feldmann and Eitan Geva [63]. The focus on the Carnot cycle was inspired by a comment by Peter during a visit to Jerusalem, pointing out a discrepancy in the quantum heat engine studies between the many studies of Otto cycles and the few studies of Carnot and other cycles. The reason was the lack of an adequate thermodynamically consistent master equation for the isothermal strokes. It took four years to address the issue. The first step is the development of the non adiabatic master equation by R.D. and Amikam Levy [34]. R.D. continued the development setting it on firm theoretical grounds with the inertial theorem. Returning to the original objective, a shortcut to equilibration was developed [91]. The present study incorporates a broad perspective as well as many new results on the role of coherence in entropy production in the Carnot and Otto cycles. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by the Adams Fellowship Program of the Israel Academy of Sciences and Humanities and the Israel Science Foundation, grant number 2244/14.

**Acknowledgments:** We would like to thank our many contributors to this study: Tova Feldmann, Eitan Geva, Lajos Diosi, Jose P Palao, Yair Rezek, Karl Heinz Hoffmann, Amikam Levy, Raam Uzdin, James Nulton, Frank Boldt, Ander Tobalina, and Andrea Insinga.

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

## **Abbreviations**

The following abbreviations are used in this manuscript:

GKLS Gorini, Kossakowski, Linblad, Sudarshan; Master equation

CPTP Completely Positive Trace Preserving map

## **Appendix A. Representations of the Qubit State**

The qubit state can be described in many alternative ways. Each representation highlights a certain aspect of the engine. We will now summarize the different approaches and the relation between them. The basic construction relies on a set of orthogonal operators that form a closed Lie Algebra {*A*ˆ}

$$\mathbb{E}\left[\mathcal{A}\_{i\prime}\mathcal{A}\_{j}\right] = \sum\_{k} \mathbb{C}\_{k}^{ij} \mathcal{A}\_{k} \tag{A1}$$

where *Cij k* is the structure tensor of the algebra. The orthogonality relation:

$$\text{tr}\{\mathcal{A}\_i^\dagger \mathcal{A}\_j\} = \delta\_{ij} \quad , \tag{A2}$$

where the identity ˆ*I* is part of the set and all other operators are therefore traceless. Under these conditions the state *ρ*ˆ can be expanded as a linear combination of the set {*A*ˆ}

$$
\hat{\rho} = \frac{1}{N}\hat{I} + \sum\_{\hat{j}} \alpha\_{\hat{j}} \hat{A}\_{\hat{j}} \quad , \tag{A3}
$$

where *αj* = *A*ˆ*j* and *N* is the size of Hilbert space. An alternative formulation includes representing the state in terms of a generalized Gibbs state [101]:

$$\hat{\rho} = \frac{1}{Z} \exp\left(\sum\_{j} \lambda\_{j} \hat{A}\_{j}\right) \quad . \tag{A4}$$

The generalized Gibbs state is the maximum entropy state subject to the constrains of the expectation value *A*ˆ*j* = tr{*ρ*ˆ*A*ˆ*j*}, this leads to a set of non linear equations which determines the Lagrange multipliers *λj*. The forms (A3) and (A4) are unique once the expectation values *A*ˆ*j* are known.

It is convenient to express the generalized Gibbs state as a product form [159]

$$\hat{\rho} = \frac{1}{Z} \prod\_{k} \exp(\gamma\_{\hat{j}} \hat{A}\_{\hat{j}}) \ . \tag{A5}$$

This form is not unique since it depends on the order of operators. Once the order is set the coefficients {*γj*} are determined from the expectation values {*A*ˆ*j*}.

Specifically, for the qubit we employ the *SU*(2) algebra and three sets of orthogonal bases. The first basis set *s* = {*S*ˆ *x*, *S*ˆ *y*, *S*ˆ *z*}*T* represents the static polarization, where *S*ˆ *j* are the spin operators with the commutation relation of the *SU*(2) algebra [*S*ˆ *i*, *S*ˆ *j*] = *ih*¯ *ijkS*ˆ *k* and *σ*ˆ*j*. In terms of the Pauli operators *σ*ˆ*j* they are expressed as *S*ˆ *j* = *h*¯ 2 *σ*ˆ*j*. An arbitrary state is expressed as a linear combination of these operators in Equation (4). A geometric interpretation uses this set as a Cartesian basis in 3D Figure 1. A time-dependent rotation around the *S*ˆ *y* axis leads, up to a scaling, to the dynamical basis set *v* (*t*) = {*H*ˆ , *L*ˆ, *C*ˆ}*T*, Equation (32). The *SU*(2) algebra in the polarization basis defines a rotation in Liouville space

$$\mathcal{R}\_{\mathcal{Y}}(\boldsymbol{\phi}) = \exp\left(\frac{i}{\hbar}[\boldsymbol{\hat{S}}\_{\mathcal{Y}}, \boldsymbol{\bullet}]\boldsymbol{\phi}\right) = e^{\frac{i}{\hbar}\boldsymbol{\hat{S}}\_{\mathcal{Y}}\boldsymbol{\phi}} \bullet e^{-\frac{i}{\hbar}\boldsymbol{\hat{S}}\_{\mathcal{Y}}\boldsymbol{\phi}}\,,\tag{A6}$$

and the relation between the two basis sets

$$
\vec{v} = \Omega(t)\mathcal{R}\_{\mathcal{Y}}(\phi(t))\vec{s} \tag{A7}
$$

where *φ* = arccos(*ω*/Ω). The explicit dependence of Ω(*t*) and *φ*(*t*) means that the dynamical basis set *v* is time-dependent.

In terms of the dynamical basis set, the linear form of the state Equation (33) is defined. Since the Hamiltonian *H*ˆ is part of *v*, it is natural to define the generalized Gibbs state

$$\rho = \frac{1}{Z} \exp\left(-(\beta \hat{\mathcal{H}} + \lambda \hat{\mathcal{L}} + \gamma \hat{\mathcal{C}})\right) \tag{A8}$$

The standard Gibbs state is obtained when *γ* = *λ* = 0 then *β* = 2 *h*¯ Ω tanh−1 2 *H*ˆ *h*¯ Ω .

The third basis set is obtained from the eigenoperators of the free propagator {*χ*ˆ, *σ*ˆ, *σ*ˆ †}*T*, Equation (52). In terms of this set the linear form of the state in the interaction representation is given by Equation (58). In addition, the state can be expressed in a generalized Gibbs state form

$$\bar{\rho} = Z^{-1} \exp \left( -\left( \not\mathbb{R} \hat{\chi} + \not\mathbb{R} + \gamma^\* \mathcal{O}^\dagger \right) \right) \;. \tag{A9}$$

An equivalent product form has been previously been employed in Ref. [91]

$$
\overline{\rho} = Z^{-1} e^{\gamma \circ} e^{\overline{\beta} \circ} e^{\gamma \circ \circ^{\star}} \, \, \, \, \, \, \, \, \tag{A10}
$$

where *Z*˜ ≡ *Z*˜ (*t*) = tr(*ρ*˜*S* (*t*)) is the partition function, with time-dependent parameters *γ*˜ (*t*) and *β*˜.

In Section 5.4 we modify the basis set *g* = {*χ*ˆ, *σ*ˆ, *σ*ˆ †}*T* to a hemitian basis: *g* = {*χ*ˆ, *σ*ˆ*x*, *σ*ˆ*y*}*T* where *σ*ˆ*x* = √1 2 (*σ*ˆ + *σ*ˆ †) and *σ*ˆ*y* = √*i* 2 (*σ*ˆ − *σ*ˆ †). The new representation allows relating the basis *v* to *g* by a scaling and rotation around the *L*ˆ axis (cf. Figure 1)

$$
\vec{\mathbf{g}}' = \frac{\sqrt{2}}{\hbar \Omega} \mathcal{R}\_L(\vec{\xi}) \vec{\mathbf{v}}\ \text{ \(A11\)}
$$

where *ξ* = arccos(1/1 + *μ*2) and R*L*(*ξ*) = exp *i h*¯ Ω [*L*ˆ, •]*ξ* .

When studying the entropy production rate in Section 5.4 we take advantage of the form

$$\bar{\rho} = \bar{Z}^{-1} \exp \left( - \left( \bar{\beta} \hat{\chi} + \bar{\gamma}\_x \mathfrak{d}\_x + \bar{\gamma}\_y \mathfrak{d}\_y \right) \right) \tag{A12}$$

The {*β*¯, *γ*¯ *x*, *γ*¯ *y*} time-dependent parameters are defined by the eigenoperators expectation values. The relations are given by

$$
\langle \xi \rangle\_{\rm int} = f\left(r\right) \bar{\beta} \quad ; \quad \langle \mathfrak{d}\_{\rm x} \rangle\_{\rm int} = f\left(r\right) \bar{\gamma}\_{\rm x} \; ; \quad \langle \mathfrak{d}\_{\rm y} \rangle\_{\rm int} = f\left(r\right) \bar{\gamma}\_{\rm y} \; . \tag{A13}
$$

where • = tr(•*ρ*˜), *f* (*r*) = − √1 2*r* tanh √*r* 2 with *r* = *β*¯2 + *γ*¯ 2 *x* + *γ*¯ 2 *y*. In the large temperature limit *f* (*r*) ≈ −1 2 , which leads to simple relations between the thermodynamic fluxes and forces in the classical regime. Alternatively, the parameters can be expressed in terms of the expectation values

*β*¯ = *s*(*k*)*χ*ˆ *int* ; *γ*˜ *x* = *s*(*k*)*σ*ˆ*x int* ; *γ*˜ *y* = *s*(*k*)*σ*ˆ*y int* , (A14) where *s*(*k*) = log 1− √ √ 2*k* 2*k*+1 / √2*k* with *k* = *σ*ˆ*x* 2 *int* + *σ*ˆ*y* 2 *int* + *χ*ˆ 2 *int*

## **Appendix B. Explicit Expressions**

The transition matrix from the basis operators *v* to the basis operators *g*. The matrix appears in the inertial solution of the qubit, Equation (51):

$$P = \left[ \begin{array}{cccc} 1 & -\mu & -\mu \\ 0 & i\kappa & -i\kappa \\ \mu & 1 & 1 \end{array} \right] \tag{A15}$$

**Table A1.** *Cont.*

## **References**

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).


*Article* **Comparative Genomic Analysis of a Panton–Valentine Leukocidin-Positive ST22 Community-Acquired Methicillin-Resistant** *Staphylococcus aureus* **from Pakistan**

**Nimat Ullah 1,2 , Samavi Nasir 1 , Zaara Ishaq 1 , Farha Anwer 1 , Tanzeela Raza 1 , Moazur Rahman 3 , Abdulrahman Alshammari 4,\*, Metab Alharbi 4 , Taeok Bae 2 , Abdur Rahman 1 and Amjad Ali 1,\***

**Abstract:** *Staphylococcus aureus* (*S. aureus*) ST22 is considered a clinically important clone because an epidemic strain EMRSA-15 belongs to ST22, and several outbreaks of this clone have been documented worldwide. We performed genomic analysis of an *S. aureus* strain Lr2 ST22 from Pakistan and determined comparative analysis with other ST22 strains. The genomic data show that Lr2 belongs to *spa*-type t2986 and harbors staphylococcal cassette chromosome *mec* (SCC*mec*) type IVa(2B), one complete plasmid, and seven prophages or prophage-like elements. The strain harbors several prophage-associated virulence factors, including Panton–Valentine leukocidin (PVL) and toxic shock syndrome toxin (TSST). The single nucleotide polymorphism (SNPs)-based phylogenetic relationship inferred from whole genome and core genome revealed that strain Lr2 exhibits the nearest identities to a South African community-acquired methicillin-resistant *S. aureus* (CA-MRSA) ST22 strain and makes a separate clade with an Indian CA-MRSA ST22 strain. Although most ST22 strains carry *blaZ*, *mecA*, and mutations in *gyrA*, the Lr2 strain does not have the *blaZ* gene but, unlike other ST22 strains, carries the antibiotic resistance genes *erm(C)* and *aac(6*0 *)-Ie-aph(2*00*)-Ia*. Among ST22 strains analyzed, only the strain Lr2 possesses both PVL and TSST genes. The functional annotation of genes unique to Lr2 revealed that mobilome is the third-largest Cluster of Orthologous Genes (COGs) category, which encodes genes associated with prophages and transposons. This possibly makes methicillin-resistant *S. aureus* (MRSA) Lr2 ST22 strain highly virulent, and this study would improve the knowledge of MRSA ST22 strains in Pakistan. However, further studies are needed on a large collection of MRSA to comprehend the genomic epidemiology and evolution of this clone in Pakistan.

**Keywords:** comparative genomic analysis; CA-MRSA; EMRSA-15; PVL-positive; ST22

#### **1. Introduction**

*Staphylococcus aureus*(*S. aureus*) has long been considered an important human pathogen that causes both hospital- and community-acquired (HA and CA) infections. This pathogen is also well-known for acquiring resistance to a variety of antibiotics [1,2]. Particularly, methicillin resistance is becoming more common in *S. aureus*, posing a growing public health threat with substantial mortality and morbidity. The methicillin-resistant *S. aureus*

**Citation:** Ullah, N.; Nasir, S.; Ishaq, Z.; Anwer, F.; Raza, T.; Rahman, M.; Alshammari, A.; Alharbi, M.; Bae, T.; Rahman, A.; et al. Comparative Genomic Analysis of a Panton–Valentine Leukocidin-Positive ST22 Community-Acquired Methicillin-Resistant *Staphylococcus aureus* from Pakistan. *Antibiotics* **2022**, *11*, 496. https://doi.org/10.3390/ antibiotics11040496

Academic Editor: Teresa V. Nogueira

Received: 2 March 2022 Accepted: 6 April 2022 Published: 8 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

(MRSA) produces a low-affinity penicillin-binding protein (PBP2a), which enables it to confer resistance to almost all beta-lactam antibiotics [3,4]. The PBP2a is encoded by the gene *mecA*, which is carried by a large mobile genetic element known as staphylococcal cassette chromosome *mec* (SCC*mec*) [5]. In the past, MRSA caused mostly hospital-acquired infections; however, now, it also causes community-acquired infections, including necrotizing pneumonia and skin and soft tissue infection (SSTI) [6]. Apart from humans, MRSA can also colonize other animal species, particularly livestock, and CA-MRSA infections can also be caused by livestock-associated MRSA (LA-MRSA) [7].

Compared with HA-MRSA, the CA-MRSA strains belong to a diverse lineage with smaller SCC*mec*, for example, SCC*mec* IV, and harbors distinct virulence factors, particularly Panton–Valentine leukocidin (PVL) [8]. PVL is a two-component toxin encoded by prophage-associated genes *lukF-PV* and *lukS-PV* [9,10]. The toxin targets phagocytic leucocytes and triggers leukocyte lysis and/or apoptosis by forming pores. Therefore, PVLpositive CA-MRSA strains cause various infections, particularly necrotizing pneumonia and SSTI [11]. Although not all CA-MRSA contain the *pvl* gene, it can be considered as a molecular marker for CA-MRSA along with the SCC*mec* type IV because PVL is associated with the global spread of CA-MRSA [12]. The well-known clones of PVL-positive CA-MRSA include ST8, frequently found in the United States [13]; ST80, a Europe clone [10]; ST59, reported in Taiwan [14]; ST30, a global strain [15]; and ST22 (SCC*mec* IV/PVL-positive), which emerged in the United Kingdom and is now spreading worldwide [16]. *S. aureus* ST22 is considered a clinically important clone because an epidemic strain EMRSA-15 belongs to ST22. The EMRSA-15 emerged in the United Kingdom in the 1990s [17] and very quickly expanded all over Europe, Australia, and Asia [18]. EMRSA-15 has SCC*mec* type IV, a deletion of 2268 bp region in fibronectin-binding protein (FnBP) locus, and a point mutation in the *ureC* gene, and it is resistant to fluoroquinolone and macrolide. The genomic analysis of ST22 strains revealed that the epidemic clone EMRSA-15 makes a distinct clade (ST22-A clade) [18].

Although there is a wealth of literature available on comprehensive epidemiological and molecular characterization of MRSA in the USA, Europe, and Asia, there is a scarcity of detailed genomic characterization of MRSA strains in Pakistan. According to some recent studies, Pakistan has a greater rate of MRSA infection, and the data suggest that MRSA clones are becoming more diverse [19–21]. The whole-genome analysis provides high resolution in both global and local outbreak investigations as well as further exploring genes associated with pathogenicity and antibiotic resistance [22]. Therefore, whole-genomebased approaches have become an indispensable tool for the real-time comparative genomic study of a variety of pathogens in terms of antibiotic resistance, the emergence of new virulent clones, and niche adaptation [23]. Furthermore, investigating virulence- and antibiotic resistance-associated genes can help identify risk factors for MRSA infection and develop efficient infection control programs. Given the lack of genome-based surveillance and epidemiological studies of MRSA in Pakistan, we aimed to study a CA-MRSA ST22 strain from Pakistan at the genomic level and compare the virulence and antibiotic resistance genes and evolutionary relationships to other ST22 strains. To our knowledge, this is the first genome-based analysis of PVL-positive CA-MRSA ST22 from Pakistan, and we found many genomic differences compared to other complete genomes of MRSA ST22 strains.

#### **2. Results**

#### *2.1. Preliminary Identification and Phenotypic Antibiotic Resistance Pattern*

The isolate Lr2 was found positive for coagulase and catalase while negative for oxidase. The isolate was resistant to oxacillin, methicillin, ampicillin, erythromycin, gentamicin, streptomycin, clindamycin, ciprofloxacin, linezolid, and tetracycline but susceptible to vancomycin, chloramphenicol, and rifampicin (Table 1). The amplification of the 500 bp *mecA* gene confirmed the nature of methicillin resistance.

**Table 1.** The phenotypic and genotypic resistance profile of Lr2.

R—resistant; S—susceptible; ND—not determined.

#### *2.2. Genomic Features and Epidemiological Typing*

The genomic DNA (gDNA) was successfully sequenced, and a total of 1,363,009 reads with a mean length of 565 bp (272-fold sequence coverage) were obtained. The kmerFinder 2.0 identified the isolate as *S. aureus*. The N50, N75, and L50 values are 131,746, 94,139, and 7, respectively; de novo assembly resulted in 52 contigs (>500 bp), and the size of the largest contig is 425,597 bp. The genome size of Lr2 is 2,831,239 bp with 32.7% GC content. Genome annotation predicted 2835 total genes, of which 2768 are CDSs. The number of predicted tRNA are 56 and rRNA are 7 (5S = 4, 16S = 2, 23S = 1). In silico genome-based epidemiological typing revealed that strain Lr2 belongs to ST22, *spa*-type t2986, and carries SCC*mec* type IVa(2B) (Table 2).

**Table 2.** The genomic features and characteristics of the MRSA strain Lr2 ST22.

#### *2.3. Predicted Antibiotic Resistance Determinants and Virulence Factors*

The Lr2 genome contains a methicillin resistance gene *mecA*, an aminoglycoside resistance gene *aac*(*6* 0 )*-Ie-aph*(*2* 00)*-Ia*, a point mutation in *gyrA* conferring resistance to fluo-

*Antibiotics* **2022**, *11*, x FOR PEER REVIEW 4 of 15

*2.3. Predicted Antibiotic Resistance Determinants and Virulence Factors*

roquinolones, and *erm(C)*—a plasmid-associated gene conferring resistance to lincosamide, macrolide (erythromycin), and streptogramin. However, the *blaZ* gene conferring penicillin resistance was found absent in Lr2 (Table 1 and Figure 1). mide, macrolide (erythromycin), and streptogramin. However, the *blaZ* gene conferring penicillin resistance was found absent in Lr2 (Table 1 and Figure 1).

The Lr2 genome contains a methicillin resistance gene *mecA*, an aminoglycoside resistance gene *aac*(*6′*)*-Ie-aph*(*2″*)*-Ia*, a point mutation in *gyrA* conferring resistance to fluoroquinolones, and *erm(C)*—a plasmid-associated gene conferring resistance to lincosa-

**Figure 1.** Circular visualization of Lr2 genome via CGViewer, showing virulence factors (light green) and antibiotic resistance determinants (light blue). **Figure 1.** Circular visualization of Lr2 genome via CGViewer, showing virulence factors (light green) and antibiotic resistance determinants (light blue).

The Virulence Factor Database (VFDB) predicated 62 virulence factors responsible for adherence, protein secretion, immune evasion, and toxin production. Adherence-associated genes are autolysin (*atl*); clumping factor A (*clfA*) and B (*clfB*); an elastin binding protein (*ebp*); cell-wall-associated fibronectin-binding protein (*ebh*); collagen adhesion (*cna*); fibrinogen-binding protein (*efb*); fibronectin-binding protein A (*fnbA*) and B (*fnbB*); Ser-Asp rich fibrinogen-binding protein C (*sdrC*), D (*sdrD*), and E (*sdrE*); intercellular adhesin (*icaA*, *icaB*, and *icaC*); and staphylococcal protein A (*spa*) (Figure 1). The virulome also consists of several exoenzymes genes, including hyaluronate lyase (*hysA*), staphylokinase (*sak*), cysteine protease (*sspB* and *sspC*), serine V8 protease (*sspa*), lipase (*geh* and *lip*), staphylocoagulase (*coa*), and thermonuclease (*nuc*). However, a six gene-cluster (*splA, splB, splC, splD, splE,* and *splF*) of serine was absent in Lr2 (Figure 1). The Lr2 genome also contains virulence factors associated with the host immune evasion, such as staphylococ-The Virulence Factor Database (VFDB) predicated 62 virulence factors responsible for adherence, protein secretion, immune evasion, and toxin production. Adherenceassociated genes are autolysin (*atl*); clumping factor A (*clfA*) and B (*clfB*); an elastin binding protein (*ebp*); cell-wall-associated fibronectin-binding protein (*ebh*); collagen adhesion (*cna*); fibrinogen-binding protein (*efb*); fibronectin-binding protein A (*fnbA*) and B (*fnbB*); Ser-Asp rich fibrinogen-binding protein C (*sdrC*), D (*sdrD*), and E (*sdrE*); intercellular adhesin (*icaA*, *icaB*, and *icaC*); and staphylococcal protein A (*spa*) (Figure 1). The virulome also consists of several exoenzymes genes, including hyaluronate lyase (*hysA*), staphylokinase (*sak*), cysteine protease (*sspB* and *sspC*), serine V8 protease (*sspa*), lipase (*geh* and *lip*), staphylocoagulase (*coa*), and thermonuclease (*nuc*). However, a six gene-cluster (*splA*, *splB*, *splC*, *splD*, *splE*, and *splF*) of serine was absent in Lr2 (Figure 1). The Lr2 genome also contains virulence factors associated with the host immune evasion, such as staphylococcal complement inhibitor (*scn*), IgG-binding protein (*sbi*), and chemotaxis inhibiting protein (*chp*) (Figure 1). It also carries type VII secretion system associated virulence factors, including soluble cytosolic (*esaB* and *esaG*) and membrane-associated protein A (*essA*), B

cal complement inhibitor (*scn*), IgG-binding protein (*sbi*), and chemotaxis inhibiting protein (*chp*) (Figure 1). It also carries type VII secretion system associated virulence factors, (*essB*), and C (*essC*), and A (*esxA*) but no *esxB*, *esxC*, *esxD*, *esaD*, and *esaE* were found in Lr2 genome (Figure 1).

The genome of Lr2 carries toxin associated genes, particularly alpha-hemolysin gene (*hla*), gamma hemolysin (*hlgA*, *hlgB*, and *hlgC*), delta hemolysin gene (*hld*), exotoxins (*set10*, *set13*, *set15*, and *set16*), enterotoxin-like O (*selo*), exotoxins (*set1*, *set2*, *set4*, *set7*, *set10*, *set13*, *set15*, and *set16*), enterotoxin-like M (*selm*), enterotoxin-like N (*seln*), enterotoxin B (*seg*), PVL (*lukF-PV* and *lukS-PV*), and TSST (tsst-1). However, beta-hemolysin gene (*hlb*), exfoliative toxin type A (*eta*), enterotoxin A (*sea*), enterotoxin Yent1 (*yent1*), enterotoxin-like K (*selk*), exotoxins (*set34*, *set37*, and *set39*), and leukotoxins (*lukM*, *lukD*, and *lukE*) were found absent (Figure 1).

#### *2.4. Predicted Plasmids and Prophages*

PlasmidFinder identified one plasmid of 2402 bp length with 99.75% sequence similarity to *S. aureus* strain E14 plasmid pDLK1 (GU562624.1). The pDLK1 plasmid consists of an erythromycin resistance gene (*emrC*) and a replication gene (*repL*), and it encodes no other factors.

The genome of Lr2 strain presents seven putative prophages including a complete prophage (PHAGE\_Staphy\_tp310\_3), five incomplete prophages (PHAGE\_Bacill\_IEBH, PHAGE\_Staphy\_PT1028, PHAGE\_Clostr\_phiC2, PHAGE\_Staphy\_phiPVL\_CN125, and PHAGE\_Staphy\_96), and a questionable prophage (PHAGE\_Staphy\_phiN315) (Table 3). Several virulence factors were predicted in the identified prophages, e.g., *sec*, *sell*, and *tsst* in the prophage PHAGE\_Bacill\_IEBH; the *ebp*, *lukF-PV*, and *lukS-PV* in PHAGE\_Staphy\_phiPVL\_CN125; the *seg*, *sei*, *yent2*, *selm*, *seln*, and *selo* in PHAGE\_Staphy\_96; *sak*, *chp*, and *scn* in PHAGE\_Staphy\_tp310\_3; and *cna* in PHAGE\_Staphy\_phiN315 (Table 3). However, no antibiotic resistance gene was found in any of the identified prophages.

#### **Table 3.** Characteristics of prophages present in MRSA Lr2 ST22.

#### *2.5. Comparative Phylogenetic Analysis*

The whole-genome-based single nucleotide polymorphisms (SNPs) phylogenetic analysis of the Lr2 ST22 strain with all other ST22 strains indicated that the Lr2 strain exhibits the nearest identities to South African CA-MRSA ST22 strain 71A\_S11 (CP010940) and Indian strain CA-MRSA ST22 VB31683 (CP035671) (Figure 2).

**Figure 2.** Whole-genome-based SNPs phylogenetic tree of CA-MRSA ST22 strain Lr2 and other ST22 strains. The heatmap shows presence and absence of antimicrobial resistance determinants and virulence factors. **Figure 2.** Whole-genome-based SNPs phylogenetic tree of CA-MRSA ST22 strain Lr2 and other ST22 strains. The heatmap shows presence and absence of antimicrobial resistance determinants and virulence factors.

#### *2.6. Comparative Analysis of Antibiotic Resistance Determinants and Virulence Factors 2.6. Comparative Analysis of Antibiotic Resistance Determinants and Virulence Factors*

The heatmap shows that aminoglycoside resistance gene *aac*(*6′*)*-Ie-aph*(*2″*)*-Ia* and a plasmid-associated gene *erm(C)* conferring resistance to lincosamide, macrolide (erythromycin), and streptogramin are only present in Lr2 and Indian strain VB31683. Both Lr2 and Indian strain VB31683 also have a point mutation in *gyrA,* conferring fluoroquinolones resistance, which is absent in a majority of ST22 strains. However, the penicillinresistance gene (*blaZ*) is absent in Lr2, while it is present in almost all ST22 strains. In contrast, the South African strain 71A\_S11, a sister strain to Lr2, carries only *blaZ* and *mecA* genes (Figure 2). However, South African strain 71A\_S11 has a very similar profile of virulence factors to that of the Lr2 strain since both TSST (tsst) and PVL (*lukF-PV* and *lukS-PV*) genes are only present in these two ST22 strains. In addition, Lr2 strains also harbor clumping factor B (*clfB*), which is absent in South African strain 71A\_S11 and all other ST22 strains. Interestingly, leukotoxin D (*lukD*) and E (*lukE*) were absent in all ST22 strains, including the Lr2 strain (Figure 2). The heatmap shows that aminoglycoside resistance gene *aac*(*6* 0 )*-Ie-aph*(*2* 00)*-Ia* and a plasmid-associated gene *erm(C)* conferring resistance to lincosamide, macrolide (erythromycin), and streptogramin are only present in Lr2 and Indian strain VB31683. Both Lr2 and Indian strain VB31683 also have a point mutation in *gyrA*, conferring fluoroquinolones resistance, which is absent in a majority of ST22 strains. However, the penicillin-resistance gene (*blaZ*) is absent in Lr2, while it is present in almost all ST22 strains. In contrast, the South African strain 71A\_S11, a sister strain to Lr2, carries only *blaZ* and *mecA* genes (Figure 2). However, South African strain 71A\_S11 has a very similar profile of virulence factors to that of the Lr2 strain since both TSST (tsst) and PVL (*lukF-PV* and *lukS-PV*) genes are only present in these two ST22 strains. In addition, Lr2 strains also harbor clumping factor B (*clfB*), which is absent in South African strain 71A\_S11 and all other ST22 strains. Interestingly, leukotoxin D (*lukD*) and E (*lukE*) were absent in all ST22 strains, including the Lr2 strain (Figure 2).

#### *2.7. Pan-Genome Analysis and Functional Annotation 2.7. Pan-Genome Analysis and Functional Annotation*

The pan genome of MRSA ST22 strains consists of 2941 genes, of which 2408 (81.8%) genes are part of the core genome, 321 (10.9%) genes are accessory, and 212 (7.2%) genes are unique (Figure 3A). The plot of the pan and core genome shows fluctuation in the number of pan and core genes, which suggests that pan- and core-genome sizes are not stable (Figure 3B). The pan genome of MRSA ST22 strains consists of 2941 genes, of which 2408 (81.8%) genes are part of the core genome, 321 (10.9%) genes are accessory, and 212 (7.2%) genes are unique (Figure 3A). The plot of the pan and core genome shows fluctuation in the number of pan and core genes, which suggests that pan- and core-genome sizes are not stable (Figure 3B).

*Antibiotics* **2022**, *11*, x FOR PEER REVIEW 7 of 15

**Figure 3.** The pan genome of MRSA ST22 strains. (**A**) The pie chart shows number of genes in core, accessory, and unique genomes of twenty genomes of ST22. (**B**) The pan genome vs. core genome plot of ST22 strains. **Figure 3.** The pan genome of MRSA ST22 strains. (**A**) The pie chart shows number of genes in core, accessory, and unique genomes of twenty genomes of ST22. (**B**) The pan genome vs. core genome plot of ST22 strains. accessory, and unique genomes of twenty genomes of ST22. (**B**) The pan genome vs. core genome plot of ST22 strains. The functional annotation of core genes shows that 960 (44.7%) Cluster of Orthologous

The functional annotation of core genes shows that 960 (44.7%) Cluster of Orthologous Genes (COGs) are involved in metabolism and transport; 440 (20.5%) in information, storage, and processing; and 416 (19.4%) in cellular processing and signaling; 317 (14.8%) are poorly characterized, and 14 (0.7%) are found to be associated with mobilome (Figure 4A). The functional annotation of the genes unique to MRSA strain Lr2 ST22 revealed that 18 (32.1%) COGs are involved in metabolism and transport, 18 (32.1%) in information, storage, and processing, 9 (16.1%) are associated with mobilome, 6 (10.7%) are poorly characterized (either involve in general function and/or unknown function), and 5 (8.9%) are involved in cellular processing and signaling (Figure 4B). The functional annotation of core genes shows that 960 (44.7%) Cluster of Orthologous Genes (COGs) are involved in metabolism and transport; 440 (20.5%) in information, storage, and processing; and 416 (19.4%) in cellular processing and signaling; 317 (14.8%) are poorly characterized, and 14 (0.7%) are found to be associated with mobilome (Figure 4A). The functional annotation of the genes unique to MRSA strain Lr2 ST22 revealed that 18 (32.1%) COGs are involved in metabolism and transport, 18 (32.1%) in information, storage, and processing, 9 (16.1%) are associated with mobilome, 6 (10.7%) are poorly characterized (either involve in general function and/or unknown function), and 5 (8.9%) are involved in cellular processing and signaling (Figure 4B). Genes (COGs) are involved in metabolism and transport; 440 (20.5%) in information, storage, and processing; and 416 (19.4%) in cellular processing and signaling; 317 (14.8%) are poorly characterized, and 14 (0.7%) are found to be associated with mobilome (Figure 4A). The functional annotation of the genes unique to MRSA strain Lr2 ST22 revealed that 18 (32.1%) COGs are involved in metabolism and transport, 18 (32.1%) in information, storage, and processing, 9 (16.1%) are associated with mobilome, 6 (10.7%) are poorly characterized (either involve in general function and/or unknown function), and 5 (8.9%) are involved in cellular processing and signaling (Figure 4B).

**Figure 4.** Functional annotation of core genome and genes unique to strain Lr2. (**A**) Distribution of COGs in core genome of MRSA strains ST22. (**B**) Distribution of COGs unique to strain Lr2. **Figure 4.** Functional annotation of core genome and genes unique to strain Lr2. (**A**) Distribution of COGs in core genome of MRSA strains ST22. (**B**) Distribution of COGs unique to strain Lr2. **Figure 4.** Functional annotation of core genome and genes unique to strain Lr2. (**A**) Distribution of COGs in core genome of MRSA strains ST22. (**B**) Distribution of COGs unique to strain Lr2.

The largest core-genome category of MRSA strains ST22 consists of genes with functions associated with metabolism and transport, which are further categorized as C (110 genes), involved in energy conversion and production; E (223 genes), involved in transport and metabolism of amino acid; F (75 genes), involved in transport and metabolism of nucleic acids; G (179 genes), involved in transport and metabolism of carbohydrates; H (121 genes), involved in transport and metabolism of coenzyme; I (80 genes), involved in lipid transport and metabolism; P (140 genes), involved in inorganic ions transport and metabolism; and Q (32 genes), involved in biosynthesis, transport, and catabolism of secondary metabolites (Figure 5). The second-largest COGs category with functions related to information, storage, and processing are individually categorized as A (2 genes), involved in modification and processing of RNAs; B (1 gene), involved in structure and dynamics of chromatins; J (200 genes), involved in ribosomal structure, biogenesis, and translation; K (139 genes), involved in transcription; and L (98 genes), involved in repair, recombination, and replication (Figure 5). The COGs category of core genome containing the lowest number of genes belonging to mobilome, prophages, and transposons are categorized as X, which accounts for about 0.7% of the total core genome (Figure 5). tions associated with metabolism and transport, which are further categorized as C (110 genes), involved in energy conversion and production; E (223 genes), involved in transport and metabolism of amino acid; F (75 genes), involved in transport and metabolism of nucleic acids; G (179 genes), involved in transport and metabolism of carbohydrates; H (121 genes), involved in transport and metabolism of coenzyme; I (80 genes), involved in lipid transport and metabolism; P (140 genes), involved in inorganic ions transport and metabolism; and Q (32 genes), involved in biosynthesis, transport, and catabolism of secondary metabolites (Figure 5). The second-largest COGs category with functions related to information, storage, and processing are individually categorized as A (2 genes), involved in modification and processing of RNAs; B (1 gene), involved in structure and dynamics of chromatins; J (200 genes), involved in ribosomal structure, biogenesis, and translation; K (139 genes), involved in transcription; and L (98 genes), involved in repair, recombination, and replication (Figure 5). The COGs category of core genome containing the lowest number of genes belonging to mobilome, prophages, and transposons are categorized as X, which accounts for about 0.7% of the total core genome (Figure 5).

The largest core-genome category of MRSA strains ST22 consists of genes with func-

*Antibiotics* **2022**, *11*, x FOR PEER REVIEW 8 of 15

**Figure 5.** Comparative functional categories of COG in core genome of ST22 strains and COGs unique to strain Lr2. **Figure 5.** Comparative functional categories of COG in core genome of ST22 strains and COGs unique to strain Lr2.

As expected, the largest category of COGs unique to Lr2 involved in metabolism and transport and are individually categorized as six genes in the C category (involved in energy conversion and production), four in the E category (involved in transport and metabolism of amino acid), and eight in the G category (involved in transport and metabolism of carbohydrates). No protein of F (involved in transport and metabolism of nucleic acids), H (involved in transport and metabolism of coenzyme), I (involved in lipid transport and metabolism), P (involved in inorganic ions transport and metabolism), and Q (involved in biosynthesis, transport, and catabolism of secondary metabolites) category was found (Figure 5). Interestingly, the third-largest category of COGs unique to Lr2 belongs to the X category, which encodes genes associated with prophages and transposons As expected, the largest category of COGs unique to Lr2 involved in metabolism and transport and are individually categorized as six genes in the C category (involved in energy conversion and production), four in the E category (involved in transport and metabolism of amino acid), and eight in the G category (involved in transport and metabolism of carbohydrates). No protein of F (involved in transport and metabolism of nucleic acids), H (involved in transport and metabolism of coenzyme), I (involved in lipid transport and metabolism), P (involved in inorganic ions transport and metabolism), and Q (involved in biosynthesis, transport, and catabolism of secondary metabolites) category was found (Figure 5). Interestingly, the third-largest category of COGs unique to Lr2 belongs to the X category, which encodes genes associated with prophages and transposons (Figure 5). Cellular processing and signaling is the lowest COGs category with five genes, which are individually classified as three genes in the U category (involved in intracellular traf-

ficking, secretions, and vesicular transport) and two in the T category (involved in signal transduction mechanisms) (Figure 5). nal transduction mechanisms) (Figure 5). In the pan genome of MRSA ST22 strains, 5 COGs were found to be present in all 19

(Figure 5). Cellular processing and signaling is the lowest COGs category with five genes, which are individually classified as three genes in the U category (involved in intracellular trafficking, secretions, and vesicular transport) and two in the T category (involved in sig-

In the pan genome of MRSA ST22 strains, 5 COGs were found to be present in all 19 strains while absent in the Lr2 strain. Functional annotation showed that those COGs are involved in transposase (COG3666) and putative lipase (COG5153), which is essential for the disintegration of autophagic bodies inside the vacuole. strains while absent in the Lr2 strain. Functional annotation showed that those COGs are involved in transposase (COG3666) and putative lipase (COG5153), which is essential for the disintegration of autophagic bodies inside the vacuole. The core-genome phylogenetic analysis grouped the Lr2 strain with a South African

The core-genome phylogenetic analysis grouped the Lr2 strain with a South African ST22 strain 71A S11 strain and an Indian ST22 strain VB31683 (Figure 6). The phylogenetic relationships determined from the core genome are in agreement with the phylogenetic tree from whole-genome-based SNPs (Figure 2). The strains share more than 80% of genes; however, the accessory genome of MRSA strain Lr2 contains many genes that are absent in most of the strains (Figure 6). ST22 strain 71A S11 strain and an Indian ST22 strain VB31683 (Figure 6). The phylogenetic relationships determined from the core genome are in agreement with the phylogenetic tree from whole-genome-based SNPs (Figure 2). The strains share more than 80% of genes; however, the accessory genome of MRSA strain Lr2 contains many genes that are absent in most of the strains (Figure 6).

*Antibiotics* **2022**, *11*, x FOR PEER REVIEW 9 of 15

**Figure 6.** Core-genome SNPs-based tree and heatmap of presence and absence of accessory genes in 20 strains of MRSA ST22. **Figure 6.** Core-genome SNPs-based tree and heatmap of presence and absence of accessory genes in 20 strains of MRSA ST22.

#### **3. Discussion**

**3. Discussion** The MRSA SCC*mec* type IV ST22 is a clinically important clone (EMRSA-15), and several outbreaks of ST22 CA-MRSA and/or CA-MSSA have been documented worldwide [24–26]. The ST22 clone is considered a genotype of HA-EMRSA-15; however, the presence of ST22 clones in the general population and its correlation with community settings have been reported recently [27]. The dissemination of MRSA from hospitals to the community and vice versa and the emergence of strains resistant to β-lactams is a major cause of concern worldwide. Consequently, tracking of emerging clones of MRSA at the genomic level is required to prevent further spread and to guide the development of rapid diagnostic tools and therapy. Although a small number of studies from Pakistan molecularly characterized MRSA, including SCC*mec*, PFGE, and MLST typing, as well as some studies that investigated PVL genes [19,28–30], none of them performed a genome-based analysis. This shows that there are limited epidemiological and genomic data on MRSA The MRSA SCC*mec* type IV ST22 is a clinically important clone (EMRSA-15), and several outbreaks of ST22 CA-MRSA and/or CA-MSSA have been documented worldwide [24–26]. The ST22 clone is considered a genotype of HA-EMRSA-15; however, the presence of ST22 clones in the general population and its correlation with community settings have been reported recently [27]. The dissemination of MRSA from hospitals to the community and vice versa and the emergence of strains resistant to β-lactams is a major cause of concern worldwide. Consequently, tracking of emerging clones of MRSA at the genomic level is required to prevent further spread and to guide the development of rapid diagnostic tools and therapy. Although a small number of studies from Pakistan molecularly characterized MRSA, including SCC*mec*, PFGE, and MLST typing, as well as some studies that investigated PVL genes [19,28–30], none of them performed a genomebased analysis. This shows that there are limited epidemiological and genomic data on MRSA reported from Pakistan, leaving a substantial knowledge gap in our understanding of this important human pathogen that what are the dominant clones of MRSA circulating

in Pakistan and what makes them antibiotic-resistant and virulent. Therefore, we performed an in-depth genome-based analysis of a CA-MRSA ST22 strain from Pakistan and its comparative genomic analysis with other MRSA ST22 available in the NCBI database.

Genome analysis indicated that Lr2 belongs to MLST type ST22, spa-type t2986, and SCC*mec* type IVa(2B). According to previous studies, *SCCmec* types IV is associated with CA-MRSA, whereas HA-MRSA mostly exhibits *SCCmec* types I, II, or III [22,31]. The strain was also found to be positive for PVL (*lukF-PV* and *lukS-PV*), which is commonly considered a CA-MRSA marker [22,32]. Another study reported a high rate of MRSA PVLpositive isolates from Pakistan, which were classified as community-associated [28]. The PVL-positive strains are responsible for abscesses formation, tissue necrosis, and increased inflammatory responses [33]. In addition to PVL, TSST-producing MRSA strains tend to cause more complex infections. The presence of both *tsst* and PVL (*lukF-PV* and *lukS-PV*) genes poses a concern for increased virulence of Lr2. A recent study also reported both PVL (*lukF-PV* and *lukS-PV*) genes and the *tsst* gene in MRSA isolates from this region [19]; however, the presence of PVL genes in combination with *tsst* gene appears to be extremely rare. The *tsst* and PVL (*lukF-PV* and *lukS-PV*) genes were found on PHAGE\_Staphy\_PT1028 and PHAGE\_Staphy\_phiPVL\_CN125, respectively. The genes encoding staphylococcal enterotoxins (i.e., *selm*, *seg*, *yent2*, *sei*, *selo*, and *seln*) were found on PHAGE\_Staphy\_96, which are involved in staphylococcal food poisoning and belong to the enterotoxin gene cluster (*egc*) [34]. In addition, the staphylococcal complement inhibitor (SCIN) encoding gene *scn* as well as chemotaxis inhibiting protein of staphylococcus (CHIPS) encoding gene *chp* were also found on prophages of the sequenced strain. These proteins have a significant role in the host-pathogen interaction and help the pathogen to evade the host's immune response [35]. Therefore, the prophages predicted in this strain can act as a reservoir of virulence factors and could contribute to strain evolution towards high virulency and pose a serious threat [36]. Aside from the prophage-associated virulence factors, genes involved in the *S. aureus* type VII secretion system, such as *essA*, *essB*, *essC, esxA, esaA*, *esaB*, and *esaG*, were also present in PVL-positive CA-MRSA strain Lr2 (Figure 1), which promote bacterial persistence [37]. The *esxA* gene is involved in the colonization and dissemination of *S. aureus* as well as triggering T-cell immune response [38]. The genes encoding adhesion factors, including fibrinogen-binding protein (*efb*), collagen-binding protein (*cna*), Ser-Asp rich fibrinogen-binding proteins C (*sdrC*), D (*sdrD*), and E (*sdrE*), elastin-binding protein (*ebp*), clumping factors A (*clfA*), and fibronectin-binding protein (*fnbA*) and (*fnbB*), were also found in the studied genome. These surface components have several functions, including adhesion to host cells, evasion of immunological responses, and biofilm formation [39].

In silico screening of antibiotic-resistant determinants shows that *mecA* gene, *aac*(*6* 0 )*-Ieaph*(*2* 00)*-Ia*, and a *gyrA* gene are present on the Lr2 chromosome. However, *erm(C)*, which confers resistance to lincosamide, macrolide (erythromycin), and streptogramin was found on plasmid pDLK1 (Figure 1). It is suggested that *erm(C)* is essential to this isolate for environmental adaptation [40]. It was noticed that *blaZ* gene is absent in the studied genome, and since the *blaZ* gene (β-lactamase) is carried by a transposon Tn*552* located on a large plasmid [41], it was probably eliminated due to the curing of that plasmid.

The SNPs phylogenetic relationship inferred from whole genome and core genome revealed that strain Lr2 exhibits the nearest identities to South African CA-MRSA ST22 strain 71A\_S11 (CP010940) and makes a separate clade with Indian strain CA-MRSA ST22 VB31683 (CP035671) (Figure 2). The heatmap of accessory genes comparison revealed variation in accessory genes among different strains showing gain or loss of genes as well as mobile genetic elements (MGEs) such as plasmids, prophages, and *SCCmec* element, etc. (Figure 6). This pattern of variation in genes was also observed in previous studies [22,42]. These variations can be attributed to their distinct genetic makeup, as we noticed diverse distribution of antibiotic resistance determinants and virulence factors in these strains. We observed that *blaZ*, *mecA*, and *gyrA* antibiotic resistance genes are present in most ST22 strains. However, *blaZ* was absent in strains Lr2, R20, IT1-S, and IT4-R, while *aac*(*6* 0 )*-Ieaph*(*2* 00)*-Ia* and *erm(C)* are unique to Lr2 and Indian strain VB31683. A similar pattern of

variation in virulence factors was also observed as the PVL genes (*lukF-PV* and *lukS-PV*) were only found in Lr2, VB31683, and 71A\_S11 strains. A virulence factor clumping factor B (*clfB*), was only present in Lr2 (Figure 2). These additional genes associated with antibiotic resistance and virulence are expected to make Lr2 more resistant to antibiotics and virulent.

The pan-genome analysis of MRSA ST22 strains revealed that more than 80% of genes are part of the core genome, which shows high conservancy in these strains. The core genome annotation using the Cluster of Orthologous Gene (COG) database revealed the two largest core genome categories with functions associated with metabolism and transport as well as related to information, storage, and processing. Previous studies also reported high conservancy in the *S. aureus* core genome and that the core genes are mostly associated with metabolism, replication, and information storage and processing [22,43,44]. However, genes belonging to mobilome, prophages, and transposons showed the lowest proportion. Interestingly, mobilome is the third-largest category of COGs unique to strain Lr2 (Figure 4). This abundance of mobile genetic elements in the Lr2 strain likely contributes to its increased virulence and antibiotic resistance.

## **4. Materials and Methods**

#### *4.1. Isolation and Antibiotic Susceptibility Testing*

*S. aureus* isolates were collected from hospitals in Peshawar, Rawalpindi/Islamabad, and Lahore, and 4 isolates were randomly selected for whole-genome sequencing. Of the 4 sequenced strains, P10 and R46 were found to belong to ST113 [20], Lr12 belong to a new sequence type (ST5352) [21], and Lr2 to ST22. The isolate Lr2 was collected from blood culture in a hospital in Lahore, Pakistan, in March 2019. The isolate was preliminarily identified by biochemical tests (catalase, oxidase, and coagulase) [45]. The antimicrobial susceptibility was performed by disc diffusion method as per CLSI guidelines against the following antibiotics: ampicillin (10 µg), oxacillin (1 µg), methicillin (10 µg), vancomycin (5 µg), gentamicin (10 µg), erythromycin (15 µg), streptomycin (25 µg), clindamycin (2 µg), ciprofloxacin (5 µg), chloramphenicol (30 µg), linezolid (30 µg), tetracycline (30 µg), rifampicin (5 µg), and fusidic acid (10 µg) [46]. Furthermore, resistance to methicillin was confirmed by PCR amplification of the *mecA* gene [20].

#### *4.2. Genome Sequencing, Assembly, and Annotation*

The genomic DNA (gDNA) was extracted from a fresh culture of Lr2 by Invitrogen® DNA extraction kit Cat no. K1820-01 (Thermo Fisher Scientific, Carlsbad, CA, USA). The gDNA integrity was checked by 0.7% agarose gel and quantified by Qubit 2.0 fluorometer Cat no. Q32866 (Manufactured by Tecan Austria GmbH, Grodig, Austria for Life Technologies). The gDNA libraries were prepared by Nextera XT Library Prep Kit (Illumina, San Diego, CA, USA), which was used as per instructions with a slight modification. Sequencing was performed in Illumina (HiSeq) system using a 250 bp paired-end protocol, and Trimmomatic 0.30 was used to trim adapters from raw reads [47]. The resultant reads were de novo assembled using SPAdes version 3.7 [48], and the generated contigs of Lr2 were annotated by NCBI's Prokaryotic Genome Annotation Pipeline (PGAP) [49]. Finally, the assembled contigs were reordered with *S. aureus* reference genome NCTC8325 using Mauve [50], and the genome was visualized using CGView server (http://cgview.ca/, accessed on 3 December 2021) [51].

#### *4.3. Genome-Based Characterization*

In silico epidemiological characteristics of the assembled genome were carried out using SCCmecFinder-1.2 for the identification of *SCCmec* type [52], MLST 1.8 [53] for Multilocus Sequence Typing, and *spa*Typer 1.0 [54] for *spa* typing available at Center for Genomic Epidemiology (CGE) webserver (https://cge.cbs.dtu.dk/services/, accessed on 22 August 2021).

#### *4.4. Prediction of Resistome, Virulome, and Mobilome*

The chromosomal mutations and acquired genes conferring antibiotic resistance were identified by ResFinder 4.1 at CGE (https://cge.cbs.dtu.dk/services/ResFinder/, accessed on 22 August 2021) [55]. Virulence factors in the genome were identified and annotated using the Virulence Factor Database (VFDB) at http://www.mgc.ac.cn/VFs/, accessed on 22 August 2021 [56]. The assembled genome was searched for plasmid replicons (*rep*) in PlasmidFinder 2.1 using default parameters [57]. The identified plasmid replicons (*rep*) were searched in PLSDB (Plasmid Database) for the identification of full-length plasmids, and the full-length plasmids were then BLASTed with the Lr2 genome [58]. The prophage sequences in the assembled genome were identified and annotated by PHASTER online tool [59]. The identified plasmids and prophages were also investigated for genes associated with antimicrobial resistance and virulence.

#### *4.5. Whole-Genome Single Nucleotide Polymorphism (SNP)-Based Phylogenetic Analysis*

The whole-genome SNPs phylogenetic analysis of MRSA Lr2 ST22 strain was analyzed against a set of publicly available ST22 complete genomes (*n* = 19, according to PATRIC, https://www.patricbrc.org/, accessed on 29 March 2021). The SNPs were called against the reference genome NCTC\_8325 (GenBank accession no. CP000253.1), and the maximumlikelihood tree was established using FastTree 2 tool [60] in CSI Phylogeny (https://cge.cbs. dtu.dk/services/CSIPhylogeny/, accessed on 15 April 2021) [61]. The following default settings were applied: SNP positions minimum depth 10, SNP positions relative depth 10, a minimum distance between SNPs (prune) 10, minimum quality of SNP 30, and minimum Zscore of SNP 1.96. The phylogenetic tree was visualized, and a heatmap of the presence and absence of antibiotic resistance genes and virulence factors was generated using Interactive Tree of Life (iTOL) [62].

#### *4.6. Pan-Genome and Cluster of Orthologous Genes (COGs) Analysis*

The MRSA ST22 genomes were annotated by Prokka at default parameters for estimation of pan genome [63]. The pangenome analysis and core-genome SNPs-based phylogenetic analysis were performed using the pangenome estimation module (PGM) of an in-house pipeline PanRV, which makes use of Roary [64,65]. The functional annotation of core genome of ST22 strains and COGs unique to Lr2 ST22 strain was performed using the functional annotation analysis module (FAM) of the in-house PanRV pipeline with 0.001 E-value, bit score 100, and percentage identity 70 [66,67].

#### **5. Conclusions**

This study provides important epidemiological and genomic data on a PVL-positive CA-MRSA ST22 strain Lr2 from Pakistan. The comparative analysis of the resistome shows that, unlike other ST22 strains, the strain Lr2 carries a plasmid-associated antibiotic resistance gene *erm(C)* and an aminoglycoside resistance gene *aac(6*0 *)-Ie-aph(2*00*)-Ia*. The strain also harbors several prophages with genes encoding important *S. aureus* virulence factors such as PVL and TSST. These antibiotic resistance- and virulence-associated genes possibly make this ST22 clone highly antibiotic-resistant and virulent, and these genes could be transferred to other MRSA clones through horizontal gene transfer. However, further studies are needed on a large collection of isolates to determine the molecular epidemiology, evolution, and dynamics of transmission of this clone in Pakistan.

**Author Contributions:** Conceptualization, N.U. and A.A. (Amjad Ali); formal analysis, N.U., S.N., F.A. and Z.I.; investigation, N.U. and T.R.; methodology, N.U., A.A. (Amjad Ali), A.A. (Abdulrahman Alshammari), M.A., T.B., M.R. and A.R.; project administration, A.A. (Amjad Ali); software, N.U., S.N., F.A. and Z.I.; supervision, A.A. (Amjad Ali); funding acquisition, A.A. (Amjad Ali) and A.A. (Abdulrahman Alshammari); writing—original draft, N.U.; writing—review and editing, T.B., A.A. (Abdulrahman Alshammari), M.R. and A.A. (Amjad Ali). All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was financially supported by the Higher Education Commission (HEC), Pakistan, under the NRPU project (4774) and King Saud University, Riyadh, Saudi Arabia, under the Researchers Supporting Project (RSP2022R491).

**Institutional Review Board Statement:** The study was approved by Atta-ur-Rahman School of Applied Biosciences (ASAB), NUST (Ref. # IRB-133, Dated: 6 December 2018).

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The PVL-positive CA-MRSA strain Lr2 ST22 whole-genome sequence reads has been submitted in the NCBI SRA (SRR15497842) and sequence data in NCBI GenBank (JAIGOF000000000) under the BioProject PRJNA520768.

**Conflicts of Interest:** The authors declare that they have no conflict of interest.

#### **References**