Design of a multi-epitope vaccine against six Nocardia species based on reverse vaccinology combined with immunoinformatics

Background Nocardia genus, a complex group of species classified to be aerobic actinomycete, can lead to severe concurrent infection as well as disseminated infection, typically in immunocompromised patients. With the expansion of the susceptible population, the incidence of Nocardia has been gradually growing, accompanied by increased resistance of the pathogen to existing therapeutics. However, there is no effective vaccine against this pathogen yet. In this study, a multi-epitope vaccine was designed against the Nocardia infection using reverse vaccinology combined with immunoinformatics approaches. Methods First, the proteomes of 6 Nocardia subspecies Nocardia subspecies (Nocardia farcinica, Nocardia cyriacigeorgica, Nocardia abscessus, Nocardia otitidiscaviarum, Nocardia brasiliensis and Nocardia nova) were download NCBI (National Center for Biotechnology Information) database on May 1st, 2022 for the target proteins selection. The essential, virulent-associated or resistant-associated, surface-exposed, antigenic, non-toxic, and non-homologous with the human proteome proteins were selected for epitope identification. The shortlisted T-cell and B-cell epitopes were fused with appropriate adjuvants and linkers to construct vaccines. The physicochemical properties of the designed vaccine were predicted using multiple online servers. The Molecular docking and molecular dynamics (MD) simulation were performed to understand the binding pattern and binding stability between the vaccine candidate and Toll-like receptors (TLRs). The immunogenicity of the designed vaccines was evaluated via immune simulation. Results 3 proteins that are essential, virulent-associated or resistant-associated, surface-exposed, antigenic, non-toxic, and non-homologous with the human proteome were selected from 218 complete proteome sequences of the 6 Nocardia subspecies epitope identification. After screening, only 4 cytotoxic T lymphocyte (CTL) epitopes, 6 helper T lymphocyte (HTL) epitopes, and 8 B cell epitopes that were antigenic, non-allergenic, and non-toxic were included in the final vaccine construct. The results of molecular docking and MD simulation showed that the vaccine candidate has a strong affinity for TLR2 and TLR4 of the host and the vaccine-TLR complexes were dynamically stable in the natural environment. The results of the immune simulation indicated that the designed vaccine had the potential to induce strong protective immune responses in the host. The codon optimization and cloned analysis showed that the vaccine was available for mass production. Conclusion The designed vaccine has the potential to stimulate long-lasting immunity in the host, but further studies are required to validate its safety and efficacy.

Background: Nocardia genus, a complex group of species classified to be aerobic actinomycete, can lead to severe concurrent infection as well as disseminated infection, typically in immunocompromised patients. With the expansion of the susceptible population, the incidence of Nocardia has been gradually growing, accompanied by increased resistance of the pathogen to existing therapeutics. However, there is no effective vaccine against this pathogen yet. In this study, a multi-epitope vaccine was designed against the Nocardia infection using reverse vaccinology combined with immunoinformatics approaches.
Methods: First, the proteomes of 6 Nocardia subspecies Nocardia subspecies (Nocardia farcinica, Nocardia cyriacigeorgica, Nocardia abscessus, Nocardia otitidiscaviarum, Nocardia brasiliensis and Nocardia nova) were download NCBI (National Center for Biotechnology Information) database on May 1st, 2022 for the target proteins selection. The essential, virulent-associated or resistant-associated, surface-exposed, antigenic, non-toxic, and non-homologous with the human proteome proteins were selected for epitope identification. The shortlisted Tcell and B-cell epitopes were fused with appropriate adjuvants and linkers to construct vaccines. The physicochemical properties of the designed vaccine were predicted using multiple online servers. The Molecular docking and molecular dynamics (MD) simulation were performed to understand the binding pattern and binding stability between the vaccine candidate and Toll-like receptors (TLRs). The immunogenicity of the designed vaccines was evaluated via immune simulation.
Results: 3 proteins that are essential, virulent-associated or resistant-associated, surface-exposed, antigenic, non-toxic, and non-homologous with the human proteome were selected from 218 complete proteome sequences of the 6 Nocardia subspecies epitope identification. After screening, only 4 cytotoxic T lymphocyte (CTL) epitopes, 6 helper T lymphocyte (HTL) epitopes, and 8 B cell

