Although markers identified by genome-wide association studies have individually strong statistical significance, their performance in prediction remains limited. Our goal was to use animal breeding genomic prediction models to predict additive genetic contributions for systolic blood pressure (SBP) using whole genome sequencing data with different validation designs. The additive genetic contributions of SBP were estimated via linear mixed model. Rare variants (MAF<0.05) were collapsed through the k-means method to create a "collapsed single-nucleotide polymorphisms." Prediction of the additive genomic contributions of SBP was conducted using genomic Best Linear Unbiased Predictor (GBLUP) and BayesCπ. Estimates of predictive accuracy were compared using common single-nucleotide polymorphisms (SNPs) versus common and collapsed SNPs, and for prediction within and across families. The additive genetic variance of SBP contributed to 18% of the phenotypic variance (h(2) = 0.18). BayesCπ had slightly better prediction accuracies than GBLUP. In both models, within-family predictions had higher accuracies both in the training and testing set than didacross-family design. Collapsing rare variants via the k-means method and adding to the common SNPs did not improve prediction accuracies. The prediction model, including both pedigree and genomic information, achieved a slightly higher accuracy than using either source of information alone. Prediction of genetic contributions to complex traits is feasible using whole genome sequencing and statistical methods borrowed from animal breeding. The relatedness of individuals between the training and testing set strongly affected the performance of prediction models. Methods for inclusion of rare variants in these models need more development.

The salamander has the remarkable ability to regenerate its limb after amputation. Cells at the site of amputation form a blastema and then proliferate and differentiate to regrow the limb. To better understand this process, we performed deep RNA sequencing of the blastema over a time course in the axolotl, a species whose genome has not been sequenced. Using a novel comparative approach to analyzing RNA-seq data, we characterized the transcriptional dynamics of the regenerating axolotl limb with respect to the human gene set. This approach involved de novo assembly of axolotl transcripts, RNA-seq transcript quantification without a reference genome, and transformation of abundances from axolotl contigs to human genes. We found a prominent burst in oncogene expression during the first day and blastemal/limb bud genes peaking at 7 to 14 days. In addition, we found that limb patterning genes, SALL genes, and genes involved in angiogenesis, wound healing, defense/immunity, and bone development are enriched during blastema formation and development. Finally, we identified a category of genes with no prior literature support for limb regeneration that are candidates for further evaluation based on their expression pattern during the regenerative process.

The ability to quantify cellular heterogeneity is a major advantage of single-cell technologies. However, statistical methods often treat cellular heterogeneity as a nuisance. We present a novel method to characterize differences in expression in the presence of distinct expression states within and among biological conditions. We demonstrate that this framework can detect differential expression patterns under a wide range of settings. Compared to existing approaches, this method has higher power to detect subtle differences in gene expression distributions that are more complex than a mean shift, and can characterize those differences. The freely available R package scDD implements the approach.

Oscillatory gene expression is fundamental to development, but technologies for monitoring expression oscillations are limited. We have developed a statistical approach called Oscope to identify and characterize the transcriptional dynamics of oscillating genes in single-cell RNA-seq data from an unsynchronized cell population. Applying Oscope to a number of data sets, we demonstrated its utility and also identified a potential artifact in the Fluidigm C1 platform.

Defects in somitogenesis result in vertebral malformations at birth known as spondylocostal dysostosis (SCDO). Somites are formed with a species-specific periodicity controlled by the "segmentation clock," which comprises a group of oscillatory genes in the presomitic mesoderm. Here, we report that a segmentation clock model derived from human embryonic stem cells shows many hallmarks of the mammalian segmentation clock in vivo, including a dependence on the NOTCH and WNT signaling pathways. The gene expression oscillations are highly synchronized, displaying a periodicity specific to the human clock. Introduction of a point of mutation into HES7, a specific mutation previously associated with clinical SCDO, eliminated clock gene oscillations, successfully reproducing the defects in the segmentation clock. Thus, we provide a model for studying the previously inaccessible human segmentation clock to better understand the mechanisms contributing to congenital skeletal defects.

