Clustering gene expression time series data using an infinite Gaussian process mixture model

Single nucleus RNA-sequencing integrated into risk variant colocalization discovers 17 cell-type-specific abdominal obesity genes for metabolic dysfunction-associated steatotic liver disease

BACKGROUND

Abdominal obesity increases the risk for non-alcoholic fatty liver disease (NAFLD), now known as metabolic dysfunction-associated steatotic liver disease (MASLD).

METHODS

To elucidate the directional cell-type level biological mechanisms underlying the association between abdominal obesity and MASLD, we integrated adipose and liver single nucleus RNA-sequencing and bulk cis-expression quantitative trait locus (eQTL) data with the UK Biobank genome-wide association study (GWAS) data using colocalization. Then we used colocalized cis-eQTL variants as instrumental variables in Mendelian randomization (MR) analyses, followed by functional validation experiments on the target genes of the cis-eQTL variants.

FINDINGS

We identified 17 colocalized abdominal obesity GWAS variants, regulating 17 adipose cell-type marker genes. Incorporating these 17 variants into MR discovers a putative tissue-of-origin, cell-type-aware causal effect of abdominal obesity on MASLD consistently with multiple MR methods without significant evidence for pleiotropy or heterogeneity. Single cell data confirm the adipocyte-enriched mean expression of the 17 genes. Our cellular experiments across human adipogenesis identify risk variant -specific epigenetic and transcriptional mechanisms. Knocking down two of the 17 genes, PPP2R5A and SH3PXD2B, shows a marked decrease in adipocyte lipidation and significantly alters adipocyte function and adipogenesis regulator genes, including DGAT2, LPL, ADIPOQ, PPARG, and SREBF1. Furthermore, the 17 genes capture a characteristic MASLD expression signature in subcutaneous adipose tissue.

INTERPRETATION

Overall, we discover a significant cell-type level effect of abdominal obesity on MASLD and trace its biological effect to adipogenesis.

FUNDING

