Analysis of genomic alterations in cancer associated human pancreatic stellate cells

Pancreatic stellate cells (PSCs) constitute important cells of the pancreatic microenvironment and their close interaction with cancer cells is important in pancreatic cancer. It is currently not known whether PSCs accumulate genetic alterations that contribute to tumor biology. Our aim was to analyze genetic alterations in cancer associated PSCs. PSC DNA was matched to DNA isolated from pancreatic cancer patients’ blood (n = 5) and analyzed by Next-Generation Sequencing (NGS). Bioinformatic analysis was performed using the GATK software and pathogenicity prediction scores. Sanger sequencing was carried out to verify specific genetic alterations in a larger panel of PSCs (n = 50). NGS and GATK analysis identified on average 26 single nucleotide variants in PSC DNA as compared to the matched blood DNA that could be visualized with the Integrative Genomics Viewer. The absence of PDAC driver mutations (KRAS, p53, p16/INK4a, SMAD4) confirmed that PSC isolations were not contaminated with cancer cells. After filtering the variants, using different pathogenicity scores, ten genes were identified (SERPINB2, CNTNAP4, DENND4B, DPP4, FGFBP2, MIGA2, POLE, SNRNP40, TOP2B, and ZDHHC18) in single samples and confirmed by Sanger sequencing. As a proof of concept, functional analysis using control and SERPINB2 knock-out fibroblasts revealed functional effects on growth, migration, and collagen contraction. In conclusion, PSC DNA exhibit a substantial amount of single nucleotide variants that might have functional effects potentially contributing to tumor aggressiveness.

www.nature.com/scientificreports/ These activated PSCs are able to produce large amounts of collagens, e.g. type type I and III, laminin and fibronectin 5 and activation of stellate cells occurs in growth factor-, cytokine-and/or oxidative stress dependent manner 6 .
Tumor cells themselves may produce various growth factors such as transforming growth factor βs (TGF-βs), vascular endothelial growth factor (VEGF) and platelet-derived growth factor (PDGF), which are known to activate PSCs 6 . PSCs not only interact with tumor cells, but they also interact with many different cells in the microenvironment like islet cells, neuronal cells, endothelial cells and immune cells 1 . They can increase the proliferation and duct formation of endothelial cells. Furthermore, PSCs can also affect the immune system by reducing T-cell activity in tumors, increasing mast cell activity and elevating the myeloid suppressor cell migration into the tumor. All of these functions directly or indirectly influence tumor cells. The production of massive amounts of extracellular matrix proteins contribute to the generation of a hypoxic environment that can also support tumor growth 1,6 .
A detailed characterization and analysis of the stellate cells is crucial for the understanding of the mechanisms of how tumor and stroma cells interact with each other 3 .
It is currently not clear whether PSCs functions are only regulated at a transcriptional, epigenetic level, or whether there are also genetic alterations that might provide PSCs with an advantage in the tumor environment and therefore indirectly influence tumor growth. Building upon our previously published work showing that cancer-associated PSC display transcriptional differences as compared to those originating from chronic pancreatitis tissues 7 , this study aims to identify potential genetic alterations in cancer-associated human PSCs which might have an effect on PSC function in pancreatic cancer.

