Computational protocol: Multi-omics Analysis Sheds Light on the Evolution and the Intracellular Lifestyle Strategies of Spotted Fever Group Rickettsia spp.

Similar protocols

Protocol publication

[…] Herein, we studied four rickettsiae, including the Dermacentor-transmitted species R. slovaca strain 13-BT (CSUR R154, Fournier et al., ) and R. raoultii strain KhabarovskT (CSUR R3, ATCC VR-1596, El Karkouri et al., ) that cause SENLAT, and the Rhipicephalus-associated R. conorii strain Malish 7T (CSUR R41, ATCC VR-613, Ogata et al., ) and the R. massiliae strain MTU5 (CSUR R132, Blanc et al., ) causing MSF. For each disease, the former species caused a more severe infection. The four species were obtained from the French “Collection de Souches de l'Unité des Rickettsies” (CSUR).Genomic sequences of the four species were downloaded from the NCBI FTP server ( To avoid potential biases across the originally published data, including unpredicted Open Reading Frames from pseudogenes in the GenBank database that were generated by different gene identification and annotation systems, all genomes were subjected to a standard re-annotation, including CDSs (coding sequences) prediction with the AMIGene software (Bocs et al., ). The automatic assignment of protein functions was performed against the RickBase (Blanc et al., ) and non-redundant NR databases using PipRick (an in-house annotation pipeline written in Perl language) and BLASTp algorithm (Altschul et al., ). The annotations were then curated and genes that were either complete or altered (split or fragment) were distinguished (Blanc et al., ). Functional classification of gene families (COG ID and Letters) was searched using COGsoft software against the Clusters of Orthologs Groups (COG) database (Kristensen et al., ). The pan-genome between the four species was constructed by subjecting predicted proteomes to a reciprocal best BLAST hit (BBH) algorithm with all-against-all search (coverage of the query length ≥60% and E-value < 10−10) using COGsoft software. Each cluster of orthologous groups of rickettsial genes from chromosomes and plasmids was named cRIGs and pRIGs, respectively. The Venn diagrams of pan-genome and pan-proteome were constructed using the Jvenn Javascript library (Bardou et al., ). To identify putative virulence factors, a BLASTp search was performed against the virulence factor database, VFDB (Chen et al., ). Multiple sequence alignments of the core genes were carried out using MAFFT software (Katoh et al., ). Phylogenetic analysis was performed with the maximum likelihood (ML) method under the JTT amino acid substitution matrix, the Nearest-Neighbor-Interchange (NNI), the gamma (Γ) distribution of parameter α to account for substitution rate heterogeneity among sites and complete deletion and the rectangular tree using MEGA software (Tamura et al., ). Moreover, a neighbor-joining (NJ) tree was constructed from a gene content distance matrix calculated according to the pan-genome data and Jaccard's dissimilarity coefficient (El Karkouri et al., ). The robustness of the nodes in both ML and NJ trees was estimated through Bootstrap (BP) analyses of 100 and 500 replicates with MEGA software and the PHYLIP package (Website:, respectively. To assess the genomic differences in predicted proteins of the core genes between the four species, the percentage of the amino acid identities and the numbers of non-synonymous mutations and insertions/deletions (InDels) were computed using the Smith-Waterman BLASTp search algorithm. For this four-way analysis, the E-value cutoff < 10−10 was used and the false positive matches were removed. [...] Raw MSE data from each biological sample were processed using the ProteinLynx Global SERVER v3.0.1 (PLGS, Waters, Guyancourt, France) for protein identification and protein quantification. Noise reduction thresholds for low energy scan ion, high energy scan ion and peptide intensity were set at 1,000, 100, and 800 counts, respectively. A mass correction was applied to all spectra using the Leucine Enkephalin lock mass calibrant at 785.8426 m/z. For each sample, a single mass spectrum file was generated by merging the mass spectra from the seven fractions. To perform the protein sequence database search, we constructed four distinct databases, each containing annotated protein sequences of one Rickettsia species. For each database, we also included protein sequences of the X. laevis species downloaded from the NCBI database, the common contaminants and the yeast ADH sequence (accession number: P00330|ADH1) from the universal protein Knowledgebase, UniProt (UniProt, ). The default parameters used for global peptide and protein identifications were set as follows (see Gopinath et al., ): at least one fragment ion match per peptide, at least three fragment ion matches per protein, at least one peptide matches for protein identification, mass tolerances was set to automatic with a window of 10 ppm for precursor ions and a window of 20 ppm for fragment ions, at least one positive charge, oxidation of methionine (M) as variable modification and carbamidomethylation (C) of cysteine as the fixed modification, and the trypsin was selected as the enzyme with up to one miss-cleavage. The initial protein false discovery rate (FDR) of the identification algorithm was set at 4% with a randomized database, leading to a peptide FDR that was typically smaller than 1% (Brioschi et al., ; Gopinath et al., ). Protein quantities were evaluated in the injected solution using the combined intensity of the three most abundant peptides per protein compared to the quantitatively added yeast ADH digest (Hi3 absolute quantification, PLGS, Waters, Guyancourt, France; Silva et al., ).Only proteins identified by at least two matched peptides were considered for proteomic analysis (Treitz et al., ). As the three replicates of each sample showed high reproducibility with high significant linear correlations (Spearman's and Pearson's correlations, P < 0.001), we grouped each three replicates and calculated the average abundances (fmol or fmol μg−1) of each protein. The fmol averages of each sample were also normalized by the median using the Perseus Software (v1.5.1.6;, Treitz et al., ). Moreover, using each protein's molecular weight from the database, the fmol quantity of each protein was converted to nanograms for the triplicates, and by summing the average ng of all proteins we obtained the total average ng in each sample (Saka et al., ). To calculate the species protein abundance for each sample, the fmol average for each protein was divided by the average ng sum for only the sample of that species, and scaled by 1,000 to yield fmol μg−1 (Saka et al., ). The ratios describing the condition virulent/milder species of two comparative models (R. slovaca/R. raoultii and R. conorii/R. massiliae) were then log2 transformed. Orthologous proteins with stringently defined fold change (FC ≥ 2) were considered to represent up-regulation (log2 ratio ≥ 1) or down-regulation (log2 ratio ≤ −1), whereas orthologous proteins with a FC < 2 were considered as equally regulated (Son et al., ). Proteins with FC = 2 included a minimal and a well-represented abundance value of 0.5 fmol. Thus, orthologous proteins with stringent abundance values (≥0.5 fmol) in one species, but with zero abundance values in the second species, were also comprehensively annotated as up-regulated or down-regulated. In contrast, proteins with these abundance values (≥0.5 fmol) in only one species (i.e., the gene is absent in the second species) were annotated as specific or unique (Gopinath et al., ). Other proteins that did not match these criteria were annotated as un-classified. […]

Pipeline specifications

Software tools PLGS, Perseus, MaxQuant
Application MS-based untargeted proteomics
Organisms Homo sapiens, Rickettsia slovaca, Rickettsia conorii, Xenopus laevis
Diseases Infection, Familial Mediterranean Fever