NIH grants R01HG010505, R01DK132775, and R01HL170604; the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant No. 802825), Academy of Finland (Grants Nos. 333021), the Finnish Foundation for Cardiovascular Research the Sigrid Jusélius Foundation and the Jane and Aatos Erkko Foundation; American Association for the Study of Liver Diseases (AASLD) Advanced Transplant Hepatology award and NIH/NIDDK (P30DK41301) Pilot and Feasibility award; NIH/NIEHS F32 award (F32ES034668); Finnish Diabetes Research Foundation, Kuopio University Hospital Project grant (EVO/VTR grants 2005-2021), the Academy of Finland grant (Contract no. 138006); Academy of Finland (Grant Nos 335443, 314383, 272376 and 266286), Sigrid Jusélius Foundation, Finnish Medical Foundation, Finnish Diabetes Research Foundation, Novo Nordisk Foundation (#NNF20OC0060547, NNF17OC0027232, NNF10OC1013354) and Government Research Funds to Helsinki University Hospital; Orion Research Foundation, Maud Kuistila Foundation, Finish Medical Foundation, and University of Helsinki.

Alternative splicing coupled to nonsense-mediated decay coordinates downregulation of non-neuronal genes in developing mouse neurons

BACKGROUND

The functional coupling between alternative pre-mRNA splicing (AS) and the mRNA quality control mechanism called nonsense-mediated decay (NMD) can modulate transcript abundance. Previous studies have identified several examples of such a regulation in developing neurons. However, the systems-level effects of AS-NMD in this context are poorly understood.

RESULTS

We developed an R package, factR2, which offers a comprehensive suite of AS-NMD analysis functions. Using this tool, we conducted a longitudinal analysis of gene expression in pluripotent stem cells undergoing induced neuronal differentiation. Our analysis uncovers hundreds of AS-NMD events with significant potential to regulate gene expression. Notably, this regulation is significantly overrepresented in specific functional groups of developmentally downregulated genes. Particularly strong association with gene downregulation is detected for alternative cassette exons stimulating NMD upon their inclusion into mature mRNA. By combining bioinformatic analyses with CRISPR/Cas9 genome editing and other experimental approaches we show that NMD-stimulating cassette exons regulated by the RNA-binding protein PTBP1 dampen the expression of their genes in developing neurons. We also provided evidence that the inclusion of NMD-stimulating cassette exons into mature mRNAs is temporally coordinated with NMD-independent gene repression mechanisms.

CONCLUSIONS

Our study provides an accessible workflow for the discovery and prioritization of AS-NMD targets. It further argues that the AS-NMD pathway plays a widespread role in developing neurons by facilitating the downregulation of functionally related non-neuronal genes.

Mediator kinase inhibition impedes transcriptional plasticity and prevents resistance to ERK/MAPK-targeted therapy in KRAS-mutant cancers

Acquired resistance remains a major challenge for therapies targeting oncogene activated pathways. KRAS is the most frequently mutated oncogene in human cancers, yet strategies targeting its downstream signaling kinases have failed to produce durable treatment responses. Here, we developed multiple models of acquired resistance to dual-mechanism ERK/MAPK inhibitors across KRAS-mutant pancreatic, colorectal, and lung cancers, and then probed the long-term events enabling survival against this class of drugs. These studies revealed that resistance emerges secondary to large-scale transcriptional adaptations that are diverse and cell line-specific. Transcriptional reprogramming extends beyond the well-established early response, and instead represents a dynamic, evolved process that is refined to attain a stably resistant phenotype. Mechanistic and translational studies reveal that resistance to dual-mechanism ERK/MAPK inhibition is broadly susceptible to manipulation of the epigenetic machinery, and that Mediator kinase, in particular, can be co-targeted at a bottleneck point to prevent diverse, cell line-specific resistance programs.

A chromatin code for limb segment identity in axolotl limb regeneration

The salamander limb correctly regenerates missing limb segments because connective tissue cells have segment-specific identities, termed "positional information". How positional information is molecularly encoded at the chromatin level has been unknown. Here, we performed genome-wide chromatin profiling in mature and regenerating axolotl limb connective tissue cells. We find segment-specific levels of histone H3K27me3 as the major positional mark, especially at limb homeoprotein gene loci but not their upstream regulators, constituting an intrinsic segment information code. During regeneration, regeneration-specific regulatory elements became active prior to the re-appearance of developmental regulatory elements. In the hand, the permissive chromatin state of the homeoprotein gene HoxA13 engages with the regeneration program bypassing the upper limb program. Comparison of regeneration regulatory elements with those found in other regenerative animals identified a core shared set of transcription factors, supporting an ancient, conserved regeneration program.

CLOCK evolved in cnidaria to synchronize internal rhythms with diel environmental cues

The circadian clock enables anticipation of the day/night cycle in animals ranging from cnidarians to mammals. Circadian rhythms are generated through a transcription-translation feedback loop (TTFL or pacemaker) with CLOCK as a conserved positive factor in animals. However, CLOCK's functional evolutionary origin and mechanism of action in basal animals are unknown. In the cnidarian Nematostella vectensis, pacemaker gene transcript levels, including NvClk (the Clock ortholog), appear arrhythmic under constant darkness, questioning the role of NvCLK. Utilizing CRISPR/Cas9, we generated a NvClk allele mutant (NvClkΔ), revealing circadian behavior loss under constant dark (DD) or light (LL), while maintaining a 24 hr rhythm under light-dark condition (LD). Transcriptomics analysis revealed distinct rhythmic genes in wild-type (WT) polypsunder LD compared to DD conditions. In LD, NvClkΔ/Δ polyps exhibited comparable numbers of rhythmic genes, but were reduced in DD. Furthermore, under LD, the NvClkΔ/Δ polyps showed alterations in temporal pacemaker gene expression, impacting their potential interactions. Additionally, differential expression of non-rhythmic genes associated with cell division and neuronal differentiation was observed. These findings revealed that a light-responsive pathway can partially compensate for circadian clock disruption, and that the Clock gene has evolved in cnidarians to synchronize rhythmic physiology and behavior with the diel rhythm of the earth's biosphere.

Proteomics of mouse brain endothelium uncovers dysregulation of vesicular transport pathways during aging

Age-related decline in brain endothelial cell (BEC) function contributes critically to neurological disease. Comprehensive atlases of the BEC transcriptome have become available, but results from proteomic profiling are lacking. To gain insights into endothelial pathways affected by aging, we developed a magnetic-activated cell sorting-based mouse BEC enrichment protocol compatible with proteomics and resolved the profiles of protein abundance changes during aging. Unsupervised cluster analysis revealed a segregation of age-related protein dynamics with biological functions, including a downregulation of vesicle-mediated transport. We found a dysregulation of key regulators of endocytosis and receptor recycling (most prominently Arf6), macropinocytosis and lysosomal degradation. In gene deletion and overexpression experiments, Arf6 affected endocytosis pathways in endothelial cells. Our approach uncovered changes not picked up by transcriptomic studies, such as accumulation of vesicle cargo and receptor ligands, including Apoe. Proteomic analysis of BECs from Apoe-deficient mice revealed a signature of accelerated aging. Our findings provide a resource for analysing BEC function during aging.

Bayesian profile regression for clustering analysis involving a longitudinal response and explanatory variables

The identification of sets of co-regulated genes that share a common function is a key question of modern genomics. Bayesian profile regression is a semi-supervised mixture modelling approach that makes use of a response to guide inference toward relevant clusterings. Previous applications of profile regression have considered univariate continuous, categorical, and count outcomes. In this work, we extend Bayesian profile regression to cases where the outcome is longitudinal (or multivariate continuous) and provide PReMiuMlongi, an updated version of PReMiuM, the R package for profile regression. We consider multivariate normal and Gaussian process regression response models and provide proof of principle applications to four simulation studies. The model is applied on budding yeast data to identify groups of genes co-regulated during the Saccharomyces cerevisiae cell cycle. We identify 4 distinct groups of genes associated with specific patterns of gene expression trajectories, along with the bound transcriptional factors, likely involved in their co-regulation process.

The transcriptional regulatory network modulating human trophoblast stem cells to extravillous trophoblast differentiation

During human pregnancy, extravillous trophoblasts play crucial roles in placental invasion into the maternal decidua and spiral artery remodeling. However, regulatory factors and their action mechanisms modulating human extravillous trophoblast specification have been unknown. By analyzing dynamic changes in transcriptome and enhancer profile during human trophoblast stem cell to extravillous trophoblast differentiation, we define stage-specific regulators, including an early-stage transcription factor, TFAP2C, and multiple late-stage transcription factors. Loss-of-function studies confirm the requirement of all transcription factors identified for adequate differentiation, and we reveal that the dynamic changes in the levels of TFAP2C are essential. Notably, TFAP2C pre-occupies the regulatory elements of the inactive extravillous trophoblast-active genes during the early stage of differentiation, and the late-stage transcription factors directly activate extravillous trophoblast-active genes, including themselves as differentiation further progresses, suggesting sequential actions of transcription factors assuring differentiation. Our results reveal stage-specific transcription factors and their inter-connected regulatory mechanisms modulating extravillous trophoblast differentiation, providing a framework for understanding early human placentation and placenta-related complications.

Age-dependent genes in adipose stem and precursor cells affect regulation of fat cell differentiation and link aging to obesity via cellular and genetic interactions

BACKGROUND

Age and obesity are dominant risk factors for several common cardiometabolic disorders, and both are known to impair adipose tissue function. However, the underlying cellular and genetic factors linking aging and obesity on adipose tissue function have remained elusive. Adipose stem and precursor cells (ASPCs) are an understudied, yet crucial adipose cell type due to their deterministic adipocyte differentiation potential, which impacts the capacity to store fat in a metabolically healthy manner.

METHODS

We integrated subcutaneous adipose tissue (SAT) bulk (n=435) and large single-nucleus RNA sequencing (n=105) data with the UK Biobank (UKB) (n=391,701) data to study age-obesity interactions originating from ASPCs by performing cell-type decomposition, differential expression testing, cell-cell communication analyses, and construction of polygenic risk scores for body mass index (BMI).

RESULTS

We found that the SAT ASPC proportions significantly decrease with age in an obesity-dependent way consistently in two independent cohorts, both showing that the age dependency of ASPC proportions is abolished by obesity. We further identified 76 genes (72 SAT ASPC marker genes and 4 transcription factors regulating ASPC marker genes) that are differentially expressed by age in SAT and functionally enriched for developmental processes and adipocyte differentiation (i.e., adipogenesis). The 76 age-perturbed ASPC genes include multiple negative regulators of adipogenesis, such as RORA, SMAD3, TWIST2, and ZNF521, form tight clusters of longitudinally co-expressed genes during human adipogenesis, and show age-based differences in cellular interactions between ASPCs and adipose cell types. Finally, our genetic data demonstrate that cis-regional variants of these genes interact with age as predictors of BMI in an obesity-dependent way in the large UKB, while no such gene-age interaction on BMI is observed with non-age-dependent ASPC marker genes, thus independently confirming our cellular ASPC results at the biobank level.

CONCLUSIONS

Overall, we discover that obesity prematurely induces a decrease in ASPC proportions and identify 76 developmentally important ASPC genes that implicate altered negative regulation of fat cell differentiation as a mechanism for aging and directly link aging to obesity via significant cellular and genetic interactions.

JMnorm: a novel joint multi-feature normalization method for integrative and comparative epigenomics

Combinatorial patterns of epigenetic features reflect transcriptional states and functions of genomic regions. While many epigenetic features have correlated relationships, most existing data normalization approaches analyze each feature independently. Such strategies may distort relationships between functionally correlated epigenetic features and hinder biological interpretation. We present a novel approach named JMnorm that simultaneously normalizes multiple epigenetic features across cell types, species, and experimental conditions by leveraging information from partially correlated epigenetic features. We demonstrate that JMnorm-normalized data can better preserve cross-epigenetic-feature correlations across different cell types and enhance consistency between biological replicates than data normalized by other methods. Additionally, we show that JMnorm-normalized data can consistently improve the performance of various downstream analyses, which include candidate cis-regulatory element clustering, cross-cell-type gene expression prediction, detection of transcription factor binding and changes upon perturbations. These findings suggest that JMnorm effectively minimizes technical noise while preserving true biologically significant relationships between epigenetic datasets. We anticipate that JMnorm will enhance integrative and comparative epigenomics.

Global analysis of the abundance of AU-rich mRNAs in response to glucocorticoid treatment

Glucocorticoids (GC) like dexamethasone (Dex) are potent anti-inflammatory agents with diverse cellular functions including the potentiation of the activity of AU-rich elements (AREs). AREs are cis-acting instability sequence elements located in the 3'UTRs of many inflammatory mediator mRNAs. Here, available RNA-seq data were used to investigate the effect of GCs on the ARE-mRNA-transcriptome. At a global scale, ARE-mRNAs had a tendency to be downregulated after GC-treatment of the A549 lung cancer cell-line, but with notable cases of upregulation. mRNA stability experiments indicated that not only the downregulated, but also the upregulated ARE-mRNAs are destabilized by Dex-treatment. Several of the most upregulated ARE-mRNAs code for anti-inflammatory mediators including the established GC targets DUSP1 and ZFP36; both code for proteins that target ARE-containing mRNAs for destruction. GCs are widely used in the treatment of COVID-19 patients; we show that ARE-mRNAs are more likely to regulate in opposite directions between Dex-treatment and SARS-CoV-2 infections compared to non-ARE mRNAs. The effect of GC treatment on ARE-mRNA abundance was also investigated in blood monocytes of COVID-19 patients. The results were heterogeneous; however, in agreement with in vitro observations, ZFP36 and DUSP1 were often amongst the most differentially expressed mRNAs. The results of this study propose a universal destabilization of ARE-mRNAs by GCs, but a diverse overall outcome in vitro likely due to induced transcription or due to the heterogeneity of COVID-19 patient's responses in vivo.

An internally controlled system to study microtubule network diversification links tubulin evolution to the use of distinct microtubule regulators

Diverse eukaryotic cells assemble microtubule networks that vary in structure and composition. While we understand how cells build microtubule networks with specialized functions, we do not know how microtubule networks diversify across deep evolutionary timescales. This problem has remained unresolved because most organisms use shared pools of tubulins for multiple networks, making it impossible to trace the evolution of any single network. In contrast, the amoeboflagellate Naegleria uses distinct tubulin genes to build distinct microtubule networks: while Naegleria builds flagella from conserved tubulins during differentiation, it uses divergent tubulins to build its mitotic spindle. This genetic separation makes for an internally controlled system to study independent microtubule networks in a single organismal and genomic context. To explore the evolution of these microtubule networks, we identified conserved microtubule binding proteins and used transcriptional profiling of mitosis and differentiation to determine which are upregulated during the assembly of each network. Surprisingly, most microtubule binding proteins are upregulated during only one process, suggesting that Naegleria uses distinct component pools to specialize its microtubule networks. Furthermore, the divergent residues of mitotic tubulins tend to fall within the binding sites of differentiation-specific microtubule regulators, suggesting that interactions between microtubules and their binding proteins constrain tubulin sequence diversification. We therefore propose a model for cytoskeletal evolution in which pools of microtubule network components constrain and guide the diversification of the entire network, so that the evolution of tubulin is inextricably linked to that of its binding partners.

Maresin 1 repletion improves muscle regeneration after volumetric muscle loss

The acute traumatic or surgical loss of skeletal muscle, known as volumetric muscle loss (VML), is a devastating type of injury that results in exacerbated and persistent inflammation followed by fibrosis. The mechanisms that mediate the magnitude and duration of the inflammatory response and ensuing fibrosis after VML remain understudied, and as such, the development of regenerative therapies has been limited. To address this need, we profiled how lipid mediators, which are potent regulators of the immune response after injury, varied with VML injuries that heal or result in fibrosis. We observed that non-healing VML injuries displayed increased pro-inflammatory eicosanoids and a lack of pro-resolving lipid mediators. Treatment of VML with a pro-resolving lipid mediator synthesized from docosahexaenoic acid, called Maresin 1, ameliorated fibrosis through reduction of neutrophils and macrophages and enhanced recovery of muscle strength. These results expand our knowledge of the dysregulated immune response that develops after VML and identify a novel immuno-regenerative therapeutic modality in Maresin 1.

Direct androgen receptor control of sexually dimorphic gene expression in the mammalian kidney

Mammalian organs exhibit distinct physiology, disease susceptibility, and injury responses between the sexes. In the mouse kidney, sexually dimorphic gene activity maps predominantly to proximal tubule (PT) segments. Bulk RNA sequencing (RNA-seq) data demonstrated that sex differences were established from 4 and 8 weeks after birth under gonadal control. Hormone injection studies and genetic removal of androgen and estrogen receptors demonstrated androgen receptor (AR)-mediated regulation of gene activity in PT cells as the regulatory mechanism. Interestingly, caloric restriction feminizes the male kidney. Single-nuclear multiomic analysis identified putative cis-regulatory regions and cooperating factors mediating PT responses to AR activity in the mouse kidney. In the human kidney, a limited set of genes showed conserved sex-linked regulation, whereas analysis of the mouse liver underscored organ-specific differences in the regulation of sexually dimorphic gene expression. These findings raise interesting questions on the evolution, physiological significance, disease, and metabolic linkage of sexually dimorphic gene activity.

Mediator Kinase Inhibition Impedes Transcriptional Plasticity and Prevents Resistance to ERK/MAPK-Targeted Therapy in KRAS-Mutant Cancers

Acquired resistance remains a major challenge for therapies targeting oncogene activated pathways. KRAS is the most frequently mutated oncogene in human cancers, yet strategies targeting its downstream signaling kinases have failed to produce durable treatment responses. Here, we developed multiple models of acquired resistance to dual-mechanism ERK/MAPK inhibitors across KRAS-mutant pancreatic, colorectal, and lung cancers, and then probed the long-term events enabling survival against this class of drugs. These studies revealed that resistance emerges secondary to large-scale transcriptional adaptations that are diverse and cell line-specific. Transcriptional reprogramming extends beyond the well-established early response, and instead represents a dynamic, evolved process that is refined to attain a stably resistant phenotype. Mechanistic and translational studies reveal that resistance to dual-mechanism ERK/MAPK inhibition is broadly susceptible to manipulation of the epigenetic machinery, and that Mediator kinase, in particular, can be co-targeted at a bottleneck point to prevent diverse, cell line-specific resistance programs.

Cellular clarity: a logistic regression approach to identify root epidermal regulators of iron deficiency response

BACKGROUND

Plants respond to stress through highly tuned regulatory networks. While prior works identified master regulators of iron deficiency responses in A. thaliana from whole-root data, identifying regulators that act at the cellular level is critical to a more comprehensive understanding of iron homeostasis. Within the root epidermis complex molecular mechanisms that facilitate iron reduction and uptake from the rhizosphere are known to be regulated by bHLH transcriptional regulators. However, many questions remain about the regulatory mechanisms that control these responses, and how they may integrate with developmental processes within the epidermis. Here, we use transcriptional profiling to gain insight into root epidermis-specific regulatory processes.

RESULTS

Set comparisons of differentially expressed genes (DEGs) between whole root and epidermis transcript measurements identified differences in magnitude and timing of organ-level vs. epidermis-specific responses. Utilizing a unique sampling method combined with a mutual information metric across time-lagged and non-time-lagged windows, we identified relationships between clusters of functionally relevant differentially expressed genes suggesting that developmental regulatory processes may act upstream of well-known Fe-specific responses. By integrating static data (DNA motif information) with time-series transcriptomic data and employing machine learning approaches, specifically logistic regression models with LASSO, we also identified putative motifs that served as crucial features for predicting differentially expressed genes. Twenty-eight transcription factors (TFs) known to bind to these motifs were not differentially expressed, indicating that these TFs may be regulated post-transcriptionally or post-translationally. Notably, many of these TFs also play a role in root development and general stress response.

CONCLUSIONS

This work uncovered key differences in -Fe response identified using whole root data vs. cell-specific root epidermal data. Machine learning approaches combined with additional static data identified putative regulators of -Fe response that would not have been identified solely through transcriptomic profiles and reveal how developmental and general stress responses within the epidermis may act upstream of more specialized -Fe responses for Fe uptake.

The dynamics of touch-responsive gene expression in cereals

Wind, rain, herbivores, obstacles, neighbouring plants, etc. provide important mechanical cues to steer plant growth and survival. Mechanostimulation to stimulate yield and stress resistance of crops is of significant research interest, yet a molecular understanding of transcriptional responses to touch is largely absent in cereals. To address this, we performed whole-genome transcriptomics following mechanostimulation of wheat, barley, and the recent genome-sequenced oat. The largest transcriptome changes occurred ±25 min after touching, with most of the genes being upregulated. While most genes returned to basal expression level by 1-2 h in oat, many genes retained high expression even 4 h post-treatment in barley and wheat. Functional categories such as transcription factors, kinases, phytohormones, and Ca2+ regulation were affected. In addition, cell wall-related genes involved in (hemi)cellulose, lignin, suberin, and callose biosynthesis were touch-responsive, providing molecular insight into mechanically induced changes in cell wall composition. Furthermore, several cereal-specific transcriptomic footprints were identified that were not observed in Arabidopsis. In oat and barley, we found evidence for systemic spreading of touch-induced signalling. Finally, we provide evidence that both the jasmonic acid-dependent and the jasmonic acid-independent pathways underlie touch-signalling in cereals, providing a detailed framework and marker genes for further study of (a)biotic stress responses in cereals.

Regulatory response to a hybrid ancestral nitrogenase in Azotobacter vinelandii

Biological nitrogen fixation, the microbial reduction of atmospheric nitrogen to bioavailable ammonia, represents both a major limitation on biological productivity and a highly desirable engineering target for synthetic biology. However, the engineering of nitrogen fixation requires an integrated understanding of how the gene regulatory dynamics of host diazotrophs respond across sequence-function space of its central catalytic metalloenzyme, nitrogenase. Here, we interrogate this relationship by analyzing the transcriptome of Azotobacter vinelandii engineered with a phylogenetically inferred ancestral nitrogenase protein variant. The engineered strain exhibits reduced cellular nitrogenase activity but recovers wild-type growth rates following an extended lag period. We find that expression of genes within the immediate nitrogen fixation network is resilient to the introduced nitrogenase sequence-level perturbations. Rather the sustained physiological compatibility with the ancestral nitrogenase variant is accompanied by reduced expression of genes that support trace metal and electron resource allocation to nitrogenase. Our results spotlight gene expression changes in cellular processes adjacent to nitrogen fixation as productive engineering considerations to improve compatibility between remodeled nitrogenase proteins and engineered host diazotrophs. IMPORTANCE Azotobacter vinelandii is a key model bacterium for the study of biological nitrogen fixation, an important metabolic process catalyzed by nitrogenase enzymes. Here, we demonstrate that compatibilities between engineered A. vinelandii strains and nitrogenase variants can be modulated at the regulatory level. The engineered strain studied here responds by adjusting the expression of proteins involved in cellular processes adjacent to nitrogen fixation, rather than that of nitrogenase proteins themselves. These insights can inform future strategies to transfer nitrogenase variants to non-native hosts.

Nonlinear expression patterns and multiple shifts in gene network interactions underlie robust phenotypic change in Drosophila melanogaster selected for night sleep duration

All but the simplest phenotypes are believed to result from interactions between two or more genes forming complex networks of gene regulation. Sleep is a complex trait known to depend on the system of feedback loops of the circadian clock, and on many other genes; however, the main components regulating the phenotype and how they interact remain an unsolved puzzle. Genomic and transcriptomic data may well provide part of the answer, but a full account requires a suitable quantitative framework. Here we conducted an artificial selection experiment for sleep duration with RNA-seq data acquired each generation. The phenotypic results are robust across replicates and previous experiments, and the transcription data provides a high-resolution, time-course data set for the evolution of sleep-related gene expression. In addition to a Hierarchical Generalized Linear Model analysis of differential expression that accounts for experimental replicates we develop a flexible Gaussian Process model that estimates interactions between genes. 145 gene pairs are found to have interactions that are different from controls. Our method appears to be not only more specific than standard correlation metrics but also more sensitive, finding correlations not significant by other methods. Statistical predictions were compared to experimental data from public databases on gene interactions. Mutations of candidate genes implicated by our results affected night sleep, and gene expression profiles largely met predicted gene-gene interactions.

Single-cell transcriptomic profiling of human pancreatic islets reveals genes responsive to glucose exposure over 24 hours

Disruption of pancreatic islet function and glucose homeostasis can lead to the development of sustained hyperglycemia, beta cell glucotoxicity, and ultimately type 2 diabetes (T2D). In this study, we sought to explore the effects of hyperglycemia on human pancreatic islet (HPI) gene expression by exposing HPIs from two donors to low (2.8mM) and high (15.0mM) glucose concentrations over 24 hours, assaying the transcriptome at seven time points using single-cell RNA sequencing (scRNA-seq). We modeled time as both a discrete and continuous variable to determine momentary and longitudinal changes in transcription associated with islet time in culture or glucose exposure. Across all cell types, we identified 1,528 genes associated with time, 1,185 genes associated with glucose exposure, and 845 genes associated with interaction effects between time and glucose. We clustered differentially expressed genes across cell types and found 347 modules of genes with similar expression patterns across time and glucose conditions, including two beta cell modules enriched in genes associated with T2D. Finally, by integrating genomic features from this study and genetic summary statistics for T2D and related traits, we nominate 363 candidate effector genes that may underlie genetic associations for T2D and related traits.

Learning perturbation-inducible cell states from observability analysis of transcriptome dynamics

A major challenge in biotechnology and biomanufacturing is the identification of a set of biomarkers for perturbations and metabolites of interest. Here, we develop a data-driven, transcriptome-wide approach to rank perturbation-inducible genes from time-series RNA sequencing data for the discovery of analyte-responsive promoters. This provides a set of biomarkers that act as a proxy for the transcriptional state referred to as cell state. We construct low-dimensional models of gene expression dynamics and rank genes by their ability to capture the perturbation-specific cell state using a novel observability analysis. Using this ranking, we extract 15 analyte-responsive promoters for the organophosphate malathion in the underutilized host organism Pseudomonas fluorescens SBW25. We develop synthetic genetic reporters from each analyte-responsive promoter and characterize their response to malathion. Furthermore, we enhance malathion reporting through the aggregation of the response of individual reporters with a synthetic consortium approach, and we exemplify the library's ability to be useful outside the lab by detecting malathion in the environment. The engineered host cell, a living malathion sensor, can be optimized for use in environmental diagnostics while the developed machine learning tool can be applied to discover perturbation-inducible gene expression systems in the compendium of host organisms.

Direct androgen receptor regulation of sexually dimorphic gene expression in the mammalian kidney

Mammalian organs exhibit distinct physiology, disease susceptibility and injury responses between the sexes. In the mouse kidney, sexually dimorphic gene activity maps predominantly to proximal tubule (PT) segments. Bulk RNA-seq data demonstrated sex differences were established from 4 and 8 weeks after birth under gonadal control. Hormone injection studies and genetic removal of androgen and estrogen receptors demonstrated androgen receptor (AR) mediated regulation of gene activity in PT cells as the regulatory mechanism. Interestingly, caloric restriction feminizes the male kidney. Single-nuclear multiomic analysis identified putative cis-regulatory regions and cooperating factors mediating PT responses to AR activity in the mouse kidney. In the human kidney, a limited set of genes showed conserved sex-linked regulation while analysis of the mouse liver underscored organ-specific differences in the regulation of sexually dimorphic gene expression. These findings raise interesting questions on the evolution, physiological significance, and disease and metabolic linkage, of sexually dimorphic gene activity.

Snapshot: a package for clustering and visualizing epigenetic history during cell differentiation

BACKGROUND

Epigenetic modification of chromatin plays a pivotal role in regulating gene expression during cell differentiation. The scale and complexity of epigenetic data pose significant challenges for biologists to identify the regulatory events controlling cell differentiation.

RESULTS

To reduce the complexity, we developed a package, called Snapshot, for clustering and visualizing candidate cis-regulatory elements (cCREs) based on their epigenetic signals during cell differentiation. This package first introduces a binarized indexing strategy for clustering the cCREs. It then provides a series of easily interpretable figures for visualizing the signal and epigenetic state patterns of the cCREs clusters during the cell differentiation. It can also use different hierarchies of cell types to highlight the epigenetic history specific to any particular cell lineage. We demonstrate the utility of Snapshot using data from a consortium project for ValIdated Systematic IntegratiON (VISION) of epigenomic data in hematopoiesis.

CONCLUSION

The package Snapshot can identify all distinct clusters of genomic locations with unique epigenetic signal patterns during cell differentiation. It outperforms other methods in terms of interpreting and reproducing the identified cCREs clusters. The package of Snapshot is available at GitHub: https://github.com/guanjue/Snapshot .

AngClust: Angle Feature-Based Clustering for Short Time Series Gene Expression Profiles

When clustering gene expression, it is expected that correlation coefficients of genes in the same clusters are high, and that gene ontology (GO) enrichment analysis of most clusters will be significant. However, existing short-term gene expression clustering algorithms have limitations. To address this problem, we proposed a novel clustering process based on angular features for short-term gene expression. Our method (named AngClust) uses angular features to indicate the change of trend in gene expression levels at two neighboring time points. The changes of angles at multiple time points reflects the change of trend of the overall expression levels. Such changes are used to measure whether the expression trends of different genes are similar. To obtain functionally significant clusters from the clustering results, we evaluated numbers of genes in clusters, average correlation coefficient, fluctuation, and their correlation with GO term enrichment. The efficacy of AngClust outperform two other measures, Euclidean distance (ED) and dynamic time warping of correlation (DTW), on a dataset of yeast gene expression. The ratios of GO and pathway term-enriched of clusters of AngClust is higher than or equal to that of STEM and TMixClust on human, mouse, and yeast time series of gene expression.

A Parallel Multiobjective PSO Weighted Average Clustering Algorithm Based on Apache Spark

Multiobjective clustering algorithm using particle swarm optimization has been applied successfully in some applications. However, existing algorithms are implemented on a single machine and cannot be directly parallelized on a cluster, which makes it difficult for existing algorithms to handle large-scale data. With the development of distributed parallel computing framework, data parallelism was proposed. However, the increase in parallelism will lead to the problem of unbalanced data distribution affecting the clustering effect. In this paper, we propose a parallel multiobjective PSO weighted average clustering algorithm based on apache Spark (Spark-MOPSO-Avg). First, the entire data set is divided into multiple partitions and cached in memory using the distributed parallel and memory-based computing of Apache Spark. The local fitness value of the particle is calculated in parallel according to the data in the partition. After the calculation is completed, only particle information is transmitted, and there is no need to transmit a large number of data objects between each node, reducing the communication of data in the network and thus effectively reducing the algorithm's running time. Second, a weighted average calculation of the local fitness values is performed to improve the problem of unbalanced data distribution affecting the results. Experimental results show that the Spark-MOPSO-Avg algorithm achieves lower information loss under data parallelism, losing about 1% to 9% accuracy, but can effectively reduce the algorithm time overhead. It shows good execution efficiency and parallel computing capability under the Spark distributed cluster.

Production of rAAV by plasmid transfection induces antiviral and inflammatory responses in suspension HEK293 cells

Recombinant adeno-associated virus (rAAV) is a clinically proven viral vector for delivery of therapeutic genes to treat rare diseases. Improving rAAV manufacturing productivity and vector quality is necessary to meet clinical and commercial demand. These goals will require an improved understanding of the cellular response to rAAV production, which is poorly defined. We interrogated the kinetic transcriptional response of HEK293 cells to rAAV production following transient plasmid transfection, under manufacturing-relevant conditions, using RNA-seq. Time-series analyses identified a robust cellular response to transfection and rAAV production, with 1,850 transcripts differentially expressed. Gene Ontology analysis determined upregulated pathways, including inflammatory and antiviral responses, with several interferon-stimulated cytokines and chemokines being upregulated at the protein level. Literature-based pathway prediction implicated multiple pathogen pattern sensors and signal transducers in up-regulation of inflammatory and antiviral responses in response to transfection and rAAV replication. Systematic analysis of the cellular transcriptional response to rAAV production indicates that host cells actively sense vector manufacture as an infectious insult. This dataset may therefore illuminate genes and pathways that influence rAAV production, thereby enabling the rational design of next-generation manufacturing platforms to support safe, effective, and affordable AAV-based gene therapies.

Deep longitudinal multi-omics analysis of Bordetella pertussis cultivated in bioreactors highlights medium starvations and transitory metabolisms, associated to vaccine antigen biosynthesis variations and global virulence regulation

Bordetella pertussis is the bacterial causative agent of whooping cough, a serious respiratory illness. An extensive knowledge on its virulence regulation and metabolism is a key factor to ensure pertussis vaccine manufacturing process robustness. The aim of this study was to refine our comprehension of B. pertussis physiology during in vitro cultures in bioreactors. A longitudinal multi-omics analysis was carried out over 26 h small-scale cultures of B. pertussis. Cultures were performed in batch mode and under culture conditions intending to mimic industrial processes. Putative cysteine and proline starvations were, respectively, observed at the beginning of the exponential phase (from 4 to 8 h) and during the exponential phase (18 h 45 min). As revealed by multi-omics analyses, the proline starvation induced major molecular changes, including a transient metabolism with internal stock consumption. In the meantime, growth and specific total PT, PRN, and Fim2 antigen productions were negatively affected. Interestingly, the master virulence-regulating two-component system of B. pertussis (BvgASR) was not evidenced as the sole virulence regulator in this in vitro growth condition. Indeed, novel intermediate regulators were identified as putatively involved in the expression of some virulence-activated genes (vags). Such longitudinal multi-omics analysis applied to B. pertussis culture process emerges as a powerful tool for characterization and incremental optimization of vaccine antigen production.

Transcriptome Profiling of Stem-Differentiating Xylem in Response to Abiotic Stresses Based on Hybrid Sequencing in Cunninghamia lanceolata

Cunninghamia lanceolata (C. lanceolata) belongs to Gymnospermae, which are fast-growing and have desirable wood properties. However, C. lanceolata's stress resistance is little understood. To unravel the physiological and molecular regulation mechanisms under environmental stresses in the typical gymnosperm species of C. lanceolata, three-year-old plants were exposed to simulated drought stress (polyethylene glycol 8000), salicylic acid, and cold treatment at 4 °C for 8 h, 32 h, and 56 h, respectively. Regarding the physiological traits, we observed a decreased protein content and increased peroxidase upon salicylic acid and polyethylene glycol treatment. Superoxide dismutase activity either decreased or increased at first and then returned to normal under the stresses. Regarding the molecular regulation, we used both nanopore direct RNA sequencing and short-read sequencing to reveal a total of 5646 differentially expressed genes in response to different stresses, of which most had functions in lignin catabolism, pectin catabolism, and xylan metabolism, indicating that the development of stem-differentiating xylem was affected upon stress treatment. Finally, we identified a total of 51 AP2/ERF, 29 NAC, and 37 WRKY transcript factors in C. lanceolata. The expression of most of the NAC TFs increased under cold stress, and the expression of most of the WRKY TFs increased under cold and SA stress. These results revealed the transcriptomics responses in C. lanceolata to short-term stresses under this study's experimental conditions and provide preliminary clues about stem-differentiating xylem changes associated with different stresses.

Unsupervised probabilistic models for sequential Electronic Health Records

We develop an unsupervised probabilistic model for heterogeneous Electronic Health Record (EHR) data. Utilizing a mixture model formulation, our approach directly models sequences of arbitrary length, such as medications and laboratory results. This allows for subgrouping and incorporation of the dynamics underlying heterogeneous data types. The model consists of a layered set of latent variables that encode underlying structure in the data. These variables represent subject subgroups at the top layer, and unobserved states for sequences in the second layer. We train this model on episodic data from subjects receiving medical care in the Kaiser Permanente Northern California integrated healthcare delivery system. The resulting properties of the trained model generate novel insight from these complex and multifaceted data. In addition, we show how the model can be used to analyze sequences that contribute to assessment of mortality likelihood.

Time-Series-Based Personalized Lane-Changing Decision-Making Model

In recent years, autonomous driving technology has been changing from "human adapting to vehicle" to "vehicle adapting to human". To improve the adaptability of autonomous driving systems to human drivers, a time-series-based personalized lane change decision (LCD) model is proposed. Firstly, according to the characteristics of the subject vehicle (SV) with respect to speed, acceleration and headway, an unsupervised clustering algorithm, namely, a Gaussian mixture model (GMM), is used to identify its three different driving styles. Secondly, considering the interaction between the SV and the surrounding vehicles, the lane change (LC) gain value is produced by developing a gain function to characterize their interaction. On the basis of the recognition of the driving style, this gain value and LC feature parameters are employed as model inputs to develop a personalized LCD model on the basis of a long short-term memory (LSTM) recurrent neural network model (RNN). The proposed method is tested using the US Open Driving Dataset NGSIM. The results show that the accuracy, F1 score, and macro-average area under the curve (macro-AUC) value of the proposed method for LC behavior prediction are 0.965, 0.951 and 0.983, respectively, and the performance is significantly better than that of other mainstream models. At the same time, the method is able to capture the LCD behavior of different human drivers, enabling personalized driving.

A species-specific lncRNA modulates the reproductive ability of the asian tiger mosquito

Long non-coding RNA (lncRNA) research has emerged as an independent scientific field in recent years. Despite their association with critical cellular and metabolic processes in plenty of organisms, lncRNAs are still a largely unexplored area in mosquito research. We propose that they could serve as exceptional tools for pest management due to unique features they possess. These include low inter-species sequence conservation and high tissue specificity. In the present study, we investigated the role of ovary-specific lncRNAs in the reproductive ability of the Asian tiger mosquito, Aedes albopictus. Through the analysis of transcriptomic data, we identified several lncRNAs that were differentially expressed upon blood feeding; we called these genes Norma (NOn-coding RNA in Mosquito ovAries). We observed that silencing some of these Normas resulted in significant impact on mosquito fecundity and fertility. We further focused on Norma3 whose silencing resulted in 43% oviposition reduction, in smaller ovaries and 53% hatching reduction of the laid eggs, compared to anti-GFP controls. Moreover, a significant downregulation of 2 mucins withing a neighboring (∼100 Kb) mucin cluster was observed in smaller anti-Norma3 ovaries, indicating a potential mechanism of in-cis regulation between Norma3 and the mucins. Our work constitutes the first experimental proof-of-evidence connecting lncRNAs with mosquito reproduction and opens a novel path for pest management.

Constrained Fourier estimation of short-term time-series gene expression data reduces noise and improves clustering and gene regulatory network predictions

BACKGROUND

Biological data suffers from noise that is inherent in the measurements. This is particularly true for time-series gene expression measurements. Nevertheless, in order to to explore cellular dynamics, scientists employ such noisy measurements in predictive and clustering tools. However, noisy data can not only obscure the genes temporal patterns, but applying predictive and clustering tools on noisy data may yield inconsistent, and potentially incorrect, results.

RESULTS

To reduce the noise of short-term (< 48 h) time-series expression data, we relied on the three basic temporal patterns of gene expression: waves, impulses and sustained responses. We constrained the estimation of the true signals to these patterns by estimating the parameters of first and second-order Fourier functions and using the nonlinear least-squares trust-region optimization technique. Our approach lowered the noise in at least 85% of synthetic time-series expression data, significantly more than the spline method ([Formula: see text]). When the data contained a higher signal-to-noise ratio, our method allowed downstream network component analyses to calculate consistent and accurate predictions, particularly when the noise variance was high. Conversely, these tools led to erroneous results from untreated noisy data. Our results suggest that at least 5-7 time points are required to efficiently de-noise logarithmic scaled time-series expression data. Investing in sampling additional time points provides little benefit to clustering and prediction accuracy.

CONCLUSIONS

Our constrained Fourier de-noising method helps to cluster noisy gene expression and interpret dynamic gene networks more accurately. The benefit of noise reduction is large and can constitute the difference between a successful application and a failing one.

The Time Sequence of Gene Expression Changes after Spinal Cord Injury

Gene expression changes following spinal cord injury (SCI) are time-dependent, and an accurate understanding of these changes can be crucial in determining time-based treatment options in a clinical setting. We performed RNA sequencing of the contused spinal cord of rats at five different time points from the very acute to chronic stages (1 hour, 1 day, 1 week, 1 month, and 3 months) following SCI. We identified differentially expressed genes (DEGs) and Gene Ontology (GO) terms at each time point, and 14,257 genes were commonly expressed at all time points. The biological process of the inflammatory response was increased at 1 hour and 1 day, and the cellular component of the integral component of the synaptic membrane was increased at 1 day. DEGs associated with cell activation and the innate immune response were highly enriched at 1 week and 1 month, respectively. A total of 2841 DEGs were differentially expressed at any of the five time points, and 18 genes (17 upregulated and 1 downregulated) showed common expression differences at all time points. We found that interleukin signaling, neutrophil degranulation, eukaryotic translation, collagen degradation, LGI–ADAM interactions, GABA receptor, and L1CAM-ankyrin interactions were prominent after SCI depending on the time post injury. We also performed gene–drug network analysis and found several potential antagonists and agonists which can be used to treat SCI. We expect to discover effective treatments in the clinical field through further studies revealing the efficacy and safety of potential drugs.

Long non-coding RNA HOTAIR regulates cytoskeleton remodeling and lipid storage capacity during adipogenesis

The long non-coding RNA HOTAIR is the most differentially expressed gene between upper- and lower-body adipose tissue, yet its functional significance in adipogenesis is unclear. We report that HOTAIR expression is transiently induced during early adipogenic differentiation of gluteofemoral adipose progenitors and repressed in mature adipocytes. Upon adipogenic commitment, HOTAIR regulates protein synthesis pathways and cytoskeleton remodeling with a later impact on mature adipocyte lipid storage capacity. Our results support novel and important functions of HOTAIR in the physiological context of adipogenesis.

Local euchromatin enrichment in lamina-associated domains anticipates their repositioning in the adipogenic lineage

Background

Interactions of chromatin with the nuclear lamina via lamina-associated domains (LADs) confer structural stability to the genome. The dynamics of positioning of LADs during differentiation, and how LADs impinge on developmental gene expression, remains, however, elusive.

Results

We examined changes in the association of lamin B1 with the genome in the first 72 h of differentiation of adipose stem cells into adipocytes. We demonstrate a repositioning of entire stand-alone LADs and of LAD edges as a prominent nuclear structural feature of early adipogenesis. Whereas adipogenic genes are released from LADs, LADs sequester downregulated or repressed genes irrelevant for the adipose lineage. However, LAD repositioning only partly concurs with gene expression changes. Differentially expressed genes in LADs, including LADs conserved throughout differentiation, reside in local euchromatic and lamin-depleted sub-domains. In these sub-domains, pre-differentiation histone modification profiles correlate with the LAD versus inter-LAD outcome of these genes during adipogenic commitment. Lastly, we link differentially expressed genes in LADs to short-range enhancers which overall co-partition with these genes in LADs versus inter-LADs during differentiation.

Conclusions

We conclude that LADs are predictable structural features of adipose nuclear architecture that restrain non-adipogenic genes in a repressive environment.

Supplementary Information

The online version contains supplementary material available at 10.1186/s13059-022-02662-6.

Time Series Transcriptome Analysis in Medicago truncatula Shoot and Root Tissue During Early Nodulation

In response to colonization by rhizobia bacteria, legumes are able to form nitrogen-fixing nodules in their roots, allowing the plants to grow efficiently in nitrogen-depleted environments. Legumes utilize a complex, long-distance signaling pathway to regulate nodulation that involves signals in both roots and shoots. We measured the transcriptional response to treatment with rhizobia in both the shoots and roots of Medicago truncatula over a 72-h time course. To detect temporal shifts in gene expression, we developed GeneShift, a novel computational statistics and machine learning workflow that addresses the time series replicate the averaging issue for detecting gene expression pattern shifts under different conditions. We identified both known and novel genes that are regulated dynamically in both tissues during early nodulation including leginsulin, defensins, root transporters, nodulin-related, and circadian clock genes. We validated over 70% of the expression patterns that GeneShift discovered using an independent M. truncatula RNA-Seq study. GeneShift facilitated the discovery of condition-specific temporally differentially expressed genes in the symbiotic nodulation biological system. In principle, GeneShift should work for time-series gene expression profiling studies from other systems.

ALOHA: Aggregated local extrema splines for high-throughput dose-response analysis

Computational methods for genomic dose-response integrate dose-response modeling with bioinformatics tools to evaluate changes in molecular and cellular functions related to pathogenic processes. These methods use parametric models to describe each gene's dose-response, but such models may not adequately capture expression changes. Additionally, current approaches do not consider gene co-expression networks. When assessing co-expression networks, one typically does not consider the dose-response relationship, resulting in 'co-regulated' gene sets containing genes having different dose-response patterns. To avoid these limitations, we develop an analysis pipeline called Aggregated Local Extrema Splines for High-throughput Analysis (ALOHA), which computes individual genomic dose-response functions using a flexible class Bayesian shape constrained splines and clusters gene co-regulation based upon these fits. Using splines, we reduce information loss due to parametric lack-of-fit issues, and because we cluster on dose-response relationships, we better identify co-regulation clusters for genes that have co-expressed dose-response patterns from chemical exposure. The clustered pathways can then be used to estimate a dose associated with a pre-specified biological response, i.e., the benchmark dose (BMD), and approximate a point of departure dose corresponding to minimal adverse response in the whole tissue/organism. We compare our approach to current parametric methods and our biologically enriched gene sets to cluster on normalized expression data. Using this methodology, we can more effectively extract the underlying structure leading to more cohesive estimates of gene set potency.

Identifying temporal and spatial patterns of variation from multimodal data using MEFISTO

Factor analysis is a widely used method for dimensionality reduction in genome biology, with applications from personalized health to single-cell biology. Existing factor analysis models assume independence of the observed samples, an assumption that fails in spatio-temporal profiling studies. Here we present MEFISTO, a flexible and versatile toolbox for modeling high-dimensional data when spatial or temporal dependencies between the samples are known. MEFISTO maintains the established benefits of factor analysis for multimodal data, but enables the performance of spatio-temporally informed dimensionality reduction, interpolation, and separation of smooth from non-smooth patterns of variation. Moreover, MEFISTO can integrate multiple related datasets by simultaneously identifying and aligning the underlying patterns of variation in a data-driven manner. To illustrate MEFISTO, we apply the model to different datasets with spatial or temporal resolution, including an evolutionary atlas of organ development, a longitudinal microbiome study, a single-cell multi-omics atlas of mouse gastrulation and spatially resolved transcriptomics.

MEFISTO models bulk and single-cell multi-omics data with temporal or spatial dependencies for interpretable pattern discovery and integration.

New insights into the evolution of SPX gene family from algae to legumes; a focus on soybean

BACKGROUND

SPX-containing proteins have been known as key players in phosphate signaling and homeostasis. In Arabidopsis and rice, functions of some SPXs have been characterized, but little is known about their function in other plants, especially in the legumes.

RESULTS

We analyzed SPX gene family evolution in legumes and in a number of key species from algae to angiosperms. We found that SPX harboring proteins showed fluctuations in domain fusions from algae to the angiosperms with, finally, four classes appearing and being retained in the land plants. Despite these fluctuations, Lysine Surface Cluster (KSC), and the third residue of Phosphate Binding Sites (PBS) showed complete conservation in almost all of SPXs except few proteins in Selaginella moellendorffii and Papaver sumniferum, suggesting they might have different ligand preferences. In addition, we found that the WGD/segmentally or dispersed duplication types were the most frequent contributors to the SPX expansion, and that there is a positive correlation between the amount of WGD contribution to the SPX expansion in individual species and its number of EXS genes. We could also reveal that except SPX class genes, other classes lost the collinearity relationships among Arabidopsis and legume genomes. The sub- or neo-functionalization of the duplicated genes in the legumes makes it difficult to find the functional orthologous genes. Therefore, we used two different methods to identify functional orthologs in soybean and Medicago. High variance in the dynamic and spatial expression pattern of GmSPXs proved the new or sub-functionalization in the paralogs.

CONCLUSION

This comprehensive analysis revealed how SPX gene family evolved from algae to legumes and also discovered several new domains fused to SPX domain in algae. In addition, we hypothesized that there different phosphate sensing mechanisms might occur in S. moellendorffii and P. sumniferum. Finally, we predicted putative functional orthologs of AtSPXs in the legumes, especially, orthologs of AtPHO1, involved in long-distance Pi transportation. These findings help to understand evolution of phosphate signaling and might underpin development of new legume varieties with improved phosphate use efficiency.

Not All That Glitters Is Gold: The Paradox of CO-dependent Hydrogenogenesis in Parageobacillus thermoglucosidasius

The thermophilic bacterium Parageobacillus thermoglucosidasius has recently gained interest due to its ability to catalyze the water gas shift reaction, where the oxidation of carbon monoxide (CO) is linked to the evolution of hydrogen (H₂) gas. This phenotype is largely predictable based on the presence of a genomic region coding for a carbon monoxide dehydrogenase (CODH—Coo) and hydrogen evolving hydrogenase (Phc). In this work, seven previously uncharacterized strains were cultivated under 50% CO and 50% air atmosphere. Despite the presence of the coo—phc genes in all seven strains, only one strain, Kp1013, oxidizes CO and yields H₂. The genomes of the H₂ producing strains contain unique genomic regions that code for proteins involved in nickel transport and the detoxification of catechol, a by-product of a siderophore-mediated iron acquisition system. Combined, the presence of these genomic regions could potentially drive biological water gas shift (WGS) reaction in P. thermoglucosidasius.

Single-cell transcriptomic analysis of zebrafish cranial neural crest reveals spatiotemporal regulation of lineage decisions during development

Neural crest (NC) cells migrate throughout vertebrate embryos to give rise to a huge variety of cell types, but when and where lineages emerge and their regulation remain unclear. We have performed single-cell RNA sequencing (RNA-seq) of cranial NC cells from the first pharyngeal arch in zebrafish over several stages during migration. Computational analysis combining pseudotime and real-time data reveals that these NC cells first adopt a transitional state, becoming specified mid-migration, with the first lineage decisions being skeletal and pigment, followed by neural and glial progenitors. In addition, by computationally integrating these data with RNA-seq data from a transgenic Wnt reporter line, we identify gene cohorts with similar temporal responses to Wnts during migration and show that one, Atp6ap2, is required for melanocyte differentiation. Together, our results show that cranial NC cell lineages arise progressively and uncover a series of spatially restricted cell interactions likely to regulate such cell-fate decisions.

Tatarakis et al. provide a single-cell transcriptomic timeline of cranial neural crest (NC) development in zebrafish and address long-standing questions surrounding the integration of NC cell migration and lineage specification. They find that lineages are specified mid-migration. These fate decisions correspond to shifts in Wnt signaling, and lineages rapidly segregate.

Long-term dynamics of the human oral microbiome during clinical disease progression

BACKGROUND

Oral microbiome dysbiosis is linked to overt inflammation of tooth-supporting tissues, leading to periodontitis, an oral condition that can cause tooth and bone loss. Microbiome dysbiosis has been described as a disruption in the symbiotic microbiota composition's stability that could adversely affect the host's health status. However, the precise microbiome dynamics that lead to dysbiosis and the progression of the disease are largely unknown. The objective of our study was to investigate the long-term dynamics of periodontitis progression and its connection to dysbiosis.

RESULTS

We studied three different teeth groups: sites that showed disease progression, sites that remained stable during the study, and sites that exhibited a cyclic deepening followed by spontaneous recovery. Time-series analysis revealed that communities followed a characteristic succession of bacteria clusters. Stable and fluctuating sites showed high asynchrony in the communities (i.e., different species responding dissimilarly through time) and a reordering of the communities where directional changes dominated (i.e., sample distance increases over time) in the stable sites but not in the fluctuating sites. Progressing sites exhibited low asynchrony and convergence (i.e., samples distance decreases over time). Moreover, new species were more likely to be recruited in stable samples if a close relative was not recruited previously. In contrast, progressing and fluctuating sites followed a neutral recruitment model, indicating that competition between closely related species is a significant component of species-species interactions in stable samples. Finally, periodontal treatment did not select similar communities but stabilized α-diversity, centered the abundance of different clusters to the mean, and increased community rearrangement.

CONCLUSIONS

Here, we show that ecological principles can define dysbiosis and explain the evolution and outcomes of specific microbial communities of the oral microbiome in periodontitis progression. All sites showed an ecological succession in community composition. Stable sites were characterized by high asynchrony, a reordering of the communities where directional changes dominated, and new species were more likely to be recruited if a close relative was not recruited previously. Progressing sites were characterized by low asynchrony, community convergence, and a neutral model of recruitment. Finally, fluctuating sites were characterized by high asynchrony, community convergence, and a neutral recruitment model.

Non-parametric modelling of temporal and spatial counts data from RNA-seq experiments

MOTIVATION

The negative binomial distribution has been shown to be a good model for counts data from both bulk and single-cell RNA-sequencing (RNA-seq). Gaussian process (GP) regression provides a useful non-parametric approach for modelling temporal or spatial changes in gene expression. However, currently available GP regression methods that implement negative binomial likelihood models do not scale to the increasingly large datasets being produced by single-cell and spatial transcriptomics.

RESULTS

The GPcounts package implements GP regression methods for modelling counts data using a negative binomial likelihood function. Computational efficiency is achieved through the use of variational Bayesian inference. The GP function models changes in the mean of the negative binomial likelihood through a logarithmic link function and the dispersion parameter is fitted by maximum likelihood. We validate the method on simulated time course data, showing better performance to identify changes in over-dispersed counts data than methods based on Gaussian or Poisson likelihoods. To demonstrate temporal inference, we apply GPcounts to single-cell RNA-seq datasets after pseudotime and branching inference. To demonstrate spatial inference, we apply GPcounts to data from the mouse olfactory bulb to identify spatially variable genes and compare to two published GP methods. We also provide the option of modelling additional dropout using a zero-inflated negative binomial. Our results show that GPcounts can be used to model temporal and spatial counts data in cases where simpler Gaussian and Poisson likelihoods are unrealistic.

AVAILABILITY AND IMPLEMENTATION

GPcounts is implemented using the GPflow library in Python and is available at https://github.com/ManchesterBioinference/GPcounts along with the data, code and notebooks required to reproduce the results presented here. The version used for this paper is archived at https://doi.org/10.5281/zenodo.5027066.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

A Convolutional Deep Clustering Framework for Gene Expression Time Series

The functional or regulatory processes within the cell are explicitly governed by the expression levels of a subset of its genes. Gene expression time series captures activities of individual genes over time and aids revealing underlying cellular dynamics. An important step in high-throughput gene expression time series experiment is clustering genes based on their temporal expression patterns and is conventionally achieved by unsupervised machine learning techniques. However, most of the clustering techniques either suffer from the short length of gene expression time series or ignore temporal structure of the data. In this work, we propose DeepTrust, a novel deep learning-based framework for gene expression time series clustering which can overcome these issues. DeepTrust initially transforms time series data into images to obtain richer data representations. Afterwards, a deep convolutional clustering algorithm is applied on the constructed images. Analyses on both simulated and biological data sets exhibit the efficiency of this new framework, compared to widely used clustering techniques. We also utilize enrichment analyses to illustrate the biological plausibility of the clusters detected by DeepTrust. Our code and data are available from http://github.com/tanlab/DeepTrust.

RNA binding protein RBMS3 is a common EMT effector that modulates triple-negative breast cancer progression via stabilizing PRRX1 mRNA

The epithelial-to-mesenchymal transition (EMT) has been recognized as a driving force for tumor progression in breast cancer. Recently, our group identified the RNA Binding Motif Single Stranded Interacting Protein 3 (RBMS3) to be significantly associated with an EMT transcriptional program in breast cancer. Additional expression profiling demonstrated that RBMS3 was consistently upregulated by multiple EMT transcription factors and correlated with mesenchymal gene expression in breast cancer cell lines. Functionally, RBMS3 was sufficient to induce EMT in two immortalized mammary epithelial cell lines. In triple-negative breast cancer (TNBC) models, RBMS3 was necessary for maintaining the mesenchymal phenotype and invasion and migration in vitro. Loss of RBMS3 significantly impaired both tumor progression and spontaneous metastasis in vivo. Using a genome-wide approach to interrogate mRNA stability, we found that ectopic expression of RBMS3 upregulates many genes that are resistant to degradation following transcriptional blockade by actinomycin D (ACTD). Specifically, RBMS3 was shown to interact with the mRNA of EMT transcription factor PRRX1 and promote PRRX1 mRNA stability. PRRX1 is required for RBMS3-mediated EMT and is partially sufficient to rescue the effect of RBMS3 knockdown in TNBC cell lines. Together, this study identifies RBMS3 as a novel and common effector of EMT, which could be a promising therapeutic target for TNBC treatment.

Time-Series Clustering Based on the Characterization of Segment Typologies

Time-series clustering is the process of grouping time series with respect to their similarity or characteristics. Previous approaches usually combine a specific distance measure for time series and a standard clustering method. However, these approaches do not take the similarity of the different subsequences of each time series into account, which can be used to better compare the time-series objects of the dataset. In this article, we propose a novel technique of time-series clustering consisting of two clustering stages. In a first step, a least-squares polynomial segmentation procedure is applied to each time series, which is based on a growing window technique that returns different-length segments. Then, all of the segments are projected into the same dimensional space, based on the coefficients of the model that approximates the segment and a set of statistical features. After mapping, a first hierarchical clustering phase is applied to all mapped segments, returning groups of segments for each time series. These clusters are used to represent all time series in the same dimensional space, after defining another specific mapping process. In a second and final clustering stage, all the time-series objects are grouped. We consider internal clustering quality to automatically adjust the main parameter of the algorithm, which is an error threshold for the segmentation. The results obtained on 84 datasets from the UCR Time Series Classification Archive have been compared against three state-of-the-art methods, showing that the performance of this methodology is very promising, especially on larger datasets.

The dynamic, combinatorial cis-regulatory lexicon of epidermal differentiation

Transcription factors bind DNA sequence motif vocabularies in cis-regulatory elements (CREs) to modulate chromatin state and gene expression during cell state transitions. A quantitative understanding of how motif lexicons influence dynamic regulatory activity has been elusive due to the combinatorial nature of the cis-regulatory code. To address this, we undertook multiomic data profiling of chromatin and expression dynamics across epidermal differentiation to identify 40,103 dynamic CREs associated with 3,609 dynamically expressed genes, then applied an interpretable deep-learning framework to model the cis-regulatory logic of chromatin accessibility. This analysis framework identified cooperative DNA sequence rules in dynamic CREs regulating synchronous gene modules with diverse roles in skin differentiation. Massively parallel reporter assay analysis validated temporal dynamics and cooperative cis-regulatory logic. Variants linked to human polygenic skin disease were enriched in these time-dependent combinatorial motif rules. This integrative approach shows the combinatorial cis-regulatory lexicon of epidermal differentiation and represents a general framework for deciphering the organizational principles of the cis-regulatory code of dynamic gene regulation.

Exploring the Meta-regulon of the CRP/FNR Family of Global Transcriptional Regulators in a Partial-Nitritation Anammox Microbiome

Microorganisms must respond to environmental changes to survive, often by controlling transcription initiation. Intermittent aeration during wastewater treatment presents a cyclically changing environment to which microorganisms must react. We used an intermittently aerated bioreactor performing partial nitritation and anammox (PNA) to investigate how the microbiome responds to recurring change. Meta-transcriptomic analysis revealed a dramatic disconnect between the relative DNA abundance and gene expression within the metagenome-assembled genomes (MAGs) of community members, suggesting the importance of transcriptional regulation in this microbiome. To explore how community members responded to cyclic aeration via transcriptional regulation, we searched for homologs of the catabolite repressor protein/fumarate and nitrate reductase regulatory protein (CRP/FNR) family of transcription factors (TFs) within the MAGs. Using phylogenetic analyses, evaluation of sequence conservation in important amino acid residues, and prediction of genes regulated by TFs in the MAGs, we identified homologs of the oxygen-sensing FNR in Nitrosomonas and Rhodocyclaceae, nitrogen-sensing dissimilative nitrate respiration regulator that responds to nitrogen species (DNR) in Rhodocyclaceae, and nitrogen-sensing nitrite and nitric oxide reductase regulator that responds to nitrogen species (NnrR) in Nitrospira MAGs. Our data also predict that CRP/FNR homologs in Ignavibacteria, Flavobacteriales, and Saprospiraceae MAGs sense carbon availability. In addition, a CRP/FNR homolog in a Brocadia MAG was most closely related to CRP TFs known to sense carbon sources in well-studied organisms. However, we predict that in autotrophic Brocadia, this TF most likely regulates a diverse set of functions, including a response to stress during the cyclic aerobic/anoxic conditions. Overall, this analysis allowed us to define a meta-regulon of the PNA microbiome that explains functions and interactions of the most active community members.

IMPORTANCE Microbiomes are important contributors to many ecosystems, including ones where nutrient cycling is stimulated by aeration control. Optimizing cyclic aeration helps reduce energy needs and maximize microbiome performance during wastewater treatment; however, little is known about how most microbial community members respond to these alternating conditions. We defined the meta-regulon of a PNA microbiome by combining existing knowledge of how the CRP/FNR family of bacterial TFs respond to stimuli, with metatranscriptomic analyses to characterize gene expression changes during aeration cycles. Our results indicated that, for some members of the community, prior knowledge is sufficient for high-confidence assignments of TF function, whereas other community members have CRP/FNR TFs for which inferences of function are limited by lack of prior knowledge. This study provides a framework to begin elucidating meta-regulons in microbiomes, where pure cultures are not available for traditional transcriptional regulation studies. Defining the meta-regulon can help in optimizing microbiome performance.