Material and methods
Patients. Pancreatic cancer tissue samples (n = 50) were collected from patients treated from June 2017 to December 2019 in the Department of Surgery of the Martin-Luther University Halle, Germany. Pathological diagnosis of pancreatic ductal adenocarcinoma and chronic pancreatitis in tissue sections was confirmed by haematoxylin and eosin staining. All tissues were fresh isolated in cell culture medium. 1 × 10 6 cells were stored at − 150 °C for further proceedings in 90% FBS 10% DMSO. Total DNA was extracted from all samples using the QIAamp DNA Blood Mini Kit according to the manufacturer´s instructions.
Blood samples for DNA isolation (250 µl EDTA blood) were obtained from the same patients. The study was approved by the ethical committee of the Martin-Luther University, Faculty of Medicine (2015-106 and 2019-037). All patients and persons involved in this study gave written informed consent. The study was performed according to the rules of the declaration of Helsinki.
Primary cell culture. Stellate cells were isolated from the tissues by using of the outgrowth method according to Bachem et al. 8 . Briefly, the tissue was cut directly after retrieval under sterile conditions into pieces of 3 mm and transferred to culture medium consisting of 400 ml low glucose Dulbecco's Modified Eagle's Medium (DMEM) (1000 mg/L, Sigma-Aldrich, Germany), 400 ml Ham ś F12 Medium (Sigma-Aldrich, Germany), 160 ml fetal bovine serum (FBS) (Sigma-Aldrich, Germany), 10 ml penicillin/streptomycin (10.000 Units penicillin and 10 mg streptomycin per ml in 0.9% NaCl, Thermofisher, Germany) and 10 ml amphotericin B (250 μg / ml, Sigma-Aldrich, Germany). The tissues with 5 ml medium were cultured in the incubator at 37 ºC and 5% CO 2 . Every second day the medium was changed with 5 ml fresh one. The average culture time was three to four weeks until there were enough PSCs for further experiments.
Immunostaining. Immunocytochemistry (ICC) on PSCs was performed with antibodies specific for periostin, GFAP, SERPINB2 and α-smooth muscle actin (Table 1). Briefly, 2.5 × 10 4 of PSCs were seeded on the slides and incubated overnight, followed by fixation with − 20 °C cold Acetone (Honeywell, Germany) and air-drying on the next day. Thereafter, the cells were incubated at 4 °C for 20 min in a 4:1 solution of Methanol (Sigma, Germany) and H 2 O 2 , respectively. Slides were then incubated with the specific antibodies overnight in a wet chamber and washed with (  9 . Resulting read pairs were mapped to the human genome (hg38) using the BWA aligner 10 . Mapped reads were further processed for variant calling according to the Best Practices Workflow 11,12 , suggested by the Broad Institute using Genome Analysis Toolkit (GATK, version 3.4) 11 . Somatic variant detection was done using the program Mutect2 from GATK (Version 4.0.4.0) against a panel of normals consisting of the corresponding blood samples 13 . Finally, all variants were annotated using ANNOVAR 14 .
PCR and Sanger sequencing. Amplification of the selected genes was performed with the primers presented in the Table 2 and Go Taq® polymerase according to the following conditions: hold 2 min at 95 °C, followed by 40 cycles of 30 s at 95 °C, 30 s at 60 °C and 30 s at 72 °C. Additional extension step at 72 °C for 30 s was included (Promega, Germany). The PCR was performed with a thermocycler (TRIO, Biometra). Products of PCR were separated on a 1% agarose gel and visualized with SYBR Safe (Thermo Fisher, Germany), followed by detection with ChemiDocTM Touch Gel Imaging System (Biorad, Germany).

Transfection.
A PAI-2 (SERPINB2) CRISPR/Cas9 KO plasmid (Santa Cruz, US) was used to create the knockout model (a knockout strategy was used and not a CRISPR-engineered mutation mimicking the identified premature stop codon). The day before transfection of the KO and the control plasmid, 1 × 10 5 cells per well were seeded into a 24-well plate. 1 µg of plasmid-DNA was carefully dissolved and mixed with 50 µl of Western blot. To analyse the knock down of SERPINB2 expression on the protein level, western blot was performed in addition to ICC staining. Proteins were extracted according to standard procedures from cultured cells 15 . Equal amounts of protein were loaded and separated by SDS-PAGE using a 12% SDS acrylamide gel. After transfer, PVDF Blotting Membranes (GE Healtcare Life science, Great Britain) were incubated with SerpinB2 antibody (1:500) (Thermo Fisher Scientific, Germany) and monoclonal rabbit antibody (1:10.000). α-tubulin (Thermo Fisher Scientific, Germany) was used as loading control (1:500). Bands were detected by chemiluminescence and quantified with densitometry.
Proliferation assay. Cells were seeded in 24-well plates (25 × 10 3 per well) in PSC-growth-medium as described above. After 24, 48 and 72 h of incubation at 37 °C in 5% CO2, Thiazolyl Blue Tetrazolium Bromide (dissolved in HBSS, Sigma-Aldrich, US) was added to the wells. At each time point DMSO (Merck, Germany) was added after another four hours of incubation. Absorbance was measured at 570 nm after the well plate was vigorously shaken.
Wound healing/migration assay. Cells were seeded in 12-well plates (300.000 cells per well) in PSCs growth medium, which was replaced by serum reduced medium (0,1% FCS) the next day. After medium-starving the cells overnight, a wound was scratched through the cell monolayer in each well (with a 100 µl pipette tip) the next morning. For quantification of wound closure, pictures were taken at timepoint zero, after four and after eight hours.
Contraction assay. Collagen gel contraction assay was performed using a Cell Contraction Assay kit (Cell Biolabs, US). Briefly, Collagen Gel Working Solution was prepared on ice according to the manufacturer's protocol and was added to a 2 × 10 6 cells/ml cellsuspension in a 4:1 ratio. 500 µl of the cell-collagen mixture was added per well in a 24-well plate and incubated for polymerization for one hour at 37 °C as well as 1 ml of PSCs www.nature.com/scientificreports/ growth medium. After two days of incubation the collagen gels were released from the sides of the culture dishes. Pictures were taken 1, 3 and 5 days after releasing the gels to quantify the change in size of the collagen gels.

Results
Pancreatic stellate cells (PSC) were isolated by the outgrowth method from patients operated on for pancreatic ductal adenocarcinoma or chronic pancreatitis. Immunocytochemistry for periostin, α-smooth muscle actin and immunofluorescence for GFAP as well as OilRed staining confirmed that the cells were PSCs (Fig. 1). DNA was extracted from PSC (n = 50 patients) as well as from blood samples of the same patients.

NGS analyses.
Initially, NGS analysis (n = 5 patients) showed between 2,559 to 5,816 variants per sample (PSC versus blood DNA). To select somatic point mutations, the program GATK was used and a coverage of 9 and an allele frequency of 0.1 was set as a cut-off. This analysis revealed on average 26 variants per sample (range 7-76) in the DNA of PSCs compared to the paired blood samples (suppl. Tables 1-5). All variants were confirmed manually utilizing the IGV viewer (an example is shown in Fig. 2). Based on a coverage of at least 30 and an allele frequency of at least 25%, 10 genes were chosen for further analysis: SERPINB2, CNTNAP4, DENND4B, DPP4, FGFBP2, MIGA2, POLE, SNRNP40, TOP2B and ZDHHC18 ( Table 3).

Data of the mutations of interest. The variants of greatest interest according to the MuTect2 analysis
data are presented in Table 3. Generally, GATK identified somatic point mutations. Genomic coordinates were visualized using the human genome 38 (hg38) as a reference sequence. All mutations were exonic. Mutations in SERPINB2 and POLE were predicted as stop gains. Further mutations are nonsynonymous single nucleotide variants, meaning the mutations caused a change in the amino acid sequence. In the majority of cases, the reference bases were Guanine (n = 7) or Cytosine (n = 3). The alternative bases were mostly Adenine (n = 5) and Thymine (n = 3). The coverage of the mutations (AD) had a range from 29 to 106.
Pathogenicity prediction. Next, pathogenicity prediction scores were analyzed for the selected variants (Table 4). The scores were used to select interesting and relevant mutations. For the SERPINB2 gene, there were no SIFT and Polyphen2 data. Nevertheless, a CADD Score of 35, which suggested the mutation might be pathogenic. The cutoff for the highest pathogenicity prediction call is 30. The MutationTaster had a maximum value of 1 such as in all cases except for the MIGA2 gene (0.9). The mutation in the MIGA2 gene showed a SIFT score of 0.075. The cutoff value for this score for being pathogenic is under 0.05. Referring only to this score, the mutation would not be predicted as pathogenic. But a value of 0.999 in the Polyphen2 Score is very close to the maximum of 1, that predicted this mutation to be very pathogenic. The CADD Score value was 23.6 for this mutation which   Functional analysis. The mutation of the SERPINB2 gene was investigated for functional relevance. To this end, immortalized fibroblasts were transfected with a SERPINB2 CRISPR/Cas9 KO plasmid. qPCR confirmed that the SERPINB2 gene was knocked out in the analyzed clones, but not in the empty vector control clones (Fig. 3A). Western blot analysis revealed knock out on the protein level as well (Fig. 3A, B), which was also confirmed by immunocytochemistry (Fig. 3C, D). Proliferation assays (Fig. 4A) revealed no relevant differences in growth after 24 h for KO clones versus control clones. In contrast, after 48 h and 72 h, knock-out clones displayed slightly reduced growth as compared to the control clones (p < 0.001). Similarly, scratch assays (Fig. 4B) demonstrated significantly reduced migration in SERPINB2-knock-out clones as compared to the control clones at 4 h and 8 h (an example is shown in Fig. 4D). Contraction assays (Fig. 4C) did not reveal a significant difference between knock-out and control clones at any time point (an example is shown in Fig. 4E).

Discussion
Pancreatic cancer is characterized by a generally poor prognosis with an overall 5-year survival rate of 10% 16 . With its main risk factors such as age, tobacco smoking, obesity, and low physical activity, it is projected to become the second leading cause of cancer-related deaths within this decade 17 . On the molecular level, key driver mutations have been identified. The KRAS gene is altered by activating mutations, and the CDKN2A gene by inactivating mutations or other mechanisms (both in > 90% of tumors) 1 . Furthermore, inactivating mutations in the TP53 gene (around 70%) and SMAD4 gene (close to 60%) are enriched in this aggressive tumor type. A large number of other mutations, including ARID1A, GL13, MLL3 and DNAH11, are altered less often (in around 10% of cases) 1 , but might have functional relevance in individual cases. Further, pancreatic tumors exhibit the characteristic of an extremely dense stroma constituting around 90% of the tumor mass 18 . The stroma plays an important role in the tumor biology by creating hypoxic conditions, acting as a potential mechanical barrier of tumor cell spread, but also diminishing efficacy of utilized chemotherapies. They further impact on carcinogenesis by increasing matrix stiffness in the tumor microenvironment directly activating epithelial-mesenchymal transition (EMT), tumor invasion, and metastasis of tumor cells 19 . Pancreatic stellate cells constitute the majority of cancer-associated fibroblasts (CAFs) in pancreatic cancer and play a crucial role in inducing this desmoplastic reaction trough production of huge amounts of ECM proteins 20 . Their role in the microenvironment of pancreatic cancer is complex with both tumor-promoting and tumor-restraining functions being described and observed in vitro and in vivo.
It is generally assumed that genetic alterations drive neoplastic transformation in primary tumor cells and that stromal cells react on these altered cells. It is not known whether PSCs acquire genetic alterations with functional relevance during pancreatic carcinogenesis. In the present study, we therefore analyzed cancer derived PSC versus matched blood genome in an exploratory cohort (n = 5) and confirmed findings in a larger cohort (n = 50). The absence of key pancreatic cancer mutations KRAS, P53, P16 and SMAD4 ruled out contamination with tumor cells. Importantly it also reveals that PSCs do not harbor these cardinal mutations associated with pancreatic cancer cells. Ten potential pathogenic mutations (in the SERPINB2-, CNTNAP4-, DENND4B-, DPP4-, FGFBP2-, MIGA2-, POLE-, SNRNP40-, TOP2B-and ZDHHC18 genes) could be identified and confirmed in www.nature.com/scientificreports/ individual cases, yet no mutation was observed in more than one of the 50 analysed cases (including also chronic pancreatitis). It is therefore unlikely that common driver mutations exist in PSCs in pancreatic cancer. From our data we cannot answer the question whether all CAFs in one tumor share mutations, or whether there are many clones present in the tumor mass. However, our lower limit of sequence coverage (30x), a variant allele frequency (VAF) of about 25% is based on the detection of at least seven mutated alleles. So, there is rather little room in the source data that would hide other somatic mutations for the following reasons. 1) detection of less than seven mutated alleles would gradually increase the chance to make a false positive mutation call. 2) of the few variants that fail our threshold, none of them would suggest a functional contribution as judged from the prediction of their pathogenicity. Besides avoiding false-positive mutation calls, a VAF above 25% also underscores the functional contribution of the respective variant. Because our method to outgrow the CAF avoids contamination by non-CAF cells, about half of the outgrowing cells would carry this variant allele if the VAF is 25%. So, it is more likely that our cut off helps to detect major subclones of the original CAF clone in the tumour. However, it does not exclude that there are smaller minor clones, but their detection is beyond the limits of our methodology.
The functional relevance of the observed mutations in PSCs in the microenvironment of pancreatic cancer is currently not known and is the focus of future studies. Interestingly, however, the detected SERPINB2 mutation is predicted to result in a stop and therefore the expression of a functional SERPINB2 protein would be lost in those PSCs in pancreatic cancer. A premature stop might have different effects as compared to a knockout, but both result in a non-functional protein.
A study by Harris et al. analyzed the role of SERPINB2 in stromal remodeling and local invasion in pancreatic cancer 21 . There, SERPINB2 -/mouse embryonic fibroblasts (MEFs) were compared to wild-type MEFs. It was demonstrated that wild-type MEFs were highly mobile, moving freely while SERPINB2 -/-MEFs were widely immobile. Phenotypically, SERPINB2 -/-MEFs had fewer but longer and more stable protrusions in contrast to wild-type MEFs that had numerous and short protrusions. Allograft experiments with nude mice (PDAC cells and MEFs (wild-type or SERPINB2 -/-) were co-injected) displayed PDAC tumors to be larger in SERPINB2 -/-MEFs than in wild-type MEFs. It was concluded that SERPINB2 plays a major role in the inhibition in PDAC local invasion trough regulation of stromal collagen remodeling. In line with these observations, we could show that SERPINB2 knock-out resulted in reduced growth and migration, yet not collagen contraction. Another study by Westrick et al. investigated that a deficiency SERPINB2 results in accelerated tumor growth 22 , which was not due to a specific deficiency within macrophages or other hematopoietic-derived cell population 22 , suggesting that SERPINB2 has a regulatory function in non-tumor cell types as endothelial cells and fibroblasts 22 .
Taken together, it is conceivable that the identified SERPINB2 mutation in cancer-associated PSCs results in enhanced pro-tumorigenic activity of PSCs.
While specific mutations (such as SERPINB2) might contribute to tumor biology in individual cases, it is also conceivable that most of the identified mutations represent passenger mutations that do not have functional effects in pancreatic cancer-associated PSCs.
Typically, alterations develop due to several risk factors such as aging, mutagenic chemicals, radiation, ultraviolet (UV) light, oxygen radicals and further factors 23 . Many alterations are defined as passengers with the ability to remain lifelong in normal tissues without having functional consequences or causing tumor development. A subdivision of genetic and epigenetic alterations composes genetic alterations as a result of aging, mutagenic chemicals, radiation and UV-light while epigenetic alterations occur next to aging through chronic inflammation, cigarette smoking and others 23 . Somatic mosaicism is a phenomenon derived from the presence of multiple cell clones containing distinct genotypes in the same individual 24 . For instance, endogenous factors like DNA-doublestrand break, inefficient DNA repair, DNA polymerase slippage and mobile elements lead to somatic mutations. Just like exogenous factors such as nicotine and UV exposure. Furthermore, there is a de novo locus-specific rate enabling genomic alterations (10 -6 to 10 -4 ) per locus and generation 25,26 . Further analysis, including single-cell analysis and comparison to other cell types would be necessary to investigate accumulations of somatic mutations and their heterogenic functions in cancer-associated human PSCs.
In conclusion, the existence of common driver like mutations in PSCs is unlikely. Nonetheless, the identified individual mutations in tumor-associated PSC genomes might affect PSCs themselves and the microenvironment of tumors. Further studies with larger number of patients, additional control tissues/cells (e.g. healthy pancreatic tissue), as well as functional studies of individual mutations are necessary. An in-depth knowledge of individual alterations -not only of cancer cells but also of non-neoplastic tumor associated cells-provides new options for investigations to generate better therapies against pancreatic cancer in the future.


