Weekly reads 9/2/25

Systemic programs and structured models across cancer and aging

Feb 15, 2026

This week's reads look at how systemic biological programs, not isolated events, influence cancer progression, ageing, gene regulation, and lineage dynamics. In polymetastatic breast cancer, Insua-Rodríguez et al. found a conserved immunosuppressive macrophage niche driven by the MIF-CD74 axis that operates in the brain, lung, liver, and bone. Yang et al. present a first-principles measure of transcriptional entropy that quantifies ageing as a breakdown in gene coordination on a learned manifold. Vagiaki et al. reframe trans-eQTL mapping using LIVI, an interpretable generative model that identifies genetically influenced gene programs at the single-cell level. Patel and Kundaje's ARSENAL, a compact DNA language model trained on regulatory regions with motif-scale inductive bias, challenges the "scale is everything" paradigm. Yang et al. create tumor-homing probiotics that deliver chemo-immunotherapy locally with remarkable precision. Gao et al. explicitly model mitochondrial drift with MitoDrift, which turns mtDNA lineage tracing into a probabilistic, confidence-aware framework.

Preprints/articles that I managed to read this week

The MIF-CD74 Axis Drives a Systemic Immunosuppressive Niche in Polymetastatic Breast Cancer

Insua-Rodríguez et al. bioRxiv (2026). 10.64898/2026.01.31.701004

The paper in one sentence

A conserved program driven by cancer-secreted MIF recruits CD74+ lipid-associated macrophages across metastatic sites, creating an immunosuppressive niche that fuels systemic breast cancer colonization.

Summary

This study uses a synchronous multi-organ metastasis model combined with in vivo niche labeling and single-cell RNA sequencing to map the metastatic ecosystem in brain, lung, liver, and bone. The authors identify a universal proximal niche dominated by CD74+ lipid-associated macrophages (LA-MAMs), which are instructed by tumor-derived MIF via the CD74 receptor. These macrophages exhibit a potent immunosuppressive signature, drive T-cell exhaustion, and enable metastatic outgrowth across all organs. Disrupting the MIF-CD74 axis, genetically or pharmacologically, reduces LA-MAM accumulation, restores T-cell function, and impairs multi-organ metastasis. Clinical data from over 100 patients confirm that the MIF-CD74 axis is a hallmark of human polymetastatic disease and predicts poor survival.

Personal highlights

Conserved macrophage niche across metastatic organs: using in vivo proximity labeling and scRNA-seq, the authors reveal that bone marrow-derived CD74+ macrophages are consistently enriched in the immediate vicinity of metastases in the brain, lung, liver, and bone, highlighting a universal cellular convergence in systemic disease.
MIF-CD74 as a master paracrine regulator of immunosuppression: cancer cell-secreted MIF signals through CD74 on macrophages to orchestrate a lipid-associated, immunosuppressive phenotype (LA-MAM), establishing a coordinated signaling axis that is active in all metastatic sites.
Lipid-associated macrophages drive T-cell suppression and exhaustion: LA-MAMs exhibit strong lipid metabolism and oxidative phosphorylation signatures, directly suppress T-cell proliferation in co-culture assays, and correlate with increased exhausted T-cell subsets in vivo.
Disruption of the MIF-CD74 axis reduces systemic metastasis: genetic knockdown of MIF or pharmacological inhibition with 4-IPP significantly impairs metastatic burden across all four organs, demonstrating that targeting this axis can systemically compromise colonization.
Clinical translation and prognostic relevance: in a 100-patient cohort of breast cancer metastases, MIF and CD74 protein expression are nearly universal and stratify post-metastasis survival, confirming the human relevance of this mechanism and its potential as a therapeutic target

Why should we care?

This work shifts from an organ-specific “seed and soil” model to a systemic, conserved program that cancer cells deploy to colonize multiple tissues simultaneously. By identifying the MIF-CD74 axis as a central coordinator of immunosuppressive macrophages and T-cell dysfunction, it offers a unifying therapeutic vulnerability for polymetastatic disease, a condition that currently lacks effective treatments and is managed only palliatively.

A manifold-based measure of transcriptional entropy for quantifying aging in single cells

Yang et al. bioRxiv (2026). doi:10.64898/2026.01.24.701460

The paper in one sentence

