HeLa cells (CCL-2) were cultured in Dulbecco’s modified Eagle’s medium (DMEM) with 10% FBS in a humid chamber at 37°C and 5% CO₂. After confluent growth, cells were washed twice with 1x sterile PBS, trypsinized, and counted in an automatic cell counter. After counting, the cells were plated (1.5x10⁶ cells/well) in a 6-well culture plate and grown for 24 h before infection. To minimize variations, we avoided passaging HeLa cells cultures and always used stock cultures passaged only once after obtaining from ATCC.

M. genitalium G37 (ATCC #33530) and M. pneumoniae M129 (ATCC #29343) were grown in 100 ml of SP-4 medium at 37°C in 150 cm² tissue culture flasks (Corning, NY) for three to four days or until the color of the medium changes to light orange. The adherent mycoplasmas are scraped from the culture flasks, washed with 1x sterile PBS, and suspended in PBS. The bacterial suspension was passed through a 23G needle to disperse the clumps and diluted to obtain OD_{600 =} 1.00, equivalent to 1X 10⁷ CFU/ml. This was further diluted with PBS to the required multiplicity of infection.

HeLa cells (1.5x10⁶ cells/well) in each well of 6-well culture plate were infected with M. genitalium G37 or M. pneumoniae at an MOI of 1:10 (1 cell: 10 bacteria) or uninfected (control cells) and incubated for 4 h at 37°C in a CO₂ incubator. This MOI and incubation time were based on our preliminary microscopic analysis. At this condition, only minimal or no necrosis/cell death of Mg or Mp infected HeLa cells was observed; hence we thought this condition was ideal for comparing the gene expression altered by these species. After the incubation, the cells were collected, washed with PBS, and RNA isolated using EZ-10 DNAaway RNA Mini-Preps Kit (Bio Basic, Toronto, Canada). Agilent Technologies 4200 TapeStation was used to analyze RNA integrity of the preparation, and RNA quantification was done using Nanodrop. 4ug of total RNA with >9 RIN (RNA Integrity Number) was used for RNA-seq library preparation.

RNA-sequencing libraries were prepared and sequenced at Novogene Corporation Inc. (Sacramento, CA, USA). After sample preparation, sample quality control was carried out by Nanodrop, agarose gel electrophoresis, and Agilent 2100. After the QC procedures, mRNA from eukaryotic organisms was enriched using oligo (dT) beads. First, the mRNA is fragmented randomly by adding a fragmentation buffer. Then the cDNA is synthesized by using an mRNA template and random hexamers primer, after which a custom second-strand synthesis buffer (Illumina), dNTPs, RNase H, and DNA polymerase I are added to initiate the second-strand synthesis. Second, after a series of steps, including terminal repair, A ligation, and sequencing adaptor ligation, the double-stranded cDNA library is completed through size selection and PCR enrichment. Quality control of the library consists of three steps: Quibit 2.0, Agilent 2100, and Q-PCR. Finally, qualified libraries are fed into Illumina sequencers (NovaSeq 6000) after pooling according to their effective concentration and expected data volume.

QC metrics at different timepoints throughout the analysis were applied to ensure good standards ( Data Sheet 1 ; Supplementary Figures 1, 2 ). Raw fastq sequencing reads quality was assessed using the bioinformatic program FastQC (Babraham Bioinformatics). The average quality per read was 36 (Phred score) ( Data Sheet 1 ; Supplementary Figure 1 ).

The quality of the transcriptome assembly was evaluated using the Hisat2 alignment metrics, FeatureCounts assignment metrics, and the bioinformatics program MultiQC (34). Data Sheet 1 ; Supplementary Table 1 shows the alignment statistics for all the biological samples.

Clustering of the sequencing samples was carried out employing Principal Component Analysis using the R statistical software (version 4.0.3) with DESeq2 package applying the function plotPCA() for each condition ( Data Sheet 1 ; Supplementary Figures 2B, 3B ).

Gene expression quality was measured by using the MA plot and volcano plot. To this end, the MA plot was generated using the R package DESeq2 calling the function plotMA() ( Data Sheet 1 ; Supplementary Figures 2C, 3C ), while the volcano plot was generated using the evaluate an expression in a data environment function (with()) ( Data Sheet 1 ; Supplementary Figures 2D, 3D ). The settings for the volcano plot include plotting log2FoldChange, -log10(p-value), and highlighting specific points that pass the padj<0.05.

Differentially expressed genes obtained from DESeq and Figure 1A were normalized using scaled TPM (transcripts per million) values obtained from the bioinformatic program StringTie gene abundance output. The heat map was plotted using R statistical software with the R package pheatmap. ( Data Sheet 1 ; Supplementary Figure 4 ).

Sequencing and mapping of NGS paired-end reads (150 bp) to the human transcriptome (hg38) with annotation from Ensembl release version 101 (August 2020) using the HISAT2 aligner (35). SAM files were converted into BAM files, then sorted and indexed using Samtools (36). Gene abundance output was generated using standard settings by employing the bioinformatics program StringTie (37). Read counts were computed using the alignment files and the FeatureCounts program from the subread-2.0.1 package (38) with reference annotation from Ensembl version 101. In the initial annotation, a total of 60,670 genomic features are considered. Counts normalization and differential gene expression analyses were calculated using the DESeq2 R package (39) and edgeR (40). Input files for DESeq2 are a metadata file with normal/experimental conditions and the read counts obtained from FeatureCounts, and two replicates were used for each condition. Settings for DESeq2 DE analysis include: study design=~ 0+Condition, rowSums(cpms)>0, estimateSizeFactors(dds), counts(sf, normalized=T), log2(norm), DESeq(dds, parallel=FALSE), results(dds, contrast = c(“Condition”, “Experimental”, “Normal”), format=‘DataFrame, pAdjustMethod =‘BH’, independentFiltering = FALSE, parallel=FALSE). Filtering criteria for differentially expressed genes of interest were only protein-coding genes obtained from the R package biomaRt annotation querying Ensembl database and a DESeq2 log2FoldChange greater than 1.5-fold change. The heat map was generated in R statistical software (version 4.0.3) using the package plotly.

Biotype classification and composition analyses were done using the Ensembl database (REST API) and biomaRt R package (41) on filtered differentially expressed genes. Gene Ontology (GO) was carried out by querying the bioinformatic Database for Annotation, Visualization, and Integrated Discovery (DAVID) (42). The gene list of interest was uploaded and queried on the online database using the Ensembl Gene ID identifier and selecting the list type as a gene list. The following annotation categories were included in the analysis: KEGG Pathways, GO Biological Processes, and UCSC TFBS. Settings for the functional annotation charts include a threshold count of 2, EASE score 0.1, Display Benjamini, and # of Records 1000. The results file was downloaded ( Data Sheets 2–14 ; Supplementary Files ) and parsed using R statistical software (version 4.0.3); subsequently, the data was visualized and plotted using the R packages ggplot2 and plotly for dot plots and heatmaps, respectively. Normalized Enrichment Scores (NES) were calculated using fgsea (Bioconductor) (43) package for fast pre-ranked gene set enrichment analysis (GSEA). Different gene datasets were obtained from the Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org) to query genes in the pathway ( Data Sheets 2–14 ; Supplementary Files ): Hallmark Pathways, KEGG Pathways, Molecular functions, and Biological processes. Settings for the analysis include a number of permutations = 1000, and stats = ranks. Heat maps were generated using plotly (R package) using the Normalized Enrichment Score (NES).

Network interaction for protein-coding genes was generated using the bioinformatic database STRING-DB (44) version (11.0b, October 2020 release). Default settings were used to query STRING-DB. The confidence score cutoff was set to 0.40. The number of indirect interactions was enhanced to obtain a complete pathway analysis (10-40 additional indirect interactions). Network visualization was done using the open-source software Cytoscape (version 3.8.2) (45), expression data (log2FoldChange) obtained from the study was included in the network using the Style/Node/Fill Color settings, as well as the edge weight, was set to the STRING-DB confidence score. Enrichment of gene ontology was added to the network analysis to supplement using the Reactome FI (46) database and identify potential relationships between differentially expressed genes and biological processes. The Cluster FI Network algorithm was applied to the network, and different modules/clusters were generated. The settings for the Reactome FI were set to standard settings using an FDR filter of 0.05, then top pathways from each module were identified and annotated.

To understand host genes regulated by human pathogenic mycoplasmas, gene expression analyses in HeLa cells infected, independently, with M. genitalium (Mg) and M. pneumoniae (Mp) were performed by directional PolyA+ RNA-seq as detailed in the methods section. A heat map ( Figure 1A and Data Sheet 1 ; Supplementary Figure 4 ) of the expressed protein-coding genes revealed that several genes in infected cells are differentially expressed compared to non-infected cells. Also, it showed significant qualitative and quantitative differences in gene expression between the Mg and Mp infected cells. In cells infected with Mg, a total of 1289 and 646 genes were up-and downregulated, respectively, in relation to non-infected cells. A similar trend was also seen with the HeLa cells infected with Mp but with a relatively lower number of up-regulated (742) and downregulated (548) genes ( Figure 1B ). To gain better insights on this, the number of genes altered or affected by each mycoplasma species was analyzed by an intersectional Venn diagram ( Figure 1B ). This indicated that many genes, 339 of the up-regulated and 221 of the downregulated genes, are commonly targeted by both mycoplasmas. Simultaneously, this displayed that each mycoplasma specifically targets several host genes. The pathogen Mg specifically up-and downregulated the expression of 918 and 410 genes, respectively, while Mp specifically up-and downregulated the expression of 388 and 295 genes. A breakdown of all the regulated genes in Mg and Mp infected cells is shown in Supplementary Files ; Data Sheet 11 . The gene showing the highest expression in Mg-infected cells was found to be the gene associated with chromosome 2 of humans and encoding the protein C2orf72, for which the function is currently unknown.

Conversely, the highest up-regulated gene in Mp infected cells was noticed to be the gene encoding the protein PECAM-1 or CD31, a protein associated with cell adhesion. Nonetheless, the highest up-regulated gene in both Mg and Mp infected cells appears to be the gene encoding the protein PPFIA4, a protein related to the disassembly of cell adhesion. Thus, overall, gene expression analyses of mycoplasma infected cells indicate that Mg and Mp target some common and species-specific genes in the host cells.

Further to understand the pathways regulated or affected by Mg and Mp in the host cells, we performed KEGG pathways enrichment analysis of the differentially expressed host genes. Data shown in Figure 2A demonstrates that Mg and Mp target distinct pathways in host cells, although there are some overlaps ( Data Sheet 1 ; Supplementary Table 3 ). The up-regulated genes with Mg infection show enrichment of genes predominantly for pathways related to nitrogen metabolism, vitamin digestion and absorption, mineral absorption, and pentose phosphate utilization ( Figure 2B ). Additionally, enrichment of genes for pathways related to inflammation and infection such as Staphylococcus aureus infection, inflammatory bowel disease (IBD), leishmaniasis infection, and tight junction are also noted. Interestingly, up-regulated genes upon Mp infection also show enrichment of genes for the latter pathways ( Figure 2C ). However, this encompasses a large number of other pathways related to infection and inflammation. Some notable pathways of this category are legionella infection, salmonella infection, malaria infection, NF-кB, Toll-like receptor signaling, cytokine, cytokine receptor interaction, chemokine signaling, and others ( Figure 2C ). Possibly this may suggest that Mp has more ability to alter the pathways related to the immune response in the host than Mg.

With regard to downregulated genes, both Mg and Mp have significant enrichment for Jak-STAT signaling, Hippo signaling, and cell adhesion pathways in the host cells ( Figures 3A, B ). In addition, however, the downregulated genes by Mg are highly enriched for critical signaling pathways such as Ras, Fox, PI3K-Akt, Rap, chemokine receptors, and TGF-β. This is very significant because these pathways are not affected by genes regulated by Mp infection.

We also analyzed how Mg and Mp infections affected the host biological processes by subjecting the differentially expressed genes to GO (Gene Ontology) analysis. As noticed for KEGG pathways, the GO analysis also showed distinct enrichment patterns for Mg and Mp regulated genes ( Figure 4A ). Biological processes related to transport (dehydroascorbic acid, bicarbonate), metal ions (cadmium, zinc), metabolism (glycolysis, gluconeogenesis, retinoid, cholesterol), macrophage differentiation, and signal transduction are the ones that show enrichment in Mg up-regulated genes ( Figure 4B ). On the contrary, enrichment of genes for biological processes such as toll-like receptor signaling, antigen processing and presentation, chemokine-mediated signaling, neutrophil chemotaxis, inflammatory response, and immune response are seen with Mp up-regulated genes suggesting that Mp primarily targets the immune process ( Figure 4C ). Besides, genes for metabolic process, namely carbohydrate biosynthetic process and sterol metabolic process ( Figure 4C ), are also enriched with Mp infection. Similar to that of up-regulated genes, the genes downregulated by Mg and Mp infection show enrichment for diverse biological processes ( Figure 5 ). While genes downregulated by Mg show enrichment for several biological processes, including signaling pathways, the genes downregulated by Mp are enriched for processes associated with ion transport (sodium, potassium, chloride) ( Figure 5 ). It should be noted that enrichment of genes for ion transport is noticed with up-regulated genes with Mg infection.

We further interrogated whether the specific genes targeted by each mycoplasma are controlled by the same set of transcription factors or different transcription factors. This revealed that Mg and Mp infections target genes regulated by distinct transcription factors ( Figure 6A ). Mg shows enrichment of up-regulated genes for transcription factors such as NFKAPPAB50, MAZR, SP1, and NMYC, which govern the expression of genes involved in inflammation, cell proliferation, and immune response ( Figure 6B ). By contrast, the pathogen Mp shows enrichment of transcription factors EGR3, EGR1, CREL, AP2GAMMA, and GATA3, which control gene expression in muscle and memory development, NFKB signaling, and inflammation ( Figure 6C ). A noteworthy observation is that both Mg and Mp have enrichment of genes transcription factor NFKAPPAB50 and related factors which play a significant role in regulating immune response and inflammation during infection. Enrichment of transcription factors for the downregulated genes of both Mg and Mp also shows some specificity. The transcription factors LHX3, FOXD3, BACH1, and CHX10, which affect the expression of genes needed for early development and immunity ( Figure 6D ), are highly enriched in the Mg downregulated genes. Conversely, the transcription factors FOXD3, NFKAPPAB, and HSF1 play roles in development, metabolism, and inflammation ( Figure 6E ), enriched in the Mp downregulated genes.

Network analysis was done to determine functional cross-talk between genes regulated by each mycoplasma infection ( Figure 7 ). Mg infection encompassed 67 nodes, of which 22 are associated with up-regulated or higher expression and 15 with downregulated or lower expression. Some noted regulated pathways include E-cadherin- and Wnt-signaling and protein-coding genes associated with phosphorylation (downregulation of MDGA1) ( Figure 7A ). Additionally, the Mg-regulated network comprised of pathways related to post-translational protein phosphorylation (VWA1-upregulated), Type I diabetes mellitus (PTPRRN- and HLA-DQB1-upregulated), endocrine resistance (BIK-upregulated), and regulation of lipid metabolism by PPAR-alpha (ANGPTL4-upregulated) ( Figure 7A ).

With MP infection, there were 65 nodes in which 14 were highly overexpressed, and the expression of 11 was attenuated. Some affected pathways include IL-4 and IL-3 signaling (JAK3- and I-1A-upregulated), IL-10 signaling (IL-1A-upregulated), malaria infection (PECAM1-upregulated), inflammatory bowel disease (IL-1A-upregulated), pathways in cancer (JAK3- and RUNX1T1-upregulated), and JAK-STAT signaling pathway (JAK3-upregulated) ( Figure 7B ). Moreover, the second cluster of genes depicted altered expression in the following pathways: Chemokine signaling pathway (CCL20- and CXCL10-upregulated), G-alpha signaling (CCL20- and CXCL10-upregulated), class A/1 receptors (EDN1-downregulated, as well as CCL20- and CXCL10-upregulated, and cytokine-cytokine receptor interaction (CCL20- and CXCL10-upregulated). In addition, network analysis showed increased interactions between ZIC1 and EN2 genes with no direct interaction with additional genes from the dataset ( Figure 7B ).

Interestingly, two gene networks involved in COVID-19 infection were found to be present indirectly in the study for Mg (CXCL8, TNF, IL-2, C3, IL-6) and Mp (CSF3, CSF2, TNF, MAPK8, CCL2, MAPK, STAT3, IL-2, IL-6, and IL-1B) ( Figure 7 ).

Overall, Mg infection led to more overexpression of protein-coding genes, affecting cellular processes to a larger extent than Mp infection.