Spatial Transcriptomic Analysis Reveals Associations between Genes and Cellular Topology in Breast and Prostate Cancers

Simple Summary Tissues consist of various cell types in complex spatial arrangements. These tissues require signaling between the cells and structural proteins to maintain their organization. Spatial arrangements or the organization of cells in tissue can be readily studied via imaging analysis of dissected tissue. In cancer tissues, there is, additionally, signaling between immune cells and between immune and cancer cells. The signaling molecules as well as the structural molecules can be readily measured using molecular profiling technology. Until now, there have not been ways to have both high-resolution imaging and high-resolution molecular profiling for the same locations in tissue. Now that the technology exists to measure both molecular profile and corresponding image, we can study the relationship between molecules in the tissue and the resultant tissue architecture in order to determine functional relationships between the two in cancer. We find that image features correspond well with extracellular matrix and other mechanisms that are important in tumor structure and aggressiveness. Abstract Background: Cancer is the leading cause of death worldwide with breast and prostate cancer the most common among women and men, respectively. Gene expression and image features are independently prognostic of patient survival; but until the advent of spatial transcriptomics (ST), it was not possible to determine how gene expression of cells was tied to their spatial relationships (i.e., topology). Methods: We identify topology-associated genes (TAGs) that correlate with 700 image topological features (ITFs) in breast and prostate cancer ST samples. Genes and image topological features are independently clustered and correlated with each other. Themes among genes correlated with ITFs are investigated by functional enrichment analysis. Results: Overall, topology-associated genes (TAG) corresponding to extracellular matrix (ECM) and Collagen Type I Trimer gene ontology terms are common to both prostate and breast cancer. In breast cancer specifically, we identify the ZAG-PIP Complex as a TAG. In prostate cancer, we identify distinct TAGs that are enriched for GI dysmotility and the IgA immunoglobulin complex. We identified TAGs in every ST slide regardless of cancer type. Conclusions: These TAGs are enriched for ontology terms, illustrating the biological relevance to our image topology features and their potential utility in diagnostic and prognostic models.

Introduction
In the United States, cancer is one of the primary public health issues because of its high prevalence [1]. According to the World Health Organization (WHO), cancer is considered the most common cause of death across the globe, within which breast cancer and prostate cancer are the two most widespread cancer types affecting women and men, respectively [2]. Breast and prostate cancers are the most common sex-specific cancers. In 2020 alone, breast cancer affected more than 2.25 million individuals worldwide and claimed the lives of more than 680,000 women [2]. Prostate cancer affected more than 1.4 million men worldwide in 2020 [2].
In the biology of cancer, the tumor microenvironment (TME) consists of resident tissues, various inflammatory and stromal cells, and the ECM [3]. Tumors usually contain elevated levels of ECM markers [3]. Using the levels of ECM as a tumor indicator could aid in detecting cancer in its early stages and help in treatment as well [4]. The ECM has been identified as a mediator of drug resistance, therefore, attempting to modulate the ECM could improve the efficacy of cancer treatments [3]. Histopathology can be analyzed with measures of cell morphology and organization, like topological data analysis (TDA), in order to understand cellular interactions in the TME [5].
TDA [6,7] is a useful mathematical framework that assumes data points have an inherent relational structure. This makes it amenable to high-dimensional datasets like those used in biomedical research and has previously been used for image analysis, proteomics, and other omics techniques, but bioinformatic applications are relatively recent, so validated links from topological features to molecular function are still lacking [8]. Due to its efficiency in analyzing molecular profiles in cancer cells, advanced bioinformatics methods could be promising for cancer diagnosis and treatments within the field of oncology [9]. While most existing methods focus on analyzing the topology of gene expression data, our focus is on topology arising from the spatial layout of cells observed in an image (i.e., H&E-stained pathological specimens). Such image topological features have been found relevant to clinical outcomes [10,11].
To the best of our knowledge, this is the first publication linking topological features from cancer histology to spatial transcriptomics. Previously, imaging has been leveraged in spatial transcriptomics to infer gene expression at higher resolutions [12]. In prostate cancer, a convolutional neural network pre-trained on prostate needle biopsies was applied in order to extract image features from prostate spatial transcriptomic samples. These features were correlated with a factor analysis of gene expression, demonstrating how morphology colocalizes with certain gene sets in prostate spatial transcriptomic data [13]. Many methods use image features extracted from neural networks: for example, canonical correlation analysis to link bulk gene expression to embeddings from auto-encoders [14], and using a large amount of spatial transcriptomic samples to train and predict the spatial gene expression of several genes from only histopathological images [15].
In these analyses, we used correlation analysis to identify gene and ITF clusters in a similar fashion to gene co-expression module detection methods [16]. The association between gene expression and image topological features was determined by generating correlation matrices and visualization using heatmaps. Specifically, we applied methods to generate correlation matrices that allowed us to find correlations between gene expression and ITFs in breast and prostate cancer. Genes that correlate well with ITFs are called TAGs, i.e., topology-associated genes. Functional enrichment analysis is a useful technique for distinguishing protein and gene families related to both normal and abnormal body pathways [17]. We applied this technique to TAGs in order to find gene ontology terms that are significantly associated with ITFs from breast cancer and prostate cancer.

Datasets
We downloaded a total of six ST slides from the 10x Genomics Visium platform repository. Samples were either preserved by Formalin-Fixed Paraffin-Embedding (FFPE) or methanol fixation. They contain gene expression information, spatial information, and a corresponding tissue image. From these images, we extracted 1400 image topology features using topological data analysis techniques [7]. The slides that we used are FFPE Human Breast Cancer, Parent Human Breast Cancer, Parent Visium Human Breast Cancer, Human Prostate Cancer, Prostate Acinar Cell Carcinoma, and Visium FFPE Human Normal Prostate. There are between 2000 and 4000 ST spots in each sample. Breast cancer samples were taken from the breast tissue of women older than 18 years old and the slides each contained more than 4000 genes after filtering out low-expression genes. Prostate cancer spatial transcriptomic samples contained between 2543 and 4371 genes after filtering. Summary of correlation analysis and clustering results across all datasets and image feature types can be found in Table S1.

Statistical and Correlation Analyses
R software was used for all statistical analyses. The correlation coefficients between image topology features and gene expression were used to identify TAGs. Of the genes with no more than 25% zero-read entries, the top 150 most variable genes in the ST data were selected for analysis. We selected 150 genes because they represent more than 50% of the total variance in the dataset, and this number was also large enough to provide statistical power in the functional enrichment software. Functional enrichment can help determine which proteins and genes have relationships with cancer [17]. We calculated the log-transform of gene expression and topology image feature values and then calculated the Pearson correlation coefficients between the ITFs and the genes. The pheatmap package [18] was used to visualize correlations between genes and ITFs.