During vertebrate development, progenitor cells give rise to tissues and organs through a complex choreography that commences at gastrulation. A hallmark event of gastrulation is the formation of the primitive streak, a linear assembly of cells along the anterior-posterior (AP) axis of the developing organism. To examine the primitive streak at a single-cell resolution, we measured the transcriptomes of individual chick cells from the streak or the surrounding tissue (the rest of the area pellucida) in Hamburger-Hamilton stage 4 embryos. The single-cell transcriptomes were then ordered by the statistical method Wave-Crest to deduce both the relative position along the AP axis and the prospective lineage of single cells. The ordered transcriptomes reveal intricate patterns of gene expression along the primitive streak.

Immune checkpoint inhibitors targeting the PD-1/PD-L1 axis lead to durable clinical responses in subsets of cancer patients across multiple indications, including non-small cell lung cancer (NSCLC), urothelial carcinoma (UC) and renal cell carcinoma (RCC). Herein, we complement PD-L1 immunohistochemistry (IHC) and tumor mutation burden (TMB) with RNA-seq in 366 patients to identify unifying and indication-specific molecular profiles that can predict response to checkpoint blockade across these tumor types. Multiple machine learning approaches failed to identify a baseline transcriptional signature highly predictive of response across these indications. Signatures described previously for immune checkpoint inhibitors also failed to validate. At the pathway level, significant heterogeneity is observed between indications, in particular within the PD-L1+ tumors. mUC and NSCLC are molecularly aligned, with cell cycle and DNA damage repair genes associated with response in PD-L1- tumors. At the gene level, the CDK4/6 inhibitor CDKN2A is identified as a significant transcriptional correlate of response, highlighting the association of non-immune pathways to the outcome of checkpoint blockade. This cross-indication analysis reveals molecular heterogeneity between mUC, NSCLC and RCC tumors, suggesting that indication-specific molecular approaches should be prioritized to formulate treatment strategies.

With improvements in next-generation sequencing technologies and reductions in price, ordered RNA-seq experiments are becoming common. Of primary interest in these experiments is identifying genes that are changing over time or space, for example, and then characterizing the specific expression changes. A number of robust statistical methods are available to identify genes showing differential expression among multiple conditions, but most assume conditions are exchangeable and thereby sacrifice power and precision when applied to ordered data.

Reprogramming of somatic cells to induced pluripotent stem cells (iPSCs) represents a profound change in cell fate. Here, we show that combining ascorbic acid (AA) and 2i (MAP kinase and GSK inhibitors) increases the efficiency of reprogramming from fibroblasts and synergistically enhances conversion of partially reprogrammed intermediates to the iPSC state. AA and 2i induce differential transcriptional responses, each leading to the activation of specific pluripotency loci. A unique cohort of pluripotency genes including Esrrb require both stimuli for activation. Temporally, AA-dependent histone demethylase effects are important early, whereas Tet enzyme effects are required throughout the conversion. 2i function could partially be replaced by depletion of components of the epidermal growth factor (EGF) and insulin growth factor pathways, indicating that they act as barriers to reprogramming. Accordingly, reduction in the levels of the EGF receptor gene contributes to the activation of Esrrb. These results provide insight into the rewiring of the pluripotency network at the late stage of reprogramming.

Early life respiratory viral infections and atopic characteristics are significant risk factors for the development of childhood asthma. It is hypothesized that repeated respiratory viral infections might induce structural remodeling by interfering with the normal process of lung maturation; however, the specific molecular processes that underlie these pathological changes are not understood. To investigate the molecular basis for these changes, we used an established Sendai virus infection model in weanling rats to compare the post-infection transcriptomes of an atopic asthma susceptible strain, Brown Norway, and a non-atopic asthma resistant strain, Fischer 344. Specific to this weanling infection model and not described in adult infection models, Sendai virus in the susceptible, but not the resistant strain, results in morphological abnormalities in distal airways that persist into adulthood. Gene expression data from infected and control lungs across five time points indicated that specific features of the immune response following viral infection were heightened and prolonged in lungs from Brown Norway rats compared with Fischer 344 rats. These features included an increase in macrophage cell number and related gene expression, which then transitioned to an increase in mast cell number and related gene expression. In contrast, infected Fischer F344 lungs exhibited more efficient restoration of the airway epithelial morphology, with transient appearance of basal cell pods near distal airways. Together, these findings indicate that the pronounced macrophage and mast cell responses and abnormal re-epithelialization precede the structural defects that developed and persisted in Brown Norway, but not Fischer 344 lungs.