Introduction
The genus Nocardia, a group of obligate aerobic gram-positive actinomycetes, is a common etiological agent of nocardiosis that widely exists in the soil, decaying plants, and other organic organisms (1). In 1888, Edmond Nocard isolated an acid-resistant filamentous pathogen from the diseased cattle with symptoms of primary suppurative pulmonary infection and granulomatous abscess (2). In 1889, Trevisan dubbed this strain "Nocardia farcinica" and introduced the Nocardia genus in detail (2). To date, 120 strains have been identified, of which 54 can cause nocardiosis, and based on the widely accepted Brown-Elliott classification, the Nocardia genus has been divided into 9 groups (3,4). When immunity is significantly compromised, Nocardia can invade the body via trauma and the respiratory tract, leading to suppurative skin infection and lung infection, or even spreading to organs such as brain and kidney through hematogenous dissemination to result in a systemic infection (5). Notably, the case fatality rate of patients with disseminated infection is 16.2%-38.2% and may be over 50% if the central nervous system (CNS) is affected (6,7). In recent years, with the increase in the usage of immunosuppressive drugs, organ transplantation, and the incidence of chronic lung diseases, the incidence of Nocardia infection has shown an upward trend.
Nocardia has been reported to be capable of surviving as a facultative intracellular parasite within macrophages and evading damage by host neutrophils and monocytes (8,9). Some Nocardia infections are asymptomatic, combined with a long incubation period, making them difficult to identify in clinical specimens (10, 11). Moreover, the cooccurrence of bacterial and fungal infections will further complicate the diagnosis (12). In addition to immune escape and diagnostic difficulty, another issue worth noting is the heterogeneity in the genome size (ranging from 6 to 10 million base pairs (bps) among distinct Nocardia species (13,14). Nocardia varies from different species in antimicrobial susceptibility patterns and virulence. Currently, clinically relevant Nocardia species have been categorized into 13 antimicrobial susceptibility patterns (15). Nevertheless, the identification of Nocardia species necessitates cutting-edge equipment and extensive professional expertise, which is generally difficult to achieve in clinical practice (16,17). At present, antibiotic therapy is the mainstay of treatment for Nocardia, but its duration of treatment is very lengthy. It is generally recommended that patients with normal immune function should receive 6-12 months of antimicrobial therapy, but those with weakened immune systems or CNS dissemination should be treated for at least 12 months (18). Notably, the development of resistance in Nocardia to preferred antibiotics, such as trimethoprimsulfamethoxazole (TMP-SMX), amikacin, and third-generation cephalosporins, is becoming highly common (19,20). Thus, it is important to develop effective prevention strategies as quickly as possible.
Vaccination is one of the most cost-effective ways to protect people and communities against infectious diseases. However, no vaccine study on human nocardiosis has been reported yet. Traditional vaccine development is a costly and time-consuming procedure, and is fraught with uncertainty regarding the vaccine's specificity, allergenicity, and toxicity, especially for pathogens with complex species and strains (21). Owing to developments in structural biology, genomics, proteomics, and computational science, the field of vaccines has undergone revolutionary changes during the past decades. Reverse vaccinology (RV), as a new technique, is able to locate immunogenic antigens from the pathogen's entire genome sequence (22). It has many advantages over traditional vaccination methods (23, 24). First, conveniencethe whole process starts with the analysis of genome sequence, with no need to cultivate microorganisms. Second, safetyit can avoid the spread of pathogenic microorganisms. Third, comprehensivenessall protein antigens expressed by the pathogens in different periods and environments can be analyzed and used as candidate antigens, even when little is known about the pathogenic mechanism and immune response of the microorganisms. Moreover, RV can be combined with other methods such as immunoinformatics or various biophysical analyses to construct a novel multi-epitope vaccine. It is worth mentioning that the vaccines developed using this approach against several additional infections, including Brucella melitensis and Acinetobacter baumannii, had been tested in animal experiments, and encouraging experimental results were obtained (25,26).
In this study, we developed a multi-epitope vaccine against Nocardia infection by combining RV with immunoinformatics techniques. All the proteomes of 6 Nocardia subspecies (i.e., Nocardia farcinica, Nocardia cyriacigeorgica, Nocardia abscessus, Nocardia otitidiscaviarum, Nocardia brasiliensis and Nocardia nova) that have a much wider range of popularity were employed to identify the ideal T-cell (CTL and HTL) and B-cell epitopes to be used in the final vaccine construct. The binding mode with immune receptors, the immune effect, the population coverage, and the potential for mass production of the designed vaccine were then assessed using immunoinformatics methods. This study not only gives theoretical reference and data for future research on the Nocardia vaccine but also offers new ideas for vaccine development against other similar pathogens. To identify the essential proteins, all selected proteins were submitted to the BLAST v2.12.0+ software to perform sequence alignment between selected proteins and proteins in the Database of Essential Genes 10 (DEG 10) based on the filtering criteria of E-value<10-4, bit score>100 and sequence identity>30% (29). The DEG 10 (http://origin.tubic.org/deg/public/index.php) is a database containing frequently updated essential proteins of prokaryotes (30). Moreover, the virulence proteins and resistance-related proteins were identified via BLASTp from the Virulence Factor database (VFDB) (http://www.mgc.ac.cn/VFs/) and The Comprehensive Antibiotic Resistance Database (CARD) (https:// card.mcmaster.ca/), respectively, with the filtering criteria set to E-value<10-4, bit score>100 and sequence identity >30% in both cases (31,32). All of the above sequence data was downloaded from DEG, VFDB and CARD on May 1st, 2022.

Screening of human nonhomologous proteins
To avoid inducing autoimmune reactions, we checked the homology between the chosen proteins and the human proteome. First, we obtained the Homo species reference proteome (taxid: 9606) from the NCBI database and then built a local database using the local BLAST v2.12.0+ software based on the obtained data. Then, sequence alignment between selected essential proteins and the established local database was performed using BLASTp. Following the exclusion criteria of E-value<10-4 and bit score>100, the essential proteins with homology to the human proteome were excluded.

Analysis of antigenicity and subcellular localization
VaxiJen2.0 (http://www.ddg-pharmfac.net/vaxijen/VaxiJen/VaxiJen. html) was utilized to assess the antigenicity of the chosen proteins based on the physicochemical properties of protein amino acids (33). To ensure that the final epitopes included in the vaccine construct have better binding ability to lymphocyte receptors or antibodies, we increased the screening threshold for antigenicity to 0.5. Moreover, to validate the results, ANTIGENpro (http://scratch.proteomics.ics.uci.edu/) was also used to predict the antigenicity of chosen proteins (34). The membrane and secretory proteins are easily accessible to lymphocyte receptors and antibodies, indicating that they may serve the ideal vaccine targets. Hence, two web servers, PSORTb v3.0.2 (https://www.psort.org/psortb/) and Cell-PLoc 2.0 (http://www.csbio.sjtu.edu.cn/bioinf/Cell-PLoc-2/) were employed to predict the subcellular localization of the selected proteins. PSORTb v3.0.2, a bacterial localization prediction program, can be used to predict the cytoplasmic membrane, outer membrane, periplasm, extracellular and other subcellular localizations (35), and Cell-PLoc 2.0 can also be used to predict subcellular localizations of proteins in many organisms (36).

Prediction of transmembrane helix, peptide signal and physicochemical properties
The DeepTMHMM server (https://dtu.biolib.com/DeepTMHMM) was employed to predict the transmembrane helix, and the results suggested that its performance was superior to that of previous versions (37). Only the selected proteins with transmembrane helix ≤1 would proceed to the next screening step. The SignalIP 6.0 server (https:// services.healthtech.dtu.dk/service.php?SignalP-6.0), which can predict 5 types of signal peptides and their cleavage sites in bacterial and archaeal proteins, was used to predict the signal sequence of selected proteins, and only the proteins without signal sequence were considered for vaccine construction (38). Physicochemical properties such as molecular weight and instability index were calculated via the ProtParam server (https:// web.expasy.org/protparam/) (39). To facilitate protein purification during vaccine development, only proteins with appropriate molecular weight (<110KDa) and ideal stability (instability index<40) were selected for epitope identification. Further, the toxicity and allergenicity of proteins were predicted by ToxinPred2 and AllerTOP v.2.0 servers, respectively (40,41). Finally, the models of selected proteins were predicted via Robetta server (42), the quality of which was validated through PROCHECK (https://saves.mbi.ucla.edu/) (43), ERRAT (https:// saves.mbi.ucla.edu/) (44) and the ProSA-web (https://prosa.services. came.sbg.ac.at/prosa.php) (45) servers.

Epitopes screening 2.2.1 Prediction of cytotoxic T lymphocyte epitopes
The CTL epitopes were restrictive in binding to MHC class I molecules. To design vaccines with a wide population coverage, we selected human leukocyte antigen (HLA) class I alleles with a high frequency of expression in this study. Each HLA-I allele covered more than 5% of the global population, which was calculated by the Population Coverage tool in the Immune Epitope Database (IEDB) (http://tools.iedb. org/mhci/). There were high frequency HLA-I alleles (HLA-A*01:01, . The NetCTLpan1.1 server (https://services.healthtech.dtu.dk/service.php?NetCTLpan-1.1) was utilized to predict CTL epitopes by combining MHC-binding peptide prediction, proteasome cleavage site prediction and TAP transport efficiency prediction (46). In this server, the threshold for epitope identification was set to 1%, which indicated that CTL epitopes had good binding to HLA-I alleles. . In this study, the epitopes with a percentile rank<1%, which are considered stronger binders, were chosen for further analysis.

Prediction of B-cell epitopes
B-cell epitopes are essential for eliciting humoral immune responses and neutralizing antibody responses. In this study, we utilized three servers to identify B-cell epitopes, and only overlapping epitopes were chosen for further analysis. The first server is ABCpred (https://webs.iiitd.edu.in/raghava/abcpred/), which uses the artificial neural network (ANN) to predict epitopes with an accuracy of 65.3% (48, 49). The second server is Bepipred2.0 (http://tools.iedb.org/bcell/), which employs a random forest method to predict B-cell epitopes from antigen sequences (50). The third server is ElliPro (http://tools.iedb.org/ellipro/), which was used to predict potential B-cell epitopes based on the 3D structure of proteins (51).

Antigenicity, toxicity, allergenicity, and surface accessibility of selected epitopes
The antigenicity of selected epitopes was predicted using the VaxiJen 2.0 server, and the epitopes with antigenicity>0.5 were further submitted to the ToxinPred2 (https://webs.iiitd.edu.in/ raghava/toxinpred2/algo.html) and AllerTOP v.2.0 (https://www. ddg-pharmfac.net/AllerTOP/) servers to check for their toxicity and allergenicity, respectively. The screening threshold of the ToxinPred2 server was set to 0.6, and its performance appeared to be superior to that of other servers. AllerTOP v.2.0 predicts the antigenicity of epitopes using the k closest neighbors (kNN) algorithm with an accuracy of 88.7% (41). The accessibility of epitopes on protein surfaces was observed by the surface mode of the PyMOL v2.4 software (52). Based on the above prediction results, all epitopes located in signal peptides, transmembrane helices, and intramembrane regions were excluded.

Immunogenicity analysis of CTL epitopes and validation of the binding of HLA class I alleles
Immunogenicity is the capacity of an agent, such as an epitope or antigen, to stimulate an immune response within the body. In this study, the Class I Immunogenicity tool (http://tools.iedb.org/ immunogenicity/) in the IEDB was employed to evaluate the immunogenicity of CTL (the higher the immunogenicity score, the greater the probability of triggering a host immune response) (53). TepiTool (http://tools.iedb.org/tepitool/) from IEDB was applied to validate HLA class I alleles that bind to CTL epitopes using a consensus algorithm that integrates four methods (54, 55), and it is generally believed that epitopes with IC50 ≤ 500 nm have moderate binding to HLA class I molecules.

Vaccine construction
Compared to conventional vaccines, multi-epitope peptide vaccines are weakly immunogenic. Human beta-defensins (HBDs) are small molecular cationic antimicrobial peptides produced by epithelial cells that exhibit antibacterial and immunomodulatory activity, which is considered to be part of the innate and adaptive immune response (68). Hence, to enhance the immunogenicity of the vaccine candidate, b-defensin 3 was added to the vaccine's Nterminus as an adjuvant. Further, for avoiding the formation of newly connected epitopes and ensuring the independent immunogenicity of each epitope, appropriate connectors were added between the epitopes. Then, the EAAAK linker was utilized to connect adjuvants with the vaccine (69), and the CTL, HTL, and B cell epitopes were linked by AAY, GPGPG, and KK linkers, respectively. The HEYGAEALERAG linkers, acting as cleavage sites for the proteasomal and lysosomal systems, were chosen as spacer linkers to connect CTL, HTL, and B cell epitopes (70). The PADRE sequence, which can facilitate antibody responses, was linked to the C-terminal of the vaccine using the EAAAK linker (71).

Features of vaccine construct
The VaxiJen2.0 server was used to predict the antigenicity of the vaccine candidate, which was then confirmed by the ANTIGENpro server. The allergenicity and toxicity were checked by AllerTOP v.2.0 and the ToxinPred2 server, respectively. Thereafter, the ProtParam tool was employed to check the physiochemical properties of the designed vaccine, including molecule weight, instability index, the number of amino acids, theoretical isoelectric point (PI), aliphatic index, and GRAVY and half-life. The potential for soluble expression of the vaccine candidate in E. coli was evaluated by the SOLpro server (http://scratch.proteomics.ics.uci.edu/) with an accuracy of about 74% (72). Also, the DeepTMHMM and SignalP6.0 servers were employed to predict possible transmembrane helices and signal peptides in the vaccine, respectively.

Secondary and 3D structure predictions, 3D structure refinement and validation
The secondary structure of the vaccine was predicted by the PSIPRED 4.0 server (http://bioinf.cs.ucl.ac.uk/psipred/) with an accuracy of 84.2% (73). The putative 3D structure of the vaccine was generated using the Robetta server (http://www.robetta.org/), which executes comparative modeling if a suitable template is in the database and executes ab initio structure prediction otherwise. Then, the primary structures were submitted to GalaxyRefine for refinement (74). The GalaxyRefine server (https://galaxy.seoklab.org/cgi-bin/ submit.cgi?type=REFINE) can improve the structure quality by rebuilding side chains and relaxing the overall structure, and its main parameter MolProbity was used to evaluate the quality of refined models (the lower the MolProbity score, the better the quality of the model). Furthermore, the quality of refined models was validated by the PROCHECK -generated Ramachandran plot, the ERRAT score, and the Z score generated ProSA-web.

Disulfide engineering of vaccine
To assist protein folding and improve the stability of the protein structure, the Disulfide by Design (DbD) v2.13 server (http://cptweb. cpt.wayne.edu/DbD2/) was employed to conduct disulfide engineering of the designed vaccine, which could check residue pairs with the ability to form disulfide bonds if the amino acid mutated to cysteine (75). The potential residue pairs were then selected based on the criteria of chi3 value -87°or +97°± 30 and energy value< 2.2 kcal/mol. Subsequently, the PROCHECK, ERRAT, and ProSA-web servers were used to assess the quality of the model with selected disulfide bonds.

Vaccine-Toll-like receptors docking
The molecular docking of the designed vaccine with the receptor can help clear the binding mode between two molecules. Toll-like receptor 4 (TLR4), a type I transmembrane protein expressed in a variety of tissues including monocytes and macrophages, plays a key role in activating innate immunity by directly recognizing pathogenrelated molecular patterns (PAMP), such as lipopolysaccharide (LPS), lipoteichoic acid (LTA), double-stranded RNA (dsRNA), etc (76). TLR2, another member of the TLR family, can recognize various microbial components from gram-positive bacteria such as lipopeptide, peptidoglycan, and fatty acids, and regulate the phenotype and function of dendritic (DC) cells, T cells, natural killer (NK) cells, macrophages and other immune cells (76,77). Hence, the molecular docking between the multi-epitope vaccine and TLR2 (PDB ID:2Z7X) or TLR4 (PDB ID:4G8A) was realized by the ClusPro 2.0 server (https://cluspro.org/help.php), which ranks the cluster of docking complexes according to their central and lowest energy scores (78). Then, the best docking model was submitted to the HADDOCK 2.4 server (https://wenmr.science.uu.nl/haddock2.4/) for refinement, with the parameters set to default (79). Besides, the PDBsum server (http://www.ebi.ac.uk/thornton-srv/databases/cgibin/pdbsum/GetPage.pl?pdbcode=index.html) was used to analyze the interacting residues between the vaccine and TLRs (80).

Molecular dynamics simulation and analysis 2.9.1 MD simulation
The MD simulation was run using the Gromacs v2022.1 software to understand the dynamics and stability of vaccine-TLR docking complexes (81), and all MD simulations and corresponding analyses were run in triplicate for ensuring the accuracy of the results. During the initial phase, the AMBER14SB_PARMBSC1 force field was added to generate the coordinate file and topology of the vaccine-receptor complex (82). The complex was then solvated by adding transferable intermolecular potential 3P (TIP3P) water molecules, followed by the addition of counter ions (Na+ and Cl-) to neutralize the charge of the simulated system. Subsequently, energy minimization was performed until reaching the desired threshold (the maximum force ≤1000.0 kJ/ mol/nm). To achieve a stable temperature and pressure environment, NVT and NPT simulations were run for 200ps and 1ns, respectively, and the system's temperature and pressure were equilibrated to 310K and 1atm, respectively. The modified Berendsen thermostat was employed to regulate the temperature of the system (83), while the Parinello-Rahman barostat was used to equilibrate the pressure (84). Subsequently, the MD simulation was run for 100 ns with an isothermal and isobaric system. At last, the LINCS algorithm was used to constrain the hydrogen bonds (85), and the Particle-Mesh Ewald summation scheme was executed to calculate long-range electrostatic interactions (86). Finally, the VMD v1.93 software was used to generate the movie of the MD simulations (87).

MM-PBSA calculation
To calculate the binding free energy estimates between the vaccine and TLRs, the last 10 ns of the trajectory was extracted and further submitted to the gmx_MMPBSA v1.56 tool (88). Besides, the energy contribution of all residues within 4 Å of the receptor and ligand was also analyzed using the gmx_MMPBSA v1.56 tool. In this study, the MM-PBSA method was adopted to predict the binding free energy of the complex, whose mathematical explanation is:

Principal component analysis and dynamical cross-correlation analysis
The principal components (PCs) are eigenvectors that determine the direction of motion, and the extent of residual motion during the MD is defined by the corresponding eigenvalues. In this study, PCs were mainly projected to two dimensions, the principal component 1 (PC1) and the principal component 2 (PC2). Then, the DCC analysis was performed to understand the correlation between the motions of the residues of the vaccine chain and the motions of the residues of the TLRs chain (89). The last 50 ns of the trajectory was extracted and used for the analyses of PCA and DCC, which were operated using the Bio3D package in R version 4.21 (90).

Immune simulation
The immune simulation of the designed vaccine was performed using the C-ImmSim server (https://kraken.iac.rm.cnr.it/C-IMMSIM/ index.php?page=0), which can simulate the immune response generated by a mammalian thymus (T cells), bone marrow (lymph and bone marrow cells), and a lymphoid organ (91). The simulation results could help us estimate the immunogenicity of the vaccine and the immune response it might induce in the human body. According to previous studies, the administration scheme adopted in this study was 1000 vaccine units at a time, once every four weeks, for a total of three injections. The total number of steps was set to 1050, with 1, 84, and 168 time steps for the first to the third injection, respectively. The host MHC molecule combinations were HLA-A*01:01, HLA-B*07:02, and HLA-DRB1*01:01, and other parameters were kept unchanged.

Population coverage
The proportion of the population that will respond to the vaccine (population coverage) is another important indicator to consider in the vaccine design process. Distinct epitopes bind to different HLA alleles, and the distribution frequency of distinct alleles varies across different regions and ethnicities. Based on the distribution frequency of different HLA alleles and the HLA alleles bound to the vaccine, the global population coverage of the vaccine was calculated using the Population Coverage tool (http://tools.iedb.org/population/) in IEDB (92).

Codon optimization and cloning
To improve the efficiency of vaccine expression in the selected expression host, the vaccine constructs were submitted to the online tool VectorBuilder (https://www.vectorbuilder.cn/tool/codonoptimization.html) for codon optimization, and E. coli (Strain K12) was selected as the expression host in this study. The codon adaptation index (CAI) and percentage GC content, generated by the VectorBuilder server, were utilized to evaluate protein expression levels. The values of CAI vary from 0 to 1, and a high CAI value indicates a high-level expression of genes. The recommended range of GC content varies from 30% to 70%. Then, the Xhoi and BamHI restriction endonuclease sites were attached to the vaccine's C terminal and N terminal, respectively. At last, the final sequence of the vaccine was inserted into the pET-28a (+) vector using the GenSmart tool (https://www.genscript.com/gensmart-design/).

Prediction of the vaccine mRNA secondary structure
The Transcription and Translation Tool (http://biomodel.uah.es/ en/lab/cybertory/analysis/trans.htm) was applied to obtain the mRNA sequence of the vaccine. Two online servers, Mfold v2.3 (http://www.unafold.org/mfold/applications/rna-folding-form-v2. php) (93) and RNAfold (http://rna.tbi.univie.ac.at/cgi-bin/ RNAWebSuite/RNAfold.cgi) (94) were used to predict the secondary structure of the vaccine mRNA. The main output result was the minimum free energy (DG Kcal/mol), and the lower the minimum free energy, the more stable the folding structure of mRNA.

Pan-proteomics analysis of primary data
A total of 218 proteomes available for the 6 representative Nocardia subspecies were obtained from the NCBI database. The detailed information was shown in Supplementary Data Sheet 1. The results of pangenome analysis predicted by the BPGA tool revealed that 1,336 core proteins were found in the 218 proteomes (Supplementary Data Sheet 2).
The flow chart for the design of a multi-epitope vaccine against 6 Nocardia subspecies was presented in Figure 1.

Selection of essential, virulence, resistance, and non-homologous proteins
The  Sheet 7), of which 3 ideal proteins were identified after analyses of antigenicity, subcellular localization, transmembrane domains, signal sequence, allergenicity, toxicity, and physicochemical properties. From Table 1, it could be seen that all these 3 proteins were antigenic, non-allergic, non-toxic, nonhomologous with human proteins, stable and of a suitable molecular FIGURE 1 Workflow designed for 6 common Nocardia subspecies vaccines. weight. Also, all the 3 proteins were essential for bacteria, and two of them were associated with antibiotic resistance, and one with bacterial virulence. After sequence alignment with the NCBI database using BLASTp, it was further found that one of the 3 proteins was a transglycosylase domain-containing protein, one was a penicillinbinding protein, and the other was type VII secretion-associated serine protease myosin-like protein. The structures of the target proteins were shown in Supplementary Data Sheet 8; Figures S1-S3.

Epitopes screening
After the prediction of antigenicity, immunogenicity, toxicity, allergenicity, and the ability to induce the generation of INF-g, IL-4 and IL-2. Then, we identified 4 ideal CTL epitopes from 356 epitopes (Table 2), 6 desirable HTL epitopes from 305 epitopes (Table 3), and 8 appropriate B cell epitopes from 61 epitopes for designing a multiepitope vaccine ( Table 4). All selected epitopes were located on the surface of the target protein and visualized by PyMOL v2.4 ( Figure 2).

Molecular docking between T cell epitopes and MHC molecules
The docking results of Autodock Vina revealed that the binding affinity of CTL epitopes to MHC I alleles varied between -9.1 and -8.4 kcal/mol (Figure 3), whereas the binding affinity of HTL epitopes to MHC II alleles varied between -6.6 and -4.4 kcal/mol. According to relevant research, the binding energy is often used to determine the likelihood that a ligand will bind to a protein target, with a larger negative binding affinity value indicating a stronger interaction between the receptor and ligand. From the above, it was suggested that T cell epitopes had a good affinity for MHC allele molecules in this study. Furthermore, the predicted results of Ligplot+ v2.2 also indicated that there were strong interactions between T-cell epitopes and MHC allele molecules, as detailed in Supplemental File 8 Tables S1, S2. The quality verification results of MHC molecules modeled by the SWISS-MODEL server are shown in Supplementary Data Sheet 8; Table S3.

Vaccine construct and features
The final vaccine sequence is shown in Figure 4A. According to the projected results of VaxiJen2.0 and ANTIGENpro, the antigenicity scores of the vaccine were 1.1966 and 0.9672, respectively. The results from the ToxinPred2 and AllerTOP v.2.0 servers revealed that the vaccine was non-toxic and non-allergic. Analysis of physicochemical characteristics revealed that the vaccine construct was composed of 392 amino acids, with a molecular weight of 40.30 kDa, and its theoretical pI was 9.47. The instability index of the vaccine was 26.07, indicating a stable nature. The aliphatic index of the vaccine was 55.64, suggesting that the vaccine was thermostable. The GRAVY value of the vaccine was -0.677, which demonstrated a hydrophilic nature. The developed vaccine had an estimated half-life of 30 hours in mammalian reticulocytes, over 20 hours in yeasts, and over 10 hours in E. coli. According to the prediction results of the SOLpro server, the probability of the   The The plot of the binding patterns of CTL epitopes with MHC class I molecules. The MHC I allele molecules were presented in blue and the CTL epitope peptides were presented as magenta.  3.6 Secondary and tertiary structure predictions, tertiary structure refinement and validatio The secondary structure of the vaccine, as analyzed by the PSIPRED 4.0 server, was composed of 65.05% coil, 19.39% a-helix and 15.56% bstrand ( Figure 4B). The Robetta server was employed to generate the primary 3D structure of the vaccine, which was further refined by GalaxyRefine( Figure 5A). The best refined vaccine model, with a MolProbity of 1.897, was validated by the PROCHECK, ERRAT, and ProSA-web servers. The Ramachandran plot generated by the PROCHECK server revealed that the vaccine had 90.3% residues in the most favored regions, 8.8% residues in additional allowed regions, and only 0.3% residues in generously allowed regions ( Figure 5B). The ERRAT score of the vaccine was 93.017( Figure 5D), and the Z-score generated by the ProSA-web server was -7.3( Figure 5C). The PROCHECK server was employed to examine the stereochemistry of protein structures in detail. Over 90% of residues in the most favored regions suggested a good quality of the vaccine model. The ERRAT server was then utilized to distinguish correct and incorrect regions in the protein structure based on characteristic atomic interactions, and a high ERRAT score indicates a good quality of the structure. The Z score generated by the ProSA-web server was used to indicate the overall model quality, and the lower the Z score, the better the model quality. All of the above analysis results suggested that the developed vaccine had a reliable and stable 3D structure.

Disulfide engineering of vaccine
The inter-and intra-chain disulfide bonds were evaluated by the DbD2 server, which reported 52 pairs of amino acid residues with the ability to form disulfide bonds. Given that 90% of native disulfide bonds have an energy value< 2.2 kcal/mol, 2.2 kcal/mol was set as the threshold for screening disulfide bonds in this study. According to the energy value and chi3 value, only 3 residue pairs (i.e., LEU6-TYR9, GLU28-ARG43, and GLY234-ALA376) were chosen and they were not falling into any of the epitope regions, which could minimize the effect of disulfide engineering on vaccine immunogenicity (Supplementary Data Sheet 8; Figure S6A). Analyses of the Ramachandran plot, ERRAT score, and Z-score were shown in Supplementary Data Sheet 8 ( Figure S6B-S6D), which indicated that the model with disulfide was of good quality.

The analysis of molecular docking of vaccine-Toll like receptors
To assess the binding affinity between the vaccine construct and TLRs (TLR2 and TLR4), molecular docking was performed using the Cluspro2.0 server. The server generated a total of 30 candidate models of the vaccine-TLR2 complex with different binding energies, from which the top-ranked model with the lowest binding energy score of -837.3 was selected for further refinement. Also, a total of 30 model complexes of vaccine-TLR4 were generated, and the top-ranked model with the lowest binding energy score of -835.1 was chosen. The HADDOCK 2.4 server was then used to refine the vaccine-receptor complex, and the HADDOCK scores of the vaccine-TLR2 complex and vaccine-TLR4 complex were -194.9 ± 1.0 and -243.4 ± 1.7, respectively. The low HADDOCK scores revealed the good binding between the vaccine and TLRs. The analysis results of the PDBsum server revealed that there were 14 hydrogen bonds and 2 salt bridges formed between the vaccine and TLR2, and 12 hydrogen bonds and 7 salt bridges formed between the vaccine and TLR4 (Supplementary Data Sheet 8; Table S4). For the vaccine-TLR2 complex, the amino acid residues of the vaccine involved in salt bridge formation were GLU174 and ASP180, and the residues involved in hydrogen bonds formation were ARG12, PRO89, VAL92, SER94, LYS79, TYR213, GLU174, TYR105, TYR225, ARG103, and ASP180 ( Figure 6A). For the vaccine-TLR4 complex, the amino acid residues of the vaccine involved in salt bridge formation were LYS98, ASP96, ARG84, LYS79, ASP77, and GLU224, and the residues involved in hydrogen bonds formation were LYS98, ASP96, HIS223, ARG84, LYS79, ASP77, and LYS51 ( Figure 6B). The above results suggested that there were strong interactions between the designed vaccine and TLRs (TLR2 and TLR4).

MD simulation
MD simulation was performed to understand the stability of vaccinereceptor complexes at the atomic level using Gromacs v2022.1. A primary analysis of trajectories was performed, including the calculation of the Backbone root-mean-square deviation (RMSD), the radius of gyration (Rg), root mean square fluctuation (RMSF), total A B FIGURE 6 The plot of visualized analysis and the interacting amino acid residues of the docking complex. The TLR2 receptor is colored lightblue, the TLR4 receptor is colored palecyan, and the vaccine is colored limon. Frontiers in Immunology frontiersin.org solvent accessible surface area (SASA), and the number of hydrogen bonds. RMSD reflects the conformational fluctuations of vaccinereceptor complexes based on the initial structure. For the vaccine-TLR2 complex, the RMSD of the first and third simulations reached relative stability around 10 ns, but the RMSD of the second simulation reached stability around 15 ns ( Figure 7A). Overall, the RMSD curves for the three simulations have similar trends. The mean RMSD values of the three simulations were 0.826 ± 0.101nm, 0.821 ± 0.139nm, and 0.850 ± 0.140nm, respectively. For the vaccine-TLR4 complex, the mean RMSD values of the three simulations were 0.791 ± 0.085 nm,1.223 ± 0.190 nm, and 1.182 ± 0.356 nm, respectively. In the first and second simulations, the RMSD reached equilibrium around 10ns and 20ns, respectively. In the third simulation, RMSD reached the first equilibrium around 10ns, fluctuated again at 40 ns, and reached the second equilibrium around 45 ns ( Figure 7B). Rg characterizes the compactness of a protein structure. The Rg value of the vaccine-TLR2 complex decreased at first and then remained at a constant level at about 15ns in all three simulations ( Figure 7C). For the vaccine-TLR4 complex, the Rg values of the second and third simulations were lower than those of the first simulation, indicating a more compact structure of the complex during the second and third simulations ( Figure 7D). The analysis of Rg also showed that the vaccine-TLR4 complex underwent two conformational changes before reaching equilibrium in the third simulation ( Figure 7D). Overall, the results of Rg indicated that the structures of both the vaccine-TLR2 complex and vaccine-TLR4 complex were compact and stable over time during the three MD simulations. The RMSF value reflects the flexibility of a single amino acid of a protein throughout the MD simulation. RMSF analysis revealed that the residues 380-392 of the vaccine chain and the residues 27-100 of TLR2 in the vaccine-TLR2 complex were highly flexible in all three simulations (Supplementary Data Sheet 8; Figure S7A and S7C). For the vaccine-TLR4 complex, it was found that the residues 250-360 of the vaccine chain were highly flexible in the first simulation, while the RMSF of the vaccine had lower values in the other two simulations. Residues 250-360 did not interact with the receptor and the vaccine structure is not compact enough resulting in high flexibility of the residues during the first simulation (Supplementary Data Sheet 8; Figure S7B). Meanwhile, residues 27-300 of the TLR4 had higher flexibility in the third simulation and the RMSF of TLR4 was relatively lower than that of the vaccine in the other two simulations (Supplementary Data Sheet 8; Figure S7D). The analysis of the total solvent accessible surface area showed the change in the complex surface area with respect to time. The curve trend of the SASA plot was opposite to that of RMSD indicating that the better interaction between receptor and ligand led to a lower surface area of the protein complexes ( Figures 7E, F). The results of the analysis of SASA once again revealed the stable binding between the vaccine and TLRs (TLR2 and TLR4). The vaccine-TLR2 complex presented about 10 hydrogen bonds at the beginning of the simulation. After 20 seconds, the number of hydrogen bonds formed in the vaccine-TLR2 complex increased and then maintained relative stability in all three simulations ( Figure 7G). During the three simulations, the number of hydrogen bonds of the vaccine-TLR4 complex increased from 10 at the beginning to about 25 hydrogen bonds at the end of the simulation ( Figure 7H). The formation of hydrogen bonds suggested that the vaccine had a good affinity for TLRs. Finally, the movies of the MD simulations generated by VMD v1.93 were shown in the Supplementary Data Sheet 9-14.

MM-PBSA calculation
The tracks of 100ps extracted from each interval during the final 10ns simulation period were processed with the calculation of MM-PBSA and several binding interaction energies, including van der Waals, electrostatic, polar solubility, solvent accessible surface area (SASA), and binding energy (DG binding ) of vaccine-receptor complexes. The binding energy of the vaccine-TLR2 complex was -161.07 kJ/mol, -208.24 kJ/mol, and -148.16 kJ/mol, respectively. For the vaccine-TLR4 complex, the binding energy were -203.16 kJ/mol, -176.99 kJ/mol, and -240.85 kJ/mol in three simulations, respectively. The detailed energy analysis presented in Table 5 revealed that the vaccine had a good affinity  for both TLR2 and TLR4 molecules. For the vaccine-TLR2 complex,  residues (LYS79, ALA92, and ARG232) in the vaccine chain contribute significantly to vaccine binding to the TLR2 receptor, particularly residues (ARG232) (Supplementary Data Sheet 8; Figure S8A, S8C and S8E). For the vaccine-TLR4 complex, MM/PBSA decomposition analysis showed that residues (PRO69, LYS77, ARG84, PRO89, SER99, TYR105 and LYS389) in the vaccine chain contributed to the binding of the vaccine to the TLR4, with a greater contribution from the residue (ARG84) (Supplementary Data Sheet 8; Figures S8B, S8D, S8F). However, the residue (GLU224) negatively affected the binding of the vaccine to TLR2 and TLR4 receptors. Overall, the results of all MM-PBSA calculations indicated that the vaccine can stably bind to the TLR2 and TLR4.

PCA and DCC analysis
To study and evaluate the association between different protein conformations, PCA was conducted for the trajectories of all C atoms in 50-100ns in the simulated system. PCA translates atomic coordinates into a limited number of PCs (commonly 2 or 3 PCs) that signify the directions in which the structure set has the greatest collective variability. The PCA scatter plots revealed that the color of instantaneous conformations changed from blue to red over time. In the direction of PC1, the successive color shifts (from blue to white to red) indicated that the conformation in both systems gradually evolved to another state during the simulation. PC1 and PC2 characterized 40.02%, 45.89%, and 49.24% of the system's motility with respect to the vaccine-TLR2 complex ( Figure 8A Figure S10 depict the DCCM matrices of the vaccine-TLR2 and vaccine-TLR4 complexes. The color gradients from blue to white to red in the figures correspond to the correlation coefficients of -1, 0 and 1, respectively, and a deeper shade of the same color represents a stronger correlation. For example, the correlation coefficient of light blue is -0.5 and that of light red is 0.5. DCC analysis of the vaccine-TLR2 complex in three MD simulations showed that the residues surrounding 1 to 51 (550 to 600 in a plot) and 151 to 231 (700 to 780 in the plot) from the vaccine were positively correlated with the residues surrounding 200 to 548 (200 to 548 in a plot) and 300 to 400 (300 to 400 in the plot) from the TLR2 chain, respectively. ( Figure 8C; Supplementary Data Sheet 8 and Figure S10A and S10B). For the vaccine-TLR4 complex, the residues surrounding 9 to 79 (610 to 680 in the plot) and 99 to 199 (700 to 800 in the plot) were positively correlated with the residues surrounding 500 to 601(500 to 601 in the plot) and 200 to 400 (200 to 400 in the plot) from the TLR4 chain, respectively ( Figure 8F; Supplementary Data Sheet 8 and Figure S10C and S10D).

Immune simulation and population coverage
The immune stimulation of the vaccine was carried out using the C-ImmSim web server, which predicts whether the vaccine can induce significant adaptive immunity in the human body. Following the first injection, lower lgM antibody titers, a B-cell isotype lgM population, and a slight B-cell isotype lgG population were detected. After the second and third injections, higher levels of IgM antibody with IgG1, IgG1+IgG2 and lgG+lgM as well as a substantial number of memory B cells were found ( Figures 9A, B). In addition to the B-cell, the plasma B cells, the total, memory and active TH-cell (helper T cell) populations also increased during the secondary and tertiary immune responses ( Figures 9B-E). Besides, TC cells (cytotoxic T lymphocytes) were observed during the immune simulation, with the maximum TC cell count exceeding 1,150 cells per mm 3 ( Figure 9F). Further, the NK cells and macrophages were activated after exposure to the vaccine ( Figures 9G, H). Figure 9I showed that the vaccine was able to activate a wide range of cytokines, including IFN-g, IL-2, and IL-4. In summary, the designed vaccine was capable of inducing the generation of a large number of immune cells and a variety of cytokines, indicating that it can successfully induce adaptive immunity and memory immunity in the host body.
Population coverage analysis suggested that 89% of the world's population might respond to the designed vaccine. More specifically, its population coverage in Europe, the United States, and China reached 93.47%, 88.91%, and 86.63% (Supplementary Data Sheet 8; Table S5), respectively. The results of immunization simulation and population coverage demonstrated that the developed vaccine was highly immunogenic and able to protect the vast majority of the global population.

Codon optimization and cloning
To ensure high-level expression and easy production of the vaccine, codon optimization was performed via the VectorBuilder server. The cDNA length of the optimized vaccine was 1,176 bps and Both the CAI value of the vaccine (i.e., 0.89) and the GC content (i.e., 59.46%) were desirable, indicating a high gene expression potential and excellent expression ability in E. coli (K12 strain). Then, the GenSmart tool was used to clone the cDNA sequence of the vaccine into the pET28a(+) vector ( Figure 10A).

Secondary structure of vaccine mRNA
According to the RNAfold server, the best secondary structure of the vaccine mRNA had a minimum free energy of -474.30 kcal/mol ( Figure 10B), whereas the centroid secondary structure had a minimum free energy of -407.80 kcal/mol ( Figure 10C). The minimum free energy of the optimal secondary structure of the vaccine mRNA was determined to be -413.80 kcal/mol by the mFold v2.3 server (Supplementary Data Sheet 8; Figure S11).
Because mRNA is more stable when its minimal free energy is lower, the vaccine was predicted to be stable after expression in vivo.

Discussion
So far, Nocardia is still a significant threat to aquaculture animals and immunocompromised people, such as the elderly, transplant recipients, HIV-infected individuals, malignant tumor patients, and patients treated with glucocorticoids (4, 95-97). Moreover, Nocardia infections have also been found in people with normal immune function (98,99). Due to the elevated incidence and fatality rate of nocardiosis, and increased Nocardia resistance to existing antibiotics, the need of developing a vaccine to mitigate Nocardia infections has been underscored.
Because there is a large number of Nocardia species and many of them are pathogenic to humans, our goal was to develop a broad- spectrum multi-epitope vaccine that can protect against the majority of Nocardia, but it is a huge challenge for traditional vaccinology. Compared with traditional techniques, RV is capable of performing batch analysis of thousands of biological parameters within a shorter time while ensuring more reliable results, thereby opening up a novel door to the research and development of vaccines (23, 24). To enable the designed vaccine to take effect on most pathogenic Nocardia strains, 218 complete proteome sequences of 6 common Nocardia species were downloaded from NCBI, and the proteins containing 50% similar sequences from all strains were included for further analysis. Then, less significant proteins, such as the human proteome, essential genes, virulent factor proteins, antibiotic-resistant proteins, and human interacting proteins, were subtracted from several other databases. Besides recognizing the virulent and antibiotic-resistant factors, the identification of potential vaccine candidates must also focus on checking their antigenicity, allergenicity, subcellular localization, and other physicochemical properties, such as the aliphatic index, GRAVY, and instability index. Finally, 3 outer membrane proteins, including one virulence factor protein and two antibiotic resistance proteins, were selected for the identification of epitopes. Epitope selection is a critical step in the development of multiepitope vaccines. In this study, the antigenicity, toxicity, and allergenicity of all predicted epitopes were evaluated, and both the T-cell (CTL and HTL) and B-cell epitopes were included in the multiepitope construction, which can avoid unsatisfactory immune effects due to insufficient immune response that is induced by T-cell epitopes or B-cell epitopes alone while ensuring safety. It is well known that CTL cells can only recognize the antigens presented by MHC class I molecules, while HTL can only recognize the antigens presented by MHC class II molecules (100). Therefore, the capacity of the epitope to bind to HLA alleles is a crucial predictor of whether it can effectively induce an immune response or not. It is worth mentioning that we performed molecular docking between all ultimately selected T-cell epitopes and HLA alleles, and the docking results suggested that the affinity between them was strong, indicating that the epitopes could be effectively presented via MHC molecules. IL-2 is indispensable for the proliferation as well as differentiation of naive T cells into effector T cells (58). Both IL-4 and INF-g are important activators of T helper cells, and IFN-g also participates in the activation of macrophages (56,57). The potential of HTL epitopes to induce the IL-2, IL-4, and INF-g generation was examined in this study. With all that mentioned, the eventually-selected epitope has the potential to construct an effective multi-epitope vaccine that is characterized by few side effects.
Adjuvants can not only improve the immunogenicity of multiepitope vaccines but also help induce a more durable immune response (101). In this study, the b-defensin, which can attenuate toxicity and enhance systemic antibody responses when administered together with antigens, was incorporated into the vaccine construct as an adjuvant (102). It is well known that TLRs can recognize different PAMPs and play an indispensable role in the innate and adaptive immune response (77). The binding affinity with TLRs is an important reference factor in evaluating the immune effect of a multi-epitope vaccine candidate. From the docking results, we observed stable interactions and high affinity between the multi-epitope vaccine and TLRs (TLR4 and TLR2), indicating that the vaccine candidate can effectively trigger the body's immune response. Furthermore, the docked complexes were molecularly simulated to assess the stability and agility of the vaccine-TLR complexes in a simulated biological environment. In the three MD simulations of the vaccine-TLR2 complex, the RMSDs all remained stable in the last 80% of the simulation period, showing that the structure of the vaccine-TLR2 complex was relatively stable. For the trend of RMSD of vaccine-TLR4 complex, although there were differences in the three simulations, RMSDs all reached equilibrium at 50ns in the three simulations. The Rg values of two vaccine-TLR complex systems also tended to stabilize after 50 ns, which further confirmed the stability of the complex. In all MD simulations, there were quite negative correlations between RMSD and SASA, which verifies the good binding between vaccine and receptor (TLR2 and TLR4) on the other hand. In all MD simulations, more hydrogen bonds were formed in both the vaccine-TLR2 and vaccine-TLR4 complexes, demonstrating high affinity of the vaccine chain to the TLR2 and TL42 chains. The MM-PBSA analysis indicated that very little energy was required to keep the two complexes stable. The DCC analysis revealed the existence of positively correlated contacts between vaccine-TLR2 amino acid residues and vaccine-TLR4 amino acid residues, indicating that the associations between the vaccine and TLRs (TLR2 and TLR4) were good. According to the above analysis, the designed vaccine is able to firmly bind to TLRs and successfully transport into the body in the natural environment, thereby inducing effective immunity. In the process of immune simulation, we observed the proliferation of TC, TH and B cells, and the production of lgA, lgG, INF-g, IL-2 and other cytokines with the injection of the vaccine. In addition, memory B cells were also observed following repeated exposures to the vaccine, which suggested that the vaccine was capable of inducing not only cellular and humoral immunity but also long-term immunity against Nocardia. Finally, to evaluate the feasibility of mass production, codon optimization and cloning were implemented for the vaccine construct, and the values of CAI and GC contents suggested that the designed vaccine could be well expressed in the E. coil K12 strain. Besides, the lower minimum free energy scores of the secondary structure of vaccine mRNA suggested that the vaccine was expected to behave stable in the body.

Conclusion
Nocardiosis poses a major public health threat due to the expansion of susceptible populations, underdiagnosis and lack of adequate awareness. In this study, we designed an antigenic, non-allergic and non-toxic multi-epitope vaccine against Nocardia using RV and immunoinformatics techniques. The results of molecular docking and MD simulation show that the designed multi-epitope vaccine has a strong affinity for different innate immune receptors of the host, and that the vaccine-receptor complex is dynamically stable in the natural environment. In addition, the designed vaccine has the potential to induce humoral immunity and immunity. Finally, the codon optimization and cloned analysis show that the vaccine is feasible for mass production. To sum up, the designed vaccine may have a good protective effect on Nocardia, but the proposed candidate vaccine needs to be tested both in vivo and in vitro to determine its efficacy and safety.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.