ITF Generation
Image features were generated by a collaborating lab that specializes in computational pathology and TDA. They were extracted from 350 × 350-pixel image patches centered around the ST spots that corresponded to the same locations from which the gene expression was measured on the ST slide. The 350-pixel size was chosen to tile the image into nonoverlapping squares so that features from one ST spot did not overlap with any others. We first identified cell locations from H&E images. Next, we used TDA methods to generate topological features. The cell locations can be visualized as points on a 2-dimensional scatter plot. The diameter of these points grows until they touch one another. As the diameter grows, more and more circles touch. When this happens, lines drawn from the circle centers can create shapes like triangles, squares, and rings. This is referred to as the "birth" of the shape. Eventually, the circle diameter is so large that all points are connected, drowning out shapes. The lifespan ("persistence") of these shapes as the circles continue to grow is known as "persistent homology". By scanning the connectedness of data at many different resolutions, persistent homology is a way to identify robust patterns in complex data from a very small to a very large scale. Persistent homology can be visualized in a persistence diagram, which can then be converted to a vector for analysis [19]. We focused on 0-dimensional and 1-dimensional topological components, which are clusters and rings, respectively.
Topology also has a distinct advantage over traditional convolutional neural networks (CNN): the features extracted by CNN kernels are task-specific, whereas topologies are not. For example, a CNN that classifies humans, dogs, and cats, may not provide features that are useful for histopathological images. Topology, instead, is a generalizable mathematical framework used to describe the shape of data and can be used to study the spatial relationships of cells on a local and global/tissue-wide level. Traditional convolutional neural networks cannot be used to process an entire whole slide image, so cannot directly represent tissue-wide architecture. Finally, because topology is not task-specific it allows us to easily compare results across cancers. The complementary functions of topology and convolutions have led to them being jointly implemented for the analysis of histopathology in cancer [20].

Clustering Correlated Gene Expression and Topology to Identify Sets of TAGs
As mentioned in 2.2, we identified TAGs by first calculating the Pearson correlation coefficient (PCC) between gene expression with topology across the entire spatial transcriptomics sample. To identify clusters of tightly correlated genes and topology features, we used the R package pheatmap [18] to identify clusters for rows (gene expression) and columns (image topology features) with a baseline of 10 clusters, then iteratively found the optimal number of clusters using the "elbow-method". K-means and the "elbow method" are commonly used to find the optimal number of clusters [21]. The optimal number of clusters was visualized with pheatmap. These clusters of topology-associated genes were used for functional enrichment, providing the link between topology (i.e., tissue architecture) and the transcriptional signatures of specific biological phenomena.

Functional Enrichment Analysis
The R package gprofiler2 is a useful tool for performing functional analysis of gene sets [22]. This package provides a list of enrichment terms for each gene set and corresponding p-values. After adjusting p-values using the holm method, we used a significance level of 0.05 to select the gene sets with the lowest p-values, which are the highest significant terms.

Integrative ITF Analysis in FFPE Human Breast Cancer ST Data
Based on a subset of image topology features identified through the previous correlation and functional enrichment analyses, we evaluated the relationships of the extracellular matrix ITFs with immune signaling. First, we identified a subset of 12 ITFs and 19 genes that identify ST spots related to immune signaling. This was achieved by looking at the common expression of immune-related on the ST array (e.g., major histocompatibility complex class I (MHC-I), CD8, and T-cell receptor (TCR) genes), and pathologist-defined tissue labels in order to identify ITFs most specific to these immune regions. We then used our recently developed tool Diagnostic Evidence Gauge of Single-cells (DEGAS) [23] with single-cell RNA sequencing data (scRNA-seq) from the triple-negative breast cancer (TNBC) samples [24] to estimate T-cell enrichment across the whole ST array. We also used DEGAS to infer ITFs in the single cells from the same TNBC dataset in order to identify which subsets of cells were more associated with specific ITFs. Finally, in the ST data, the 12 ITFs from each spot were used to train XGBoost [25] and LightGBM [26] classifiers to predict the individual expression of the 19 immune-signaling genes at that spot using leaveone-out cross-validation. The accuracy of the gene expression prediction was evaluated using the Pearson correlation coefficient between the predicted and true gene expression in the hold-out group ( Figure S1).

1-Dimensional ITFs Correlate with Gene Sets and Functional Enrichment
We applied the workflow outlined in Figure 1. The 1-dimensional ITFs were found to have both positive and negative correlations to the 150 highest TAGs across all six of the ST slides used in this analysis ( Figure 2). Furthermore, the ITFs tended to form distinct clusters presumably due to the ITFs describing similar topologies in the image. Based on the varying levels of correlation between ITF clusters and TAG clusters in the ST arrays, these similar topology ITFs can be described in terms of the associated TAG clusters. We identified several such TAG clusters in the breast cancer ST slides.   From the Parent Visium Human Breast Cancer ST slide, we found that ITF clusters 1 and 2 were more correlated with TAG cluster 6 ( Figure 2A, Table S2). ITF cluster 3 was the most similar to TAG cluster 2, representing the enrichment of the Mitochondrial inheritance, and ITF clusters 4 and 5 were the most similar to TAG cluster 4, representing the enrichment of the Collagen type I trimer (Figure 2A, Table S2). Lastly, ITF clusters 6 and 7 were the most similar to TAG cluster 5, representing the enrichment of the ZAG-PIP complex (Figure 2A, Table S2). In the Parent Human Breast Cancer slide, we found that ITF cluster 1 was the most similar to TAG cluster 3, representing the enrichment of the Dystroglycan binding ( Figure 2B, Table S3). ITF cluster 2 was the most similar to TAG cluster 5, and ITF clusters 4 and 5 were more correlated with TAG cluster 2, representing the enrichment of Diabetic cardiomyopathy ( Figure 2B, Table S3). Lastly, ITF cluster 5 was the most similar to TAG cluster 4, representing the enrichment of the ZAG-PIP complex ( Figure 2B, Table S3). The FFPE Human Breast Cancer slide revealed that all of the ITFs were the most similar to TAG cluster 2 ( Figure 2C, Table S4). It is also worth noting that the ECM TAG clusters were found in every one of the breast cancer ST slides and though it was not the highest ranked in terms of the PCC there were ITF clusters both positively and negatively correlated with it ( Figure 2A,C). In addition to the breast cancer ST slides, there were also clear associations with TAG clusters in the prostate cancer slides.
In the Human Prostate Cancer slide, ITF clusters 1 and 2 were the most similar clusters to TAG cluster 1, representing the enrichment of the Extracellular exosome ( Figure 2D, Table S5). ITF clusters 3 and 4 were more correlated with TAG cluster 3, representing the enrichment of the Smooth Muscle Contraction ( Figure 2D, Table S5). Likewise, in the Prostate Acinar Carcinoma slide, ITF cluster 3 was the most similar to TAG cluster 1, representing the enrichment of the Extracellular exosome, and ITF clusters 1, 4, and 5 were the most similar to TAG cluster 4, representing the enrichment of the Smooth Muscle Contraction ( Figure 2E, Table S6). ITF cluster 2 was the most similar to TAG cluster 5, representing the enrichment of the Insulin-like growth factor binding ( Figure 2E, Table S6). In the normal prostate slide, Visium FFPE Human Prostate, and ITF clusters 1 and 4 were the most similar clusters to TAG cluster 4, representing the enrichment of the Smooth Muscle Contraction ( Figure 2F, Table S7). ITF cluster 2 was most correlated with TAG cluster 2, representing the enrichment of the Collagen-containing extracellular matrix, and ITF cluster 3 was most correlated with TAG cluster 1, representing the enrichment of the Extracellular exosome ( Figure 2F, Table S7).
We noted that the most significant association among all prostate ST slides was that the 1-dimensional image topology features (ITF81-ITF231) were related to the extracellular matrix. The ECM in the Human Prostate Cancer slide was gene cluster 1, which was the most similar to ITF cluster 1 (ITF1-ITF97 and ITF197-264) and ITF cluster 2 (ITF98-ITF196). In the Prostate Acinar Cell Carcinoma slide, the ECM elements were in TAG clusters 1 and 2, but only TAG cluster 1 had a high similarity with ITF cluster 3 (ITF98-ITF159). The ECM ontology terms were found in the Visium FFPE Human Normal Prostate slide TAG clusters 1, 2, and 3 as well. ITF cluster 3 was the most similar to TAG cluster 1 and ITF cluster 2 was the most similar to TAG cluster 2. Clearly, there was much overlap in 1-dimensional ITFs as they relate to the ECM.