Although immune checkpoint inhibitors (ICIs) are established as effective cancer therapies, overcoming therapeutic resistance remains a critical challenge. Here we identify interleukin 6 (IL-6) as a correlate of poor response to atezolizumab (anti-PD-L1) in large clinical trials of advanced kidney, breast, and bladder cancers. In pre-clinical models, combined blockade of PD-L1 and the IL-6 receptor (IL6R) causes synergistic regression of large established tumors and substantially improves anti-tumor CD8+ cytotoxic T lymphocyte (CTL) responses compared with anti-PD-L1 alone. Circulating CTLs from cancer patients with high plasma IL-6 display a repressed functional profile based on single-cell RNA sequencing, and IL-6-STAT3 signaling inhibits classical cytotoxic differentiation of CTLs in vitro. In tumor-bearing mice, CTL-specific IL6R deficiency is sufficient to improve anti-PD-L1 activity. Thus, based on both clinical and experimental evidence, agents targeting IL-6 signaling are plausible partners for combination with ICIs in cancer patients.

Suncus etruscus is one of the world's smallest mammals, with an average body mass of about 2 grams. The Etruscan shrew's small body is accompanied by a very high energy demand and numerous metabolic adaptations. Here we report a chromosome-level genome assembly using PacBio long read sequencing, 10X Genomics linked short reads, optical mapping, and Hi-C linked reads. The assembly is partially phased, with the 2.472 Gbp primary pseudohaplotype and 1.515 Gbp alternate. We manually curated the primary assembly and identified 22 chromosomes, including X and Y sex chromosomes. The NCBI genome annotation pipeline identified 39,091 genes, 19,819 of them protein-coding. We also identified segmental duplications, inferred GO term annotations, and computed orthologs of human and mouse genes. This reference-quality genome will be an important resource for research on mammalian development, metabolism, and body size control.

The axolotl (Ambystoma mexicanum) has long been the subject of biological research, primarily owing to its outstanding regenerative capabilities. However, the gene expression programs governing its embryonic development are particularly underexplored, especially when compared to other amphibian model species. Therefore, we performed whole transcriptome polyA+ RNA sequencing experiments on 17 stages of embryonic development. As the axolotl genome is unsequenced and its gene annotation is incomplete, we built de novo transcriptome assemblies for each stage and garnered functional annotation by comparing expressed contigs with known genes in other organisms. In evaluating the number of differentially expressed genes over time, we identify three waves of substantial transcriptome upheaval each followed by a period of relative transcriptome stability. The first wave of upheaval is between the one and two cell stage. We show that the number of differentially expressed genes per unit time is higher between the one and two cell stage than it is across the mid-blastula transition (MBT), the period of zygotic genome activation. We use total RNA sequencing to demonstrate that the vast majority of genes with increasing polyA+ signal between the one and two cell stage result from polyadenylation rather than de novo transcription. The first stable phase begins after the two cell stage and continues until the mid-blastula transition, corresponding with the pre-MBT phase of transcriptional quiescence in amphibian development. Following this is a peak of differential gene expression corresponding with the activation of the zygotic genome and a phase of transcriptomic stability from stages 9-11. We observe a third wave of transcriptomic change between stages 11 and 14, followed by a final stable period. The last two stable phases have not been documented in amphibians previously and correspond to times of major morphogenic change in the axolotl embryo: gastrulation and neurulation. These results yield new insights into global gene expression during early stages of amphibian embryogenesis and will help to further develop the axolotl as a model species for developmental and regenerative biology.

DNA methyltransferase 1 (DNMT1) is required for embryogenesis but roles in late forming organ systems including the prostate, which emerges from the urethral epithelium, have not been fully examined. We used a targeted genetic approach involving a Shhcre recombinase to demonstrate requirement of epithelial DNA methyltransferase-1 (Dnmt1) in mouse prostate morphogenesis. Dnmt1 mutant urethral cells exhibit DNA hypomethylation, DNA damage, p53 accumulation and undergo cell cycle arrest and apoptosis. Urethral epithelial cells are disorganized in Dnmt1 mutants, leading to impaired prostate growth and maturation and failed glandular development. We evaluated oriented cell division as a mechanism of bud elongation and widening by demonstrating that mitotic spindle axes typically form parallel or perpendicular to prostatic bud elongation axes. We then deployed a ShhcreERT allele to delete Dnmt1 from a subset of urethral epithelial cells, creating mosaic mutants with which to interrogate the requirement for cell division in specific prostatic bud epithelial populations. DNMT1- cell distribution within prostatic buds is not random as would be expected in a process where DNMT1 was not required. Instead, replication competent DNMT1 + cells primarily accumulate in prostatic bud margins and tips while replication impeded DNMT1- cells accumulate in prostatic bud cores. Together, these results highlight the role of DNMT1 in regulating epithelial bud formation by maintaining cell cycle progression and survival of rapidly dividing urethral epithelial cells, which can be extended to the study of other developing epithelial organs. In addition, our results show that prostatic buds consist of two epithelial cell populations with distinct molecular and functional characteristics that could potentially contribute to specialized lineages in the adult prostate.

While the aging process is central to the pathogenesis of age-dependent diseases, it is poorly understood at the molecular level. We identified a mouse mutant with accelerated aging in the retina as well as pathologies observed in age-dependent retinal diseases, suggesting that the responsible gene regulates retinal aging, and its impairment results in age-dependent disease. We determined that a mutation in the transmembrane 135 (Tmem135) is responsible for these phenotypes. We observed localization of TMEM135 on mitochondria, and imbalance of mitochondrial fission and fusion in mutant Tmem135 as well as Tmem135 overexpressing cells, indicating that TMEM135 is involved in the regulation of mitochondrial dynamics. Additionally, mutant retina showed higher sensitivity to oxidative stress. These results suggest that the regulation of mitochondrial dynamics through TMEM135 is critical for protection from environmental stress and controlling the progression of retinal aging. Our study identified TMEM135 as a critical link between aging and age-dependent diseases.

Genome-wide association studies have revealed that many low-penetrance breast cancer susceptibility loci are located in non-protein coding genomic regions; however, few have been characterized. In a comparative genetics approach to model such loci in a rat breast cancer model, we previously identified the mammary carcinoma susceptibility locus Mcs1a. We now localize Mcs1a to a critical interval (277 Kb) within a gene desert. Mcs1a reduces mammary carcinoma multiplicity by 50% and acts in a mammary cell-autonomous manner. We developed a megadeletion mouse model, which lacks 535 Kb of sequence containing the Mcs1a ortholog. Global gene expression analysis by RNA-seq revealed that in the mouse mammary gland, the orphan nuclear receptor gene Nr2f1/Coup-tf1 is regulated by Mcs1a. In resistant Mcs1a congenic rats, as compared with susceptible congenic control rats, we found Nr2f1 transcript levels to be elevated in mammary gland, epithelial cells, and carcinoma samples. Chromatin looping over ∼820 Kb of sequence from the Nr2f1 promoter to a strongly conserved element within the Mcs1a critical interval was identified. This element contains a 14 bp indel polymorphism that affects a human-rat-mouse conserved COUP-TF binding motif and is a functional Mcs1a candidate. In both the rat and mouse models, higher Nr2f1 transcript levels are associated with higher abundance of luminal mammary epithelial cells. In both the mouse mammary gland and a human breast cancer global gene expression data set, we found Nr2f1 transcript levels to be strongly anti-correlated to a gene cluster enriched in cell cycle-related genes. We queried 12 large publicly available human breast cancer gene expression studies and found that the median NR2F1 transcript level is consistently lower in 'triple-negative' (ER-PR-HER2-) breast cancers as compared with 'receptor-positive' breast cancers. Our data suggest that the non-protein coding locus Mcs1a regulates Nr2f1, which is a candidate modifier of differentiation, proliferation, and mammary cancer risk.