Researchers developed an unsupervised, manifold-based metric called transcriptional entropy to quantify aging in single cells by measuring the breakdown of gene expression coordination, revealing cell-type-specific vulnerabilities and molecular mechanisms across tissues.

Summary

This study introduces a first-principles framework to quantify transcriptional entropy, a measure of intrinsic gene expression noise, from single-cell RNA-seq data. Unlike supervised methods that rely on predefined markers, this approach uses the deviation of a cell’s transcriptome from a low-dimensional manifold to capture loss of transcriptional coordination with age. Applied to multiple aging datasets (Tabula Muris Senis, SenNet, kidney, liver), the method identifies stem and progenitor cells as particularly vulnerable, correlates with chromatin-based mitotic age, and disentangles two aging mechanisms: loss of expression precision and activation of stress-response programs.

Personal highlights

First-principles quantification of transcriptional noise: transcriptional entropy is derived from deviations of single-cell expression from a learned manifold, capturing intrinsic biological variation without relying on clustering or predefined gene sets, making it broadly applicable across cell types and tissues.
Distinguishes aging mechanisms at the gene level: yhe framework separately identifies genes that become more disordered with age (loss of precision) and genes whose expression correlates with cellular entropy (stress-response activation), offering a dual-axis view of transcriptional dysregulation.
Cross-modal validation with epigenetic aging: transcriptional entropy strongly correlates with chromatin-based mitotic age estimates from EpiTrace, linking transcriptomic noise to epigenetic drift and reinforcing its relevance as a fundamental aging metric.
Reveals tissue- and zone-specific aging patterns: in kidney proximal tubules, entropy increases specifically in injury-prone segments (S2/S3), while in liver hepatocytes, it peaks in the regenerative midlobular zone (zone 2), highlighting spatially structured vulnerability.
Unsupervised and marker-free, outperforming existing tools: unlike SenMayo or variance-based methods, transcriptional entropy consistently detects age-related dysregulation across diverse tissues and cell types, including in compartments where traditional senescence scores fail.

LIVI: Mapping trans-eQTLs at single-cell resolution with interpretable deep learning

Vagiaki, D. et al. bioRxiv (2026). https://doi.org/10.64898/2026.02.04.703363

The paper in one sentence

LIVI is an interpretable deep learning framework that uses a structured variational autoencoder with linear decoders to decompose single-cell gene expression into cell-state and donor-specific components, enabling scalable and statistically powerful discovery of trans-eQTLs that act across gene networks and continuous cell states.

Summary

Mapping trans-eQTLs, or genetic variants that affect gene expression on different chromosomes, has traditionally been difficult due to the massive multiple testing burden (billions of potential variant-gene pairs), small effect sizes, and the complex, context-dependent nature of these regulatory effects. Existing single-cell eQTL methods are optimised for cis effects, use predefined discrete cell types, or suffer from statistical circularity when genotypes are included during training. LIVI (Latent Interaction Variational Inference) overcomes these limitations with a purpose-built variational autoencoder architecture. The model takes raw single-cell counts as input and decomposes expression into three interpretable components: (1) canonical cell-state factors (C) capturing shared transcriptional programs across donors, (2) cell-state-specific donor factors (D×C) representing inter-individual variation that interacts with cellular context, and (3) global donor factors (V) capturing population structure. LIVI uses a linear decoder for each latent space, which maintains interpretability by directly linking each latent factor to a sparse set of genes using learned weights. This enables LIVI to operate as a factor analysis model embedded within a deep generative framework, scaling to millions of cells while maintaining biological transparency. The donor factors (D) are learned without access to genetic data, preventing circularity in subsequent association testing. After training, these factors serve as compact quantitative phenotypes (typically 500-700 dimensions, compared to ~200,000 gene-cell-type combinations) for efficient eQTL mapping. The discovered associations can then be projected back to single cells using the interaction model and decoded to identify the specific genes and pathways affected, effectively converting statistical associations into mechanistic hypotheses with cell-state resolution. When applied to the OneK1K dataset (981 donors, over 1 million PBMCs), LIVI identified more trans-eQTLs than alternative latent variable methods, recovered signals missed by conventional single-gene testing, and revealed how polygenic risk for autoimmune diseases manifests in specific cell types and gene programs.

Personal highlights