0-Dimensional ITFs Correlate with Gene Sets and Functional Enrichment
In the prostate ST slides, the 0-dimensional ITFs did not associate nearly as much with the ECM as the 1-dimensional ITFs and, instead, regularly associated with smooth musclerelated TAG clusters. The 0-dimensional ITFs showed many similarities when compared to the 1-dimensional ITFs in the breast cancer ST slides. For instance, in the breast tissue ST slides, the ZAG-PIP complex also appeared to be associated with the 0-dimensional ITF clusters.
From the Parent Visium Human Breast Cancer slide, among the most significant correlations, we found that ITF cluster 1 was the most similar to TAG cluster 3, representing the enrichment of Mitochondrial inheritance ( Figure 3A, Table S2). ITF clusters 2, 3, and 4 were most associated with TAG cluster 5, representing the enrichment of the ZAG-PIP complex ( Figure 3A, Table S2). The Parent Human Breast Cancer slide identified ITF cluster 1 as the most similar cluster to TAG cluster 3, representing the enrichment of the Dystroglycan binding. In contrast, ITF clusters 2, 3, 4, and 5 were the most correlated with TAG cluster 5, representing the ZAG-PIP complex ( Figure 3B, Table S3). In the FFPE Human Breast Cancer slide, all of the ITF clusters were the most similar to TAG cluster 2, representing the enrichment of the Centrocecal scotoma ( Figure 3C, Table S4). Unlike the similarity between 1-dimensional and 0-dimensional ITF functional assessment, the prostate cancer 0-dimensional ITFs reflected different functions than their 1-dimensional counterparts.

Integrative Analysis Reveals Immune Signaling in the FFPE Human Breast Cancer Slide
Due to the high prevalence of ECM-related terms found in the preceding analyses and the lack of top rank mappings between the ITFs and the ECM TAG clusters in breast cancer ST slides, we next evaluated the subset of 1-dimensional ITFs in the FFPE Human Breast Cancer ST dataset with respect to immune signaling. A detailed look at the accom-  In the Human Prostate Cancer slide, all four ITF clusters were the most similar to TAG cluster 3, representing the enrichment of the Smooth Muscle Contraction ( Figure 3D, Table S5). In the Prostate Acinar Cell Carcinoma slide, ITF clusters 1, 2, 3, and 4 were the most similar to TAG cluster 5, representing the enrichment of the Smooth Muscle Contraction ( Figure 3E, Table S6). ITF cluster 5 was most correlated with TAG cluster 4, representing the enrichment of the IgA immunoglobulin complex ( Figure 3E, Table S6). Like the Human Prostate Cancer slide, the Visium FFPE Human Normal Prostate slide was readily associated with the smooth muscle because all of the ITF clusters were the most similar to TAG cluster 5 which represented Smooth Muscle Contraction ( Figure 3F, Table S7). In prostate ST tissues, the 1-dimensional ITFs appeared to be more informative of biology, evaluated by the number of unique TAG ontology terms, than the 0-dimensional ITFs.

Integrative Analysis Reveals Immune Signaling in the FFPE Human Breast Cancer Slide
Due to the high prevalence of ECM-related terms found in the preceding analyses and the lack of top rank mappings between the ITFs and the ECM TAG clusters in breast cancer ST slides, we next evaluated the subset of 1-dimensional ITFs in the FFPE Human Breast Cancer ST dataset with respect to immune signaling. A detailed look at the accompanying image to the ST slide shows that some locations in the tumor have a denser immune infiltrate than others ( Figure 4A). One mechanism by which tumor cells can be detected and attacked by the immune system is through antigen presentation on MHC-I, which is subsequently detected by TCR on CD8 T-cells ( Figure 4B). From this ST data, we estimated the relative abundance of T-cells using our own tool, DEGAS [23], and scRNA-seq from TNBC samples [24] (Figure 4C). We identified these T-cells as containing many CD8 T-cells by the abundance of CD8 in those regions ( Figure 4D). The adjacent tumor region expressed high levels of MHC-I ( Figure 4E) and the identified T-cells expressed higher levels of TCR genes compared to other lymphocyte-rich areas ( Figure 4F). These gene expression signals seem to indicate a large number of cytotoxic CD8 T-cells attacking that specific tumor region denoted by an arrow in Figure 4. We see that the topology of the cells that are measured using ITF121 is also unique in this region, possibly denoting some arrangement of cells that influences T-cell response to the HLA-I antigen presentation ( Figure 4G). We also found that when we predicted the value of ITF121 in the TNBC scRNA-seq [24] using DEGAS [23], epithelial cells expressing higher levels of MHC-I genes tended to have a higher association with ITF121 (p-value = 3.3 × 10 −5 ). From these analyses, it is clear that cellular topology, measured by ITFs, is associated with different cell types and immune signaling in the TME and the ECM. The associations between tumor-infiltrating lymphocytes (TIL) and prognosis are complicated. Tumor-infiltrating lymphocytes in Ductal Carcinoma have been studied previously. High TIL densities were associated with a higher nuclear grade, more necrosis, and a higher rate of tumor recurrence. These same authors note an opposite direction of effect in triple-negative breast cancer (the most aggressive subtype) where high levels of TILs were associated with a better prognosis [27]. The topology could be applied to identify complex patterns of tumor-infiltrating lymphocytes in cancer. This is a field of active research [5,28,29].
As shown in Figure 5, we generated a correlation matrix, clustered it on the paired ITFs and TAGs using pheatmap ( Figure 5A), then predicted the TAG expression from the ITFs by training both an AdaBoost and an XGBoost regression model (Table S8). The mean AdaBoost correlation (0.524) was 0.014 greater than the XGBoost (0.510) during leave-oneout cross-validation but was not significant (p-value = 0.07, Table S8). The ITF to TAG correlation matrix shows multiple immune-related genes ( Figure 5A). Some of these genes such as HLA-E and HLA-C are directly related to MHC-I signaling by the cancer cell. The gene STAT1 is related to T-cells and immune signaling. In Figure 5B, we further performed leave-one-out cross-validation using LightGBM in order to predict the expression levels of each of the 19 genes based on the 12 image features. A LightGBM model was trained on all but one of the spots in the ST slide, reserving the last remaining spot for testing. This was repeated for each of the 19 genes. The correlations between the predicted gene  Figure 5B). These results on the 10X Genomics Visium FFPE Human Breast Cancer slide showed that ITFs represent important cellular immune-tumor architectures that correspond to actual distinct molecular events and cell types and can be used to identify cell-cell interactions.
some arrangement of cells that influences T-cell response to the HLA-I antigen presentation ( Figure 4G). We also found that when we predicted the value of ITF121 in the TNBC scRNA-seq [24] using DEGAS [23], epithelial cells expressing higher levels of MHC-I genes tended to have a higher association with ITF121 (p-value = 3.3 × 10 −5 ). From these analyses, it is clear that cellular topology, measured by ITFs, is associated with different cell types and immune signaling in the TME and the ECM. The associations between tumor-infiltrating lymphocytes (TIL) and prognosis are complicated. Tumor-infiltrating lymphocytes in Ductal Carcinoma have been studied previously. High TIL densities were associated with a higher nuclear grade, more necrosis, and a higher rate of tumor recurrence. These same authors note an opposite direction of effect in triple-negative breast cancer (the most aggressive subtype) where high levels of TILs were associated with a better prognosis [27]. The topology could be applied to identify complex patterns of tumor-infiltrating lymphocytes in cancer. This is a field of active research [5,28,29]. As shown in Figure 5, we generated a correlation matrix, clustered it on the paired ITFs and TAGs using pheatmap ( Figure 5A), then predicted the TAG expression from the ITFs by training both an AdaBoost and an XGBoost regression model (Table S8). The mean AdaBoost correlation (0.524) was 0.014 greater than the XGBoost (0.510) during leave-one-  Table S8). The ITF to TA correlation matrix shows multiple immune-related genes ( Figure 5A). Some of these gene such as HLA-E and HLA-C are directly related to MHC-I signaling by the cancer cell. Th gene STAT1 is related to T-cells and immune signaling. In Figure 5B, we further pe formed leave-one-out cross-validation using LightGBM in order to predict the expressio levels of each of the 19 genes based on the 12 image features. A LightGBM model wa trained on all but one of the spots in the ST slide, reserving the last remaining spot fo testing. This was repeated for each of the 19 genes. The correlations between the predicte gene expression values and the true expression values across the 19 genes ranged from 0.49 to 0.58 ( Figure 5B). These results on the 10X Genomics Visium FFPE Human Brea Cancer slide showed that ITFs represent important cellular immune-tumor architecture that correspond to actual distinct molecular events and cell types and can be used to iden tify cell-cell interactions.

Discussion
All ST slides had sets of genes significantly correlated with topology image feature These gene sets were enriched for several ontology terms. Five ECM-related gene term are common to the samples from both cancers: extracellular space, extracellular exosome, e tracellular vesicle, extracellular organelle, and extracellular membrane-bounded organelle. Thes