The normalization of RNA-seq data is essential for accurate downstream inference, but the assumptions upon which most normalization methods are based are not applicable in the single-cell setting. Consequently, applying existing normalization methods to single-cell RNA-seq data introduces artifacts that bias downstream analyses. To address this, we introduce SCnorm for accurate and efficient normalization of single-cell RNA-seq data.

Integrated multi-omics evaluation of 823 tumors from advanced renal cell carcinoma (RCC) patients identifies molecular subsets associated with differential clinical outcomes to angiogenesis blockade alone or with a checkpoint inhibitor. Unsupervised transcriptomic analysis reveals seven molecular subsets with distinct angiogenesis, immune, cell-cycle, metabolism, and stromal programs. While sunitinib and atezolizumab + bevacizumab are effective in subsets with high angiogenesis, atezolizumab + bevacizumab improves clinical benefit in tumors with high T-effector and/or cell-cycle transcription. Somatic mutations in PBRM1 and KDM5C associate with high angiogenesis and AMPK/fatty acid oxidation gene expression, while CDKN2A/B and TP53 alterations associate with increased cell-cycle and anabolic metabolism. Sarcomatoid tumors exhibit lower prevalence of PBRM1 mutations and angiogenesis markers, frequent CDKN2A/B alterations, and increased PD-L1 expression. These findings can be applied to molecularly stratify patients, explain improved outcomes of sarcomatoid tumors to checkpoint blockade versus antiangiogenics alone, and develop personalized therapies in RCC and other indications.

Cellular gene expression changes throughout a dynamic biological process, such as differentiation. Pseudotimes estimate cells' progress along a dynamic process based on their individual gene expression states. Ordering the expression data by pseudotime provides information about the underlying regulator-gene interactions. Because the pseudotime distribution is not uniform, many standard mathematical methods are inapplicable for analyzing the ordered gene expression states. Here we present single-cell inference of networks using Granger ensembles (SINGE), an algorithm for gene regulatory network inference from ordered single-cell gene expression data. SINGE uses kernel-based Granger causality regression to smooth irregular pseudotimes and missing expression values. It aggregates predictions from an ensemble of regression analyses to compile a ranked list of candidate interactions between transcriptional regulators and target genes. In two mouse embryonic stem cell differentiation datasets, SINGE outperforms other contemporary algorithms. However, a more detailed examination reveals caveats about poor performance for individual regulators and uninformative pseudotimes.

Sus scrofa domesticus (pig) has served as a superb large mammalian model for biomedical studies because of its comparable physiology and organ size to humans. The derivation of transgene-free porcine induced pluripotent stem cells (PiPSCs) will, therefore, benefit the development of porcine-specific models for regenerative biology and its medical applications. In the past, this effort has been hampered by a lack of understanding of the signaling milieu that stabilizes the porcine pluripotent state in vitro. Here, we report that transgene-free PiPSCs can be efficiently derived from porcine fibroblasts by episomal vectors along with microRNA-302/367 using optimized protocols tailored for this species. PiPSCs can be differentiated into derivatives representing the primary germ layers in vitro and can form teratomas in immunocompromised mice. Furthermore, the transgene-free PiPSCs preserve intrinsic species-specific developmental timing in culture, known as developmental allochrony. This is demonstrated by establishing a porcine in vitro segmentation clock model that, for the first time, displays a specific periodicity at ∼3.7 h, a timescale recapitulating in vivo porcine somitogenesis. We conclude that the transgene-free PiPSCs can serve as a powerful tool for modeling development and disease and developing transplantation strategies. We also anticipate that they will provide insights into conserved and unique features on the regulations of mammalian pluripotency and developmental timing mechanisms.