Structured decomposition of genetic and cellular variation: LIVI explicitly disentangles gene expression into canonical cell-state factors, cell-state-specific donor effects, and global donor variation through a multi-decoder architecture. This separation, enforced by keeping cell-state factors fixed during donor factor training, eliminates a key source of non-identifiability and enables independent interpretation of cellular context versus genetic influence.
Interpretable deep learning via linear decoders with sparsity constraints: unlike black-box VAEs, LIVI employs sparse linear decoders that map each latent factor directly to a weighted set of genes. This design choice preserves the scalability of deep generative models while returning biologically transparent outputs: each D×C factor is a sparse, interpretable gene program whose activity can be tracked across cells, donors, and genetic variants.
Cell-state-continuous eQTL mapping without predefined annotations: LIVI does not require discrete cell type labels. By modeling donor effects as interactions with a learned, continuous cell-state space, the framework naturally captures genetic effects that span continuous trajectories or that manifest in subsets of cells that do not align with canonical annotations, a class of associations systematically missed by conventional pseudobulk approaches.
Statistically rigorous, non-circular association testing: donor factors are inferred without access to genetic data, and association testing is performed post hoc using linear mixed models that account for population structure. This avoids the permutation-based calibration required by methods that incorporate genotypes during training, yielding calibrated test statistics without computational overhead.
From factor-level associations to single-cell resolution effect maps: LIVI bridges donor-level statistics and single-cell biology. The interaction model (Softmax(C)A ⊙ d_y) allows estimated SNP effect sizes to be propagated back to individual cells, producing continuous, genome-inferred perturbation maps that reveal exactly where in the cell-state space a genetic variant exerts its influence—and on which genes.

Why should we care?

Trans-eQTLs represent the missing link between GWAS hits and the cellular phenotypes they influence. Yet they remain drastically underpowered in conventional analyses, buried under a mountain of statistical tests and obscured by their tendency to affect multiple genes with individually small effect sizes. LIVI offers a way out. By reframing trans-eQTL mapping as a problem of discovering genetically influenced gene programs rather than testing isolated variant–gene pairs, it collapses the search space from billions of hypotheses to thousands, while simultaneously aggregating weak signals into detectable, biologically coherent units.

ARSENAL: Short-context regulatory DNA language models with motif-discovery regularization

Patel, A., & Kundaje, A. bioRxiv (2026). https://doi.org/10.64898/2026.02.05.703637

The paper in one sentence

ARSENAL is a compact, short-context DNA language model pretrained exclusively on regulatory genomic regions with a frequency-domain Fourier regularizer that biases masked reconstructions toward motif-scale features, enabling superior zero-shot motif discovery, variant effect prediction, and transferable representations for downstream regulatory genomics tasks.

Summary

DNA language models (DNALMs) have emerged as powerful tools for learning regulatory sequence syntax from genomic data, yet most current approaches scale to millions or billions of parameters and train on whole-genome sequences, diluting their capacity for learning precise, motif-resolution regulatory features. ARSENAL takes a radically different, and surprisingly effective, contrarian stance: smaller is better, if you train on the right data with the right inductive biases. The model consists of a compact 8-layer transformer (768 embedding dimensions) trained exclusively on 350 bp windows centered on ENCODE candidate cis-Regulatory Elements (cCREs), a curated set of ~1.3 million experimentally supported regulatory regions. This targeted pretraining strategy concentrates the model’s capacity on functional non-coding sequence, avoiding dilution by the vast, information-sparse genomic background. But the core methodological innovation lies in ARSENAL’s Fourier motif-discovery regularizer. Drawing on prior work in supervised attribution priors, Patel and Kundaje adapt frequency-domain constraints to the self-supervised setting: the model is penalized when its per-base likelihood reconstructions contain excessive high-frequency (noise-like) or low-frequency (repeat-like) variation. This simple auxiliary loss biases the learned likelihood landscape to emphasize motif-scale features (6–20 bp), implicitly guiding the model toward the characteristic length scale of transcription factor binding sites without any supervised motif annotations.

Personal highlights

Regulatory-region-only pretraining concentrates model capacity: ARSENAL is pretrained exclusively on ENCODE cCREs rather than whole-genome sequence. This simple domain restriction ensures that the model’s limited capacity is allocated to learning functional regulatory syntax rather than memorizing repetitive, low-information genomic background, a form of data efficiency through biological curation.
Fourier-domain regularization induces motif-scale likelihood structure: the frequency-domain auxiliary loss penalizes reconstructions whose per-base likelihoods contain inappropriate spectral components, softly biasing the model toward learning 6–20 bp motif-scale features. This imports an inductive bias from supervised attribution priors into self-supervised learning, yielding interpretable likelihood landscapes without any supervised motif annotations.
Zero-shot variant scoring from 350 bp windows outperforms long-context models: ARSENAL achieves state-of-the-art correlation with experimentally measured dsQTL and caQTL effect sizes using only 350 bp of sequence context, substantially shorter than the 2 kb+ windows used in prior evaluations. This demonstrates that effective regulatory variant effect prediction does not require long-range context when motif-scale syntax is captured with sufficient fidelity.
Transferable embeddings improve supervised regulatory models: ARSENAL’s per-base embeddings, when substituted for one-hot encoding in ChromBPNet, yield consistent gains in chromatin accessibility prediction across five cell lines and improved counterfactual variant scoring. The self-supervised representations generalize across assay modalities and cellular contexts, providing useful inductive bias for downstream supervised tasks.
Controllable generation of cell-type-specific regulatory sequences: coupled with pretrained oracle models, ARSENAL supports objective-guided beam search to generate synthetic regulatory sequences with targeted properties, high predicted activity in HepG2, low activity in H1-hESC, or differential specificity. TF-MoDISco on the resulting sequences reveals emergent enrichment for appropriate cell-type transcription factor motifs, validating the approach for regulatory sequence design.

Why should we care?

For the past several years, the field of genomic language modeling has implicitly equated progress with scale: longer contexts, more parameters, more tokens, more FLOPs. ARSENAL suggests that this equation is, at best, incomplete. By training a compact model exclusively on regulatory regions and regularizing toward motif-scale features, the authors achieve state-of-the-art performance on zero-shot variant effect prediction and motif discovery, tasks that are supposed to require massive scale, with an 8-layer transformer and 350 bp windows.

Engineered probiotics for tumor-targeted combination chemoimmunotherapy

Yang, Z. et al. bioRxiv (2026). https://doi.org/10.64898/2026.02.04.703875

The paper in one sentence

A single strain of engineered E. coli Nissle 1917 delivers enzyme/prodrug chemotherapy, an IL-15 superagonist, and a PD-L1-blocking nanobody directly within tumors, achieving localized chemoimmunotherapy with synergistic antitumor immunity and minimal systemic toxicity.

Summary

Combining chemotherapy and immunotherapy is conceptually appealing but clinically challenging: chemotherapeutics have no tumour specificity and cause systemic toxicity, whereas immunotherapies rely on pre-existing immunity and can cause immune-related side effects. Yang and colleagues create a living bacterial platform that circumvents these limitations by confining both modalities to the tumour microenvironment. The authors begin with E. coli Nissle 1917 (EcN), a probiotic strain that selectively colonises tumours after intravenous administration. They modify it to express cytosine deaminase (CD), which converts the nontoxic prodrug 5-fluorocytosine (5-FC) into the cytotoxic chemotherapy drug 5-fluorouracil (5-FU). The enzyme is tagged with a heparin-binding peptide (pCD), which anchors it to extracellular matrix components and prevents it from leaving the tumour. However, wild-type EcN expresses a dihydropyrimidine dehydrogenase (encoded by the preTA operon) that converts 5-FU to inactive DHFU, which is the same metabolic pathway that causes 5-FU resistance in humans. Deleting the preTA operon eliminates bacterial drug catabolism, transforming marginal efficacy into effective tumour control. The immune phenotyping of the optimised enzyme/prodrug therapy reveals a double-edged effect: it activates CD8 T cells, NK cells, and NKT cells while also upregulating PD-L1 on suppressive myeloid populations and expanding activated regulatory T cells. This finding motivates further engineering: the authors co-express an IL-15 superagonist (s15), which promotes CD8 and NK proliferation while inhibiting Tregs, as well as a PD-L1-blocking nanobody (PDL1nb) from the same plasmid. The triple-engineered strain (EcNx^ΔpreTA-pCD/PDL1nb/s15) generates three payloads simultaneously. In the MC38 colorectal tumour model, a single intravenous dose followed by 5-FC administration results in complete tumour regression in a subset of animals, with no detectable body weight loss. The therapy activates dendritic cells, polarises M1 macrophages, promotes CD4 T cell proliferation, reverses exhaustion, and expands IFNγ-producing CD8 and NK cells.

Personal highlights

Tumor-selective enzyme/prodrug delivery with ECM anchoring: EcN bacteria naturally colonize tumors following intravenous injection, achieving >10⁹ CFU/g in tumor tissue with near-undetectable levels in healthy organs. The cytosine deaminase enzyme is tagged with a PlGF2-derived heparin-binding peptide that anchors it to extracellular matrix components, ensuring activated 5-FU remains localized rather than diffusing systemically.
Prevention of bacterial drug catabolism by preTA knockout: wild-type EcN expresses dihydropyrimidine dehydrogenase (encoded by preTA), the same enzyme that causes 5-FU resistance in humans, converting 5-FU to inactive DHFU. Deletion of the preTA operon eliminates this metabolic sink, increasing intratumoral 5-FU bioavailability and converting a marginal therapeutic effect into robust tumor control.
Mechanism-guided rational combination design: Their data-driven approach identifies IL-15 superagonist and PD-L1 blockade as logical partners to counteract the therapy’s immunosuppressive side effects while amplifying its immunostimulatory potential.
Single-strain co-delivery of three orthogonal payloads: the engineered bacteria simultaneously produce a prodrug-converting enzyme, a cytokine superagonist, and a checkpoint-blocking nanobody from a single stabilized plasmid. This demonstrates that living therapeutics can coordinate multi-agent combinations with precise temporal and spatial control, delivering chemotherapy, immunotherapy, and immunomodulation from a single intravenous injection.
Complete tumor regressions with no observable toxicity: in the MC38 model, the triple-engineered strain achieves complete regression in a subset of tumors following intratumoral injection and durable growth suppression following intravenous administration. No body weight loss or other signs of systemic toxicity are observed, a striking contrast to conventional 5-FU chemotherapy, which causes significant weight loss at efficacious doses.

Why should we care?

Cancer therapy is constrained by a persistent trade-off: effective treatments often lack specificity, and specific treatments are often ineffective. Chemotherapy kills tumors but damages healthy tissue; immunotherapy can produce durable responses but only in a minority of patients; combining them risks additive toxicities without guaranteed synergy. This work offers a way out by outsourcing drug delivery to a living system that does what no synthetic formulation can: actively home to tumors, sense its environment, and produce multiple therapeutic agents on-site, on-demand.

MitoDrift: Modeling mitochondrial inheritance enables high-precision single-cell lineage tracing in humans

Gao, T. et al. bioRxiv (2026). https://doi.org/10.64898/2026.02.12.705660

The paper in one sentence

MitoDrift is a probabilistic framework that models mitochondrial DNA heteroplasmy drift as a Wright-Fisher process, enabling confidence-refined lineage trees that accurately recover clonal relationships in primary human tissues without experimental barcoding.

Summary

The authors develop MitoDrift, a probabilistic framework that views mtDNA lineage tracing as an intracellular population-genetic process observed using noisy single-cell measurements. The model combines a discrete Wright-Fisher drift process along lineage edges and a binomial observation model at the leaves to compute tree likelihood using message passing on a hidden Markov tree. The parameters are learned using expectation maximisation, and the posterior clade support is estimated using Metropolis-Hastings MCMC sampling over tree topologies. Branches with low confidence are collapsed, resulting in a refined tree rich in accurate clonal relationships. MitoDrift is validated against orthogonal ground truth in two complementary settings. First, they use lentiviral barcoding (LARRY) in primary human HSCs, where exogenous barcodes provide precise clone assignments. MitoDrift achieves 75% clone recovery (Jaccard ≥ 0.5) and 77% clade precision, outperforming standard phylogenetic methods (NJ: 55%, UPGMA: 28%). Second, they compare their results to whole-genome sequencing of single colonies from eight healthy donors, with nuclear SNV-based phylogenies serving as ground truth. MitoDrift achieves superior precision-recall, recovering approximately 13% of clades with 50% precision and ~10% with 70% precision. This is a lower bound given the challenges of detecting low-VAF heteroplasmy in WGS.When MitoDrift is applied to native human haematopoiesis, it reveals age-related declines in clonal diversity with cell-type-specific patterns: myeloid, B, and erythroid compartments show significant reductions, while T cells maintain diversity, consistent with long-lived memory cells and lineage-biased output from dominant HSC clones. MitoDrift detects heritable regulatory programs in purified HSCs, including AP-1/stress response, stemness/lymphoid priming, and chromatin organization, with significant phylogenetic autocorrelation across longitudinal sampling. In aged donors, AP-1 transcriptional activity correlates with clone size, indicating a link between inflammatory programs and clonal expansion. MitoDrift resolves therapy-associated clonal remodelling in multiple myeloma that would otherwise be undetectable by copy number analysis. Post-treatment tumours in a deep responder (MM1) exhibit increased clonal diversity, indicating the eradication of dominant clones and outgrowth from a polyclonal reservoir. Within a therapy-resistant 1q-gain subclone, phylogeny-state analysis identifies CD44+ adhesion/migratory cells as the most closely related to post-treatment persisters, indicating a potential resistance program.

Personal highlights

Wright-Fisher hidden Markov tree models mtDNA drift explicitly: MitoDrift treats heteroplasmy evolution as a discrete-state Markov chain along lineage edges, with transition probabilities derived from a Wright-Fisher process parameterized by effective population size and generations. This population-genetic foundation replaces ad hoc distance metrics with a generative model that accounts for mutation loss, fixation, and stochastic drift, the core biological processes that have confounded mtDNA lineage tracing.
Confidence-based topology refinement prioritizes accurate clades: rather than accepting a single tree topology, MitoDrift samples tree space via MCMC and computes posterior support for each clade. Collapsing low-confidence branches produces a refined tree that explicitly acknowledges uncertainty while retaining well-supported groupings. This enables downstream analyses to focus on reliable lineage structures, trading fine-scale resolution for precision in a dataset-specific, tunable manner.
Orthogonal validation against lentiviral barcoding and WGS: the authors establish rigorous ground-truth benchmarks using two independent modalities: exogenous LARRY barcodes in primary HSCs (definitive clone assignments) and nuclear SNV-based phylogenies from single-colony WGS (time-resolved lineages). MitoDrift consistently outperforms existing methods across both benchmarks, demonstrating that its performance gains are not dataset-specific and that mtDNA-based lineage tracing can achieve quantitative accuracy.
Cell-type-specific clonal diversity in aging hematopoiesis: Applying MitoDrift to healthy donors reveals that age-associated declines in clonal diversity are not uniform across lineages. Myeloid, B, and erythroid compartments show marked reductions, while T cells preserve diversity, consistent with long-lived memory T cells and lineage-biased output from dominant HSC clones. This suggests that reduced clonal complexity in myeloid/erythroid output may compromise hematopoietic redundancy and increase stress susceptibility.
Phylogeny-state analysis links heritable programs to clonal expansion and therapy resistance: Integrating MitoDrift trees with multioomic cell states enables quantitative dissection of heritable versus plastic programs. In HSCs, AP-1/stress-associated regulons show significant phylogenetic signal and associate with clone size in aged donors. In multiple myeloma, phylogeny-state analysis within a therapy-resistant subclone identifies CD44+ adhesion/migratory cells as most closely related to post-treatment persisters, nominating a candidate mechanism of cell adhesion-mediated drug resistance.

Why should we care?

Lineage tracing is at the heart of developmental biology, cancer evolution, and stem cell biology, but barcodes cannot be used experimentally in humans. Mitochondrial DNA mutations provide a natural alternative, but their utility has been limited by a fundamental mismatch: the tools we use to analyse them (standard phylogenetic methods) assume that mutations are stable, binary, and inherited cleanly, whereas in reality mtDNA variants drift, disappear, and are measured noisily. MitoDrift addresses this mismatch by incorporating the appropriate biology into the model. Instead of treating heteroplasmy as a static trait, it explicitly simulates the Wright-Fisher drift process that governs mitochondrial inheritance. The end result is not only better trees, but trees with confidence estimates, allowing us to identify which branches are reliable and which are not. This transforms mtDNA lineage tracing from a qualitative, descriptive tool into a quantitative, hypothesis-testing framework.

Other papers that peeked my interest and were added to the purgatory of my “to read” pile

Thanks for reading.

Cheers,

Seb.

Sebcentrism

Discussion about this post

Ready for more?