Discussion
All ST slides had sets of genes significantly correlated with topology image features. These gene sets were enriched for several ontology terms. Five ECM-related gene terms are common to the samples from both cancers: extracellular space, extracellular exosome, extracellular vesicle, extracellular organelle, and extracellular membrane-bounded organelle. These terms had the most significant p-values of the functional enrichments.
The ECM plays a significant role in tissue and organ organization and performs various functions such as cell homeostasis, differentiation, and many other physical and biochemical functions [30]. Cancer is driven by the unregulated growth of clonal cells which can evade immune system response [31]. The ECM has been shown to enhance cancer growth by changing cell functions, which can make tissues more fertile environments for cancer cells to proliferate [32]. The ECM plays a role in the progression of cancer, where the composition, toughness, hydration, and adherence to the ECM components of cells in the tumor and microenvironment evolve during breast cancer development [33]. Extracellular vesicles and extracellular exosomes can mediate cell communication and have been shown to contribute to the spread of prostate cancers [34]. Collagen plays an important part in cancer progression, sometimes contributing to the development of cancer cells and the enhancement of tumor growth [35]. Collagen Type I Trimer was significantly enriched in our analysis. Cancer-associated fibroblasts produce collagen type I on which cancer cells can migrate using fibronectin, and drugs targeting the ECM are an evolving therapeutic avenue [35].
When the FFPE Human Breast Cancer slide was studied in more detail for immunerelated signaling that could impact the TME, we identified image features that appear to be predictive of immune-related genes (e.g., MHC-I). Furthermore, the localization of some image features appears to closely mirror immune signaling molecules: CD8, TCR, and MHC-I expression have a very similar distribution to the 1-dimensional ITF121. It is worth noting that ITF121 was also commonly associated with ECM in the prostate ST slides, along with ITF81-ITF231. In the Prostate Acinar Cell Carcinoma sample, we found enrichment for the IgA immunoglobulin complex. The prostate excretes IgA, and changes in IgA expression could be indicative of abnormal prostate immunity [36]. Immunoglobulins can appear at high levels in some cancers [37] and have been shown to be markers for an immunosuppressive subpopulation of B-cells in prostate cancer [38].
We identified the enrichment for the ZAG-PIP complex in the topology-associated genes from the Parent Visium Human Breast Cancer and Parent Human Breast Cancer ST slides. This is a complex of Zinc alpha2-glycoprotein (ZAG) with Prolactin inducible protein (PIP). ZAG has been studied more in prostate cancer, where this secretory protein has androgen-response enhancer elements and pushes the G1/S cell cycle transition [39]. PIP is also inducible by androgens, plays a role in the adaptive immune response, and is a marker that can be used to identify metastatic foci as having originated from breast tissue [40]. In Pip deficient mice, researchers found abnormal lymphoid proliferations in prostatic tissues and hypertrophy of the thymus medulla [41].
The chief advantage of this study is the novel use of ITFs to describe the cellular topologies associated with gene expression using spatial transcriptomics. Unlike convolutional neural networks, topological features can explicitly describe tissue architecture. We use the term tissue "architecture" because topology describes the shape of networks of cells through the study of "persistent homology" (described previously), providing a complement to measuring the morphology of individual cells which our group has studied extensively previously [42,43]. The chief limitation of our study is in not having visualized topological features for interpretation by a pathologist. This is of great interest to us, as we plan to study it using the deep-transfer-learning framework DEGAS [23] to identify high-risk topologies and visualize them for histopathological interpretation.

Conclusions
In this study, we took advantage of ST in order to determine how gene expression is related to cellular topology in breast and prostate cancers. We have found TAGs that are common in both cancers and some that are unique to each cancer. ECM-related enrichment terms are the most common terms in both the breast cancer and prostate cancer samples. In the prostate ST slides, 1-dimensional Image Topological Features ITF81-ITF231 were the most tightly correlated with ECM gene clusters. Among the breast cancer slides, ZAG-PIP Complex gene expression was shared. In the prostate cancer slides, enrichment for Gastrointestinal Dysmotility gene expression and IgA immunoglobin gene sets were shared across samples. The significance of this work is in demonstrating the utility of topology features in order to describe tissue architecture, and that ITFs can be correlated with gene expression. Gene sets that are correlated with ITFs are significantly enriched for biological functions. We do not contend that these biological functions are guaranteed to be present in the tissue. Instead, we introduced the important benefit of topology as means to rigorously quantify tissue architecture. Pathological assessment is usually largely qualitative (or, at best, semi-quantitative), therefore, quantification of tissue architecture is necessary in order to take advantage of histopathological information for integrative omics research. The histopathology of cancer is crucially important for the diagnosis and prognosis of patients-the amount and configuration of reactive tumor stroma, the spatial arrangement of TILs, and the density of tumor cells all have prognostic implications for patients and underlying genetic correlates. Here, we have demonstrated that topology is a promising technique to quantify the complex, high-dimensional, and highly variable characteristics of histopathological specimens, that can be effectively integrated with other biomedical data.
We can conclude that cellular topology has clear associations with the gene expression profiles from breast and prostate cancers. Understanding the biological meaning behind image features describing cellular topology can lead to improved image feature extraction methods, predictive models based on histopathological images, and integrative omics and image analysis, which may better describe the relationships between molecular events and the spatial arrangements of cells in the tumor and microenvironment.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/cancers14194856/s1, Figure S1: Pearson correlation between XGBoost-predicted gene expression based on image topology features and measured gene expression, Table S1: Summary of correlation analysis and clustering results across all datasets and image feature types,   Data Availability Statement: All of the spatial transcriptomics data and related images are freely available on the 10X Genomics data commons.


The Ton Motor

The Ton complex is a molecular motor at the inner membrane of Gram-negative bacteria that uses a proton gradient to apply forces on outer membrane (OM) proteins to permit active transport of nutrients into the periplasmic space. Recently, the structure of the ExbB–ExbD subcomplex was determined in several bacterial species, but the complete structure and stoichiometry of TonB have yet to be determined. The C-terminal end of TonB is known to cross the periplasm and interact with TonB-dependent outer membrane transport proteins with high affinity. Yet despite having significant knowledge of these transport proteins, it is not clear how the Ton motor opens a pathway across the outer membrane for nutrient import. Additionally, the mechanism by which energy is harnessed from the inner membrane subcomplex and transduced to the outer membrane via TonB is not well understood. In this review, we will discuss the gaps in the knowledge about the complete structure of the Ton motor complex and the relationship between ion flow used to generate mechanical work at the outer membrane and the nutrient transport process.

INTRODUCTION
The bacterial outer membrane (OM) is the first line of defense for Gram-negative bacteria against its environment. The outer membrane forms a resistant barrier against toxins and environmental threats yet must allow a variety of substances to cross without compromising the membrane. The OM can be permeated through outer membrane proteins that facilitate the diffusion of small molecules and nutrients into the periplasm (Nikaido, 2003;Vergalli et al., 2019). The OM lacks a hydrolysable energy source or an electrochemical gradient; thus, some nutrients are large or at low concentrations are actively transported across the OM with the aid of the Ton complex. Ton is a multi-subunit membrane protein complex (TonB-ExbB-ExbD) that uses the proton gradient across the inner membrane as its energy source (Figure 1). At the inner membrane, ExbB and ExbD harness the proton motive force (pmf) and transfer it to the TonB subunit. TonB is anchored to the inner membrane by a single N-terminal transmembrane helix and has an ordered C-terminal domain linked by a central proline-rich periplasmic domain (Kohler et al., 2010). This linker is long enough to span the whole periplasmic space, allowing the TonB C-terminal domain to reach the OM and interact with TonB-dependent transporters (TBDTs), providing energy for nutrient import. In this review, we will discuss the gaps in the knowledge about the structure of the complete Ton motor complex (TonB-ExbB-ExbD, or TBD), as well as the molecular interactions and mechanism by which the complex uses the pmf at the inner membrane to generate mechanical work at the outer membrane for the transport of nutrients. April 2022 | Volume 13 | Article 852955 Ratliff et al. The Ton Motor

OVERVIEW OF THE TON UPTAKE SYSTEM
TonB-dependent transporters are comprised of a 22-stranded ß-barrel C-terminal domain with a N-terminal plug domain inserted into the interior of the barrel (Noinaj et al., 2010). TBDTs import a range of small molecules and nutrients based on the specificity of the transporter, which include ironsiderophores, divalent metals, carbohydrates, cobalamin, and peptides (Schauer et al., 2008;Noinaj et al., 2010;Calmettes et al., 2015;Madej et al., 2020). Upon ligand binding on the extracellular face of the TBDT, conformational changes are induced, exposing the TonB box, a short, conserved N-terminal sequence, to the periplasm. The C-terminal periplasmic domain of the TonB subunit then interacts with the TBDT TonB box, forming a stable complex that physically connects the TBDT to the inner membrane with measured affinity up to tens of nanomolar (Freed et al., 2013;Sarver et al., 2018;Josts et al., 2019). What happens after the formation of the TBDT and TonB complex is widely unknown, but it is hypothesized that the energy from the pmf is transmitted to TonB and used to alter the conformation of the TBDT, eventually opening a channel in the TBDT.

STRUCTURES OF THE TON COMPLEX
Over the past decade the stoichiometry of ExbB and ExbD has been highly disputed, but recent reports strongly support that the ExbBD subcomplex has a 5:2 ratio with ExbB forming a pentameric hydrophobic central pore encircling a dimer of ExbD single helices (Figure 1; Celia et al., 2020). The first high-resolution structure of the EcExbBD subcomplex showed EcExbB as a pentameric structure (Celia et al., 2016). Mass FIGURE 1 | A schematic representation of the Ton uptake system (left) and molecular models of the different components (right). Left panel: the TonB-dependent transporter (TBDT, green cylinder) is anchored in the outer membrane (OM). The ligand binds on the extracellular face of the TBDT and exposes a conserved domain called the TonB box to the periplasmic side. The TonB-ExbB-ExbD (TBD) complex is anchored in the inner membrane (IM) and uses the proton motive force (pmf, proton gradient across the IM, symbolized with the red arrowhead) to generate force and movement. The TBD complex is made of a pentamer of the ExbB subunit (blue, orange, grey, purple, and green) that defines a central pore in which a dimer of ExbD subunits (red and yellow) resides. The TonB subunit (gold) binds at the periphery of the ExbBD subcomplex. The elongated periplasmic domain of TonB allows its C-terminal globular domain to reach the OM and form a stable interaction with the TBDT TonB box. It is hypothesized that the ExbBD subcomplex forms the proton channel and that the energy derived from proton translocation is propagated through the TonB subunit to the TBDT, eventually opening a channel into the TBDT and allowing the bound ligand to diffuse into the periplasm. Right panel: molecular representations of known components of the Escherichia coli Ton system. The structural models are shown with ribbons and molecular surfaces. The color coding is the same than for the schematic view on the left. The TM and flexible periplasmic domain on the EcTonB subunit are not known and shown as in the schematic representation. The crystallographic structure of EcBtuB (green) in complex with the EcTonB periplasmic domain (gold) is represented in the OM (pdb 2GSK; Shultis et al., 2006). Two models of the NMR structure of the EcExbD periplasmic domains are shown in red and yellow (pdb 2PFU; Garcia-Herrero et al., 2007). The cryo-EM structure of the EcExbBD complex is shown in the IM (pdb 6TYI; Celia et al., 2019). The grey and purple ExbB subunits are not represented in order to show the ExbD TM domains. Molecular graphics have been performed with UCSF ChimeraX (Pettersen et al., 2021). spectrometry experiments performed on Escherichia coli-native membranes further supported the pentameric nature of EcExbB (Chorev et al., 2018). Cryo-EM single-particle analysis (SPA) was also used to determine high-resolution structures of the ExbBD complexes from E. coli, Pseudomonas savastanoi, and Serratia marcescens, all revealing a 5:2 ratio (Celia et al., 2019;Deme et al., 2020;Biou et al., 2021). MotAB and PomAB are the motor complexes that power the rotation of the flagellum, using the proton or sodium gradient, respectively. Both share extensive homology with ExbBD and are believed to derive from a common ancestor (Marmon, 2013;Lai et al., 2020).
MotAB, PomAB, and ExbBD share a high level of conservation in the transmembrane helices that form the central hydrophobic pore of ExbB/MotA/PomA. Recently published cryo-EM structures of MotAB and PomAB also confirmed 5:2 ratio of these motor complexes (Deme et al., 2020;Santiveri et al., 2020). For all these structures of ExbBD and MotAB determined by cryo-EM, the periplasmic domains of ExbD and MotB were not visible because of their high flexibility.
TonB appears to cycle through binding and release from TBDTs, sometimes referred to as the TonB energization cycle (Larsen et al., 1999;Kaserer et al., 2008;Gresock et al., 2015). While the structure of the EcExbBD subcomplex is now established, a detailed molecular description of a full TBD complex is not yet available. A cryo-EM map of the full P. savastanoi PSExbB complex has been reported (EMD-10897; Deme et al., 2020). The 3.8 Å resolution structure clearly shows the PsExbB and PsExbD subunits, with an additional rod-like structure that likely corresponds to a single TMH of a TonB subunit. The density traverses the micelle on the exterior of the PsExbB complex and packs against a region of TM1 of one PsExbB, in a location predicted by analysis of covariance and close to the cytoplasmic leaflet of the membrane (Deme et al., 2020). The lower density and resolution of the TonB TMH in the cryo-EM map were attributed to a partial dissociation of TonB from ExbBD upon freezing. Cryo-EM SPA of the E. coli EcTBD complex brought similar observations, showing the TonB TMH in a similar orientation compared to EcExbBD, and interacting with part of TM1 of EcExbB (H.C. unpublished results). The precise oligomeric state of TonB in the TBD complex is still not clearly known and needs to be further investigated.
Numerous studies have focused on the molecular interactions between the periplasmic domains of ExbD and TonB. The C-terminal domain of TonB can form a dimer in vivo (Sauter et al., 2003). The physiological importance of this dimer is not well-established, but it might be involved in the binding to the peptidoglycan (PG) layer, therefore localizing the TonB C-terminal domain close to the OM Postle, 2004, 2005;Kaserer et al., 2008;Postle et al., 2010). The C-terminal folded domain of EcExbD is a dimer in the EcExbBD complex (Figure 1; Gresock et al., 2015;Celia et al., 2016). The EcExbD dimerization interface might be altered during the energization cycle as some EcExbD monomers were found to be interacting with the TonB C-terminal domain in a pmf-dependent fashion Postle, 2011, 2012). This network of interactions likely reflects a dynamic interplay between ExbD and TonB during the energization process.
Both TonB and ExbD have a flexible periplasmic linker between their TM and C-terminal folded domains (Figure 1). Recent studies have highlighted the importance of the disordered, periplasmic linker domain of EcExbD. A conserved motif just upstream of the TM domain of EcExbD, V45, V47, L49, and P50 was found to be required for Ton function (Kopp and Postle, 2020). The TonB periplasmic linker is long enough to allow the C-terminal folded domain of TonB to reach the TBDTs in the OM. Most TonB sequences exhibit a conserved proline-rich domain in the linker region that is suspected to adopt a poly-proline type II helical rod conformation, conferring rigidity to the linker (Chu et al., 2007;Kohler et al., 2010). However the presence of this poly-proline-rich domain was found nonessential for energy transduction (Larsen et al., 1994).

ExbBD AND MotAB SHARE A SIMILAR PROTON CHANNEL
The Ton, Tol, and Mot complexes use the pmf to generate movement and share extensive homology (Cascales et al., 2001;Marmon, 2013;Ratliff et al., 2021). Tol is involved in the regulation of the OM integrity and cell division through interaction with the OM-associated Pal-TolB protein complex. Mot is powering the flagellum rotation. The Ton and Tol complexes are the most closely related: the TonB-ExbB-ExbD and TolA-TolQ-TolR subunits have the same topology, are highly homologous, and cross-complementation between TolA and ExbBD, and TonB and TolQR, has been observed (Braun and Herrmann, 1993;Lloubes et al., 2012). While there is no reported structure of TolQR, it is expected to be very similar to ExbBD (for more information on the Tol system, see Szczepaniak et al., 2020; and the dedicated review in this issue of Frontiers in Microbiology).
The MotAB complex is homologous to ExbBD/TolQR but lacks a TonB/TolA-like subunit. It associates with the motility apparatus and uses the pmf to generate torque, driving the rotation of the flagellum (Minamino and Imada, 2015;Lai et al., 2020). Like ExbD and TolR, MotB has a single TM domain, followed by a flexible periplasmic linker and a folded C-terminal domain. Between the TM and the periplasmic linker, MotB has a conserved sequence that associates with MotA and acts as a plug that prevents the flux of protons through the complex (Hosking et al., 2006). It is believed that upon association of the MotAB complex to the flagellar apparatus, the C-terminal domains of MotB bind the PG and the plug domains dissociate from MotA, allowing the flux of proton and torque generation (Kojima et al., 2009(Kojima et al., , 2018O'Neill et al., 2011).
Several high-resolution cryo-EM structures of MotAB have been reported, and all show the same 5:2 architecture as ExbBD (Celia et al., 2019;Deme et al., 2020;Santiveri et al., 2020;Biou et al., 2021). The transmembrane region that forms the pore encircling the ExbD/MotB TM helices is particularly conserved between MotAB and ExbBD (Deme et al., 2020;Ratliff et al., 2021). The Taylor lab reported three distinct structures of the Campylobacter jejuni CjMotAB complex: CjMotAB in the plugged state, a construct lacking the CjMotB plug domain, and a construct lacking the CjMotB plug domain plus a substitution of the conserved essential Asp22 residue in the CjMotB TM domain into Asn, so as to mimic the protonated state of Asp22 (Santiveri et al., 2020). Using Mole 2.5, a software that detects tunnels and cavities in macromolecules, a potential proton channel that connects the conserved Asp22 on CjMotB to the periplasm was revealed (Pravda et al., 2018;Santiveri et al., 2020). Close to the opening on the periplasmic side, the side chain of residue Phe186 on CjMotA TM4 was found in two different conformations, eventually acting as a gate that would open and close the proton channel to the periplasm. Sequence alignments show that the Phe186 is a conserved residue with a bulky side chain, mostly appearing as Leu in the consensus conserved sequence of MotA, ExbB, and TolQ (see Figure 2E; Cascales et al., 2001).
We used a similar approach to probe the structure of EcExbBD (PDB 6TYI) for potential channels. Using the same parameters for Mole 2.5 and omitting in the calculation the side chain of Leu178 in TM3 of EcExbB (equivalent of Phe186 in TM4 of CjMotA), a channel was found that connects the periplasmic side of EcExbBD and the essential Asp25 on EcExbD TM that opens to the cytoplasmic cavity (Figures 2A-D). The channel opens between two EcExbB subunits (blue and green represented chains) on the periplasmic side, close to Phe156 (Figures 2C,D). In the configuration shown in Figure 2C, the channel is occluded by Leu178 side chain. As for Phe186 of CjMotA, Leu178 could act as a gate that would modulate the opening of the proton channel. Figure 2E shows the consensus sequence resulting from multiple alignments of ExbB/TolQ/ MotA and ExbD/TolR/MotB (Cascales et al., 2001). Highlighted in blue and green (EcExbB TM2-3) and red (EcExbD TM) are the residues that are lining up the channel found in the EcExbBD structure.
It is noteworthy that most of these conserved residues apparently involved in the channel formation are highly conserved among ExbB, TolQ, and MotA as shown in the consensus sequence Figure 2E. The channels for EcExbBD and CjMotAB are remarkably similar, reflecting the high homology of the two systems, and likely a similar usage of proton translocation to generate movement.
Rotary models have been proposed for the MotAB complex to harness the pmf, which are discussed elsewhere in this issue of Frontiers in Microbiology. These models suggest that the MotA pentamer rotates around the MotB dimer at 36° increments for each proton translocated (Deme et al., 2020;Santiveri et al., 2020). While it is not yet clear what prompts the movement of MotA following the protonation/deprotonation of the conserved Asp on MotB, it is likely that a rotation mechanism is responsible for torque generation. Since the MotAB and ExbBD complexes are highly homologous, it is expected that the use of the pmf is the same and that the ExbB pentamer rotates around the ExbD dimer as well during proton translocation.
Based on these models, we propose a rotary model in which the ExbD TMs rotate by increments of 36° for each proton translocated ( Figure 2F). It is hypothesized that the proton travels through the open channel from the periplasm to the highly conserved Asp on ExbD TM, inducing a change of conformation of protonated Asp that would drive the power stroke, resulting in the rotation of the two ExbD TMs by 36°. The change of conformation would close the channel on the side of the protonated Asp, while a new channel would open on the opposite ExbD TM, allowing a new cycle to take place. The rotation would then continue as long as the channels are in the open state. Is it not known how the proton channel opens, and how the open state is regulated, but it likely depends on the TonB subunit.

WRAP AND PULL MODEL MECHANISM
Several models for the mechanism of action of Ton have been proposed. They all rely on the association of the TonB C-terminal domain with the TBDT TonB box, and the application of force to displace or alter the conformation of the TBDT plug domain, either through rotation or pulling of TonB (Chimento et al., 2005;Klebba, 2016).
In the pulling model, the TonB C-terminal domain bound to the TBDT is pulled into the periplasm by the ExbBD complex and gradually unfolds the TBDT plug domain, eventually opening a channel large enough to allow the bound nutrient to diffuse toward the periplasm (Chimento et al., 2005). Molecular simulations using the BtuB TBDT in complex with TonB show that the interaction between BtuB and TonB is strong enough to sustain a pulling force perpendicular to the OM plane and would partially unfold the BtuB plug domain (Gumbart et al., 2007). This mechanism is supported by in vitro single-molecule force spectroscopy experiments on the BtuB/TonB and FhuA/ TonB complexes, providing direct evidence that the interaction between TonB and the TonB box is strong and can withstand the amount of force needed to unfold half of the TBDT plug domain before dissociation occurs (Hickman et al., 2017). Some bacteriocins bind TBDTs and hijack the Ton system to gain access and kill bacteria with high efficiency (Atanaskovic and Kleanthous, 2019). The bacteriocin pyocin S2 has its own TonB box motif and has been shown to translocate through the P. aeruginosa iron transporter FpvAI (White et al., 2017). Pyocin S2 binds FpvAI at the same binding site than the natural pyoverdin-iron siderophore and initiates the partial unfolding of the plug domain through interaction of the FpvA TonB box with TonB. The TonB box motif of pyocin S2 is then presented to TonB through the channel created, and the force exerted by the Ton system eventually unfolds the pyocin and drags it into the periplasm (White et al., 2017). This work represents the first time the translocation of a bacteriocin through a TBDT was shown, and the mechanism is in good agreement with the pulling model.
To reconciliate the pulling model with an eventual rotation of ExbBD, we propose a model in which the rotation of the ExbD dimer leads to the wrapping of the TonB linker around ExbD, leading to the pulling of the TonB C-terminal domain bound to the TBDT TonB box (Figure 3). In this model, the tethering of TonB to the OM, through the interaction between the TonB C-terminal domain and the TBDT TonB box, induces a conformational change in ExbBD that opens the proton channel. The translocation of protons through ExbBD leads to the rotation of the ExbD dimer within the ExbB pentamer ( Figure 2F) and the associated TonB. The periplasmic domain  Cascales et al. (2001). Lowercase letters represent residues present in the 60% consensus and uppercase letters for residues in the 90% consensus. Gaps are marked "-" and "*" when present in the 90 and 60% consensus, respectively (Cascales et al., 2001). The numbering corresponds to the E. coli sequences of ExbB and ExbD. The highlighted residues in blue, green, and red are the ones shown with arrows on (C). They all are in the 90% consensus range, except for Thr181. The regions highlighted in grey correspond to the last two TMs of ExbB, TolQ, and MotA, and the single TM of ExbD, TolR, and MotB. (F) Rotary model of ExbBD. The view is the same as in (B) and shows a schematic slice of the TM domains of ExbB and ExbD. The positions of the proton channel are shown with the pink circle. The cycle starts with the channel opening between the green and blue ExbB subunits. The proton travels to the conserved Asp on the red ExbD TM, inducing a conformational change resulting in the rotation of the two ExbD TMs by 36°. The conformational changes lead to the closure of the channel between the green and blue ExbB subunits, while a new channel opens between the grey and orange ExbB subunits. A second proton now travels to the Asp on the yellow ExbD TM, resulting in a new rotation of 36°. The channel between the grey and orange ExbB subunits closes, while a new channel opens between the purple and green subunits, allowing a third proton to travel to the Asp on the red ExbD TM. The rotation can proceed as long as the channels are in the open state. Molecular graphics have been performed with UCSF Chimera (Pettersen et al., 2004).
of ExbD comes into contact with the TonB periplasmic linker, which starts to wrap around ExbD while the rotation proceeds. The wrapping of the two proteins leads to a pulling motion on the C-terminal domain of TonB, which gradually unfolds the plug domain of the TBDT. With the unfolding of the plug domain proceeding, a channel through the TBDT gradually opens, allowing the ligand to move into the periplasm. At some point, the force necessary to further unfold the plug domain is greater than the interaction between the TonB box and TonB, causing TonB to dissociate from the TBDT. TonB is no longer tethered to the OM, causing the proton channel to close. The TonB box and plug domains fold back into the TBDT barrel, and the TonB linker unwraps, reverting both the TBDT and TBD complex to their resting states. Further work is needed to better understand how the Ton system functions. The high-resolution structures of ExbBD and MotAB have provided insights on how the pmf is used by these protein complexes. The high structural and sequence homologies between the two systems suggest that they share the same mechanism, and any new information gathered on one system will likely translate to the other. Some of the important unanswered questions for the Ton system are: how is the signal of the binding of TonB to the TBDT TonB box at the OM transmitted to TBD at the IM? what are the molecular events that open the proton channel in ExbBD?; and how is the potential rotation of ExbB and ExbD transferred to TonB? A combination of in vivo and in vitro approaches is needed to further understand this multicompartment system.

AUTHOR CONTRIBUTIONS
AR, SB, and HC contributed to writing the manuscript and editing. HC prepared figures. All authors contributed to the article and approved the submitted version.  (5). (6) The unfolding of the plug domain reaches a point where the opening of the channel allows the ligand to diffuse to the periplasm, and the force necessary to unfold the rest of the plug is greater than the force necessary to maintain the interaction between TonB and the TonB box. TonB then dissociates from the TBDT (E), releasing the tension on the TonB periplasmic linker and reverting the ExbBD complex to the closed state. (7) The plug domain folds back into the TBDT barrel; and (F) the TBD complex goes back to its resting state.