Optimal scaling of digital transcriptomes

Statistical normalization methods in microbiome data with application to microbiome cancer research

Mounting evidence has shown that gut microbiome is associated with various cancers, including gastrointestinal (GI) tract and non-GI tract cancers. But microbiome data have unique characteristics and pose major challenges when using standard statistical methods causing results to be invalid or misleading. Thus, to analyze microbiome data, it not only needs appropriate statistical methods, but also requires microbiome data to be normalized prior to statistical analysis. Here, we first describe the unique characteristics of microbiome data and the challenges in analyzing them (Section 2). Then, we provide an overall review on the available normalization methods of 16S rRNA and shotgun metagenomic data along with examples of their applications in microbiome cancer research (Section 3). In Section 4, we comprehensively investigate how the normalization methods of 16S rRNA and shotgun metagenomic data are evaluated. Finally, we summarize and conclude with remarks on statistical normalization methods (Section 5). Altogether, this review aims to provide a broad and comprehensive view and remarks on the promises and challenges of the statistical normalization methods in microbiome data with microbiome cancer research examples.

FDG uptake reflects an immune-enriched subtype of thyroid cancer: Clinical implications of imaging-based molecular characterization

INTRODUCTION

Iodine and FDG uptakes have been established as methods to define the biological properties of thyroid cancer. As various cells in the tumor microenvironment (TME) affect tumor metabolism, we investigated the association between glucose metabolism in thyroid cancer and the TME using transcriptomic analyses.

METHODS

We used F-18 FDG PET and RNA sequencing data of thyroid cancer to find associations between TME cell types and glucose metabolism. In addition, publicly available single-cell RNA sequencing data of papillary thyroid cancer was used to investigate glucose metabolism in cell types of the TME. The correlations between the FDG uptake and biological properties of the TME, including glucose metabolism and tumor differentiation score (TDS) were evaluated. Estimation of the proportions of immune and cancer cells (EPIC) was performed. The biological properties of each cell type were also assessed in the single-cell RNA sequencing data.

RESULTS

FDG uptake showed a positive correlation with the enrichment score of macrophages and glycolysis activity. In single-cell RNA sequencing, immune cells had both high glucose transporters (GLUTs) and glycolysis signatures, while thyrocytes including cancer cells showed relatively low GLUTs and glycolysis signatures, suggesting that FDG uptake mainly occurred in immune cells of the TME. Moreover, the high GLUTs of myeloid cells were negatively associated with TDS.

CONCLUSIONS

Our findings suggest that thyroid cancer with high FDG uptake can be mediated by enriched immune cells of the TME. We suggest that FDG uptake in thyroid cancer could be a marker for the immune-rich type and provide clinical implications for treatment stratification.

Clustering by antigen-presenting genes reveals immune landscapes and predicts response to checkpoint immunotherapy

Immune checkpoint blockade (ICB) has demonstrated efficacy by reinvigorating immune cytotoxicity against tumors. However, the mechanisms underlying how ICB induces responses in a subset of patients remain unclear. Using bulk and single-cell transcriptomic cohorts of melanoma patients receiving ICB, we proposed a clustering model based on the expression of an antigen-presenting machinery (APM) signature consisting of 23 genes in a forward-selection manner. We characterized four APM clusters associated with distinct immune characteristics, cancer hallmarks, and patient prognosis in melanoma. The model predicts differential regulation of APM genes during ICB, which shaped ICB responsiveness. Surprisingly, while immunogenically hot tumors with high baseline APM expression prior to treatment are correlated with a better response to ICB than cold tumors with low APM expression, a subset of hot tumors with the highest pre-ICB APM expression fail to upregulate APM expression during treatment. In addition, they undergo immunoediting and display infiltration of exhausted T cells. In comparison, tumors associated with the best patient prognosis demonstrate significant APM upregulation and immune infiltration following ICB. They also show infiltration of tissue-resident memory T cells, shaping prolonged antitumor immunity. Using only pre-treatment transcriptomic data, our model predicts the dynamic APM-mediated tumor-immune interactions in response to ICB and provides insights into the immune escape mechanisms in hot tumors that compromise the ICB efficacy. We highlight the prognostic value of APM expression in predicting immune response in chronic diseases.

cdev: a ground-truth based measure to evaluate RNA-seq normalization performance

Normalization of RNA-seq data has been an active area of research since the problem was first recognized a decade ago. Despite the active development of new normalizers, their performance measures have been given little attention. To evaluate normalizers, researchers have been relying on ad hoc measures, most of which are either qualitative, potentially biased, or easily confounded by parametric choices of downstream analysis. We propose a metric called condition-number based deviation, or cdev, to quantify normalization success. cdev measures how much an expression matrix differs from another. If a ground truth normalization is given, cdev can then be used to evaluate the performance of normalizers. To establish experimental ground truth, we compiled an extensive set of public RNA-seq assays with external spike-ins. This data collection, together with cdev, provides a valuable toolset for benchmarking new and existing normalization methods.

A peptide encoded within a 5' untranslated region promotes pain sensitization in mice

ABSTRACT

Translational regulation permeates neuronal function. Nociceptors are sensory neurons responsible for the detection of harmful stimuli. Changes in their activity, termed plasticity, are intimately linked to the persistence of pain. Although inhibitors of protein synthesis robustly attenuate pain-associated behavior, the underlying targets that support plasticity are largely unknown. Here, we examine the contribution of protein synthesis in regions of RNA annotated as noncoding. Based on analyses of previously reported ribosome profiling data, we provide evidence for widespread translation in noncoding transcripts and regulatory regions of mRNAs. We identify an increase in ribosome occupancy in the 5' untranslated regions of the calcitonin gene-related peptide (CGRP/Calca). We validate the existence of an upstream open reading frame (uORF) using a series of reporter assays. Fusion of the uORF to a luciferase reporter revealed active translation in dorsal root ganglion neurons after nucleofection. Injection of the peptide corresponding to the calcitonin gene-related peptide-encoded uORF resulted in pain-associated behavioral responses in vivo and nociceptor sensitization in vitro. An inhibitor of heterotrimeric G protein signaling blocks both effects. Collectively, the data suggest pervasive translation in regions of the transcriptome annotated as noncoding in dorsal root ganglion neurons and identify a specific uORF-encoded peptide that promotes pain sensitization through GPCR signaling.

Modulating innate immune activation states impacts the efficacy of specific Aβ immunotherapy

INTRODUCTION

Passive immunotherapies targeting Aβ continue to be evaluated as Alzheimer's disease (AD) therapeutics, but there remains debate over the mechanisms by which these immunotherapies work. Besides the amount of preexisting Aβ deposition and the type of deposit (compact or diffuse), there is little data concerning what factors, independent of those intrinsic to the antibody, might influence efficacy. Here we (i) explored how constitutive priming of the underlying innate activation states by Il10 and Il6 might influence passive Aβ immunotherapy and (ii) evaluated transcriptomic data generated in the AMP-AD initiative to inform how these two cytokines and their receptors' mRNA levels are altered in human AD and an APP mouse model.

METHODS

rAAV2/1 encoding EGFP, Il6 or Il10 were delivered by somatic brain transgenesis to neonatal (P0) TgCRND8 APP mice. Then, at 2 months of age, the mice were treated bi-weekly with a high-affinity anti-Aβ1-16 mAb5 monoclonal antibody or control mouse IgG until 6 months of age. rAAV mediated transgene expression, amyloid accumulation, Aβ levels and gliosis were assessed. Extensive transcriptomic data was used to evaluate the mRNA expression levels of IL10 and IL6 and their receptors in the postmortem human AD temporal cortex and in the brains of TgCRND8 mice, the later at multiple ages.

RESULTS

Priming TgCRND8 mice with Il10 increases Aβ loads and blocks efficacy of subsequent mAb5 passive immunotherapy, whereas priming with Il6 priming reduces Aβ loads by itself and subsequent Aβ immunotherapy shows only a slightly additive effect. Transcriptomic data shows that (i) there are significant increases in the mRNA levels of Il6 and Il10 receptors in the TgCRND8 mouse model and temporal cortex of humans with AD and (ii) there is a great deal of variance in individual mouse brain and the human temporal cortex of these interleukins and their receptors.

CONCLUSIONS

The underlying immune activation state can markedly affect the efficacy of passive Aβ immunotherapy. These results have important implications for ongoing human AD immunotherapy trials, as they indicate that underlying immune activation states within the brain, which may be highly variable, may influence the ability for passive immunotherapy to alter Aβ deposition.

Measurement of Organ-Specific and Acute-Phase Blood Protein Levels in Early Lyme Disease

Lyme disease results from infection of humans with the spirochete Borrelia burgdorferi. The first and most common clinical manifestation is the circular, inflamed skin lesion referred to as erythema migrans; later manifestations result from infections of other body sites. Laboratory diagnosis of Lyme disease can be challenging in patients with erythema migrans because of the time delay in the development of specific diagnostic antibodies against Borrelia. Reliable blood biomarkers for the early diagnosis of Lyme disease in patients with erythema migrans are needed. Here, we performed selected reaction monitoring, a targeted mass spectrometry-based approach, to measure selected proteins that (1) are known to be predominantly expressed in one organ (i.e., organ-specific blood proteins) and whose blood concentrations may change as a result of Lyme disease, or (2) are involved in acute immune responses. In a longitudinal cohort of 40 Lyme disease patients and 20 healthy controls, we identified 10 proteins with significantly altered serum levels in patients at the time of diagnosis, and we also developed a 10-protein panel identified through multivariate analysis. In an independent cohort of patients with erythema migrans, six of these proteins, APOA4, C9, CRP, CST6, PGLYRP2, and S100A9, were confirmed to show significantly altered serum levels in patients at time of presentation. Nine of the 10 proteins from the multivariate panel were also verified in the second cohort. These proteins, primarily innate immune response proteins or proteins specific to liver, skin, or white blood cells, may serve as candidate blood biomarkers requiring further validation to aid in the laboratory diagnosis of early Lyme disease.

New steps in mucilage biosynthesis revealed by analysis of the transcriptome of the UDP-rhamnose/UDP-galactose transporter 2 mutant

Upon imbibition, epidermal cells of Arabidopsis thaliana seeds release a mucilage formed mostly by pectic polysaccharides. The Arabidopsis mucilage is composed mainly of unbranched rhamnogalacturonan-I (RG-I), with low amounts of cellulose, homogalacturonan, and traces of xylan, xyloglucan, galactoglucomannan, and galactan. The pectin-rich composition of the mucilage and their simple extractability makes this structure a good candidate to study the biosynthesis of pectic polysaccharides and their modification. Here, we characterize the mucilage phenotype of a mutant in the UDP-rhamnose/galactose transporter 2 (URGT2), which exhibits a reduction in RG-I and also shows pleiotropic changes, suggesting the existence of compensation mechanisms triggered by the lack of URGT2. To gain an insight into the possible compensation mechanisms activated in the mutant, we performed a transcriptome analysis of developing seeds using RNA sequencing (RNA-seq). The results showed a significant misregulation of 3149 genes, 37 of them (out of the 75 genes described to date) encoding genes proposed to be involved in mucilage biosynthesis and/or its modification. The changes observed in urgt2 included the up-regulation of UAFT2, a UDP-arabinofuranose transporter, and UUAT3, a paralog of the UDP-uronic acid transporter UUAT1, suggesting that they play a role in mucilage biosynthesis. Mutants in both genes showed changes in mucilage composition and structure, confirming their participation in mucilage biosynthesis. Our results suggest that plants lacking a UDP-rhamnose/galactose transporter undergo important changes in gene expression, probably to compensate modifications in the plant cell wall due to the lack of a gene involved in its biosynthesis.

Differences between Dorsal Root and Trigeminal Ganglion Nociceptors in Mice Revealed by Translational Profiling

Nociceptors located in the trigeminal ganglion (TG) and DRG are the primary sensors of damaging or potentially damaging stimuli for the head and body, respectively, and are key drivers of chronic pain states. While nociceptors in these two tissues show a high degree of functional similarity, there are important differences in their development lineages, their functional connections to the CNS, and recent genome-wide analyses of gene expression suggest that they possess some unique genomic signatures. Here, we used translating ribosome affinity purification to comprehensively characterize and compare mRNA translation in Scn10a-positive nociceptors in the TG and DRG of male and female mice. This unbiased method independently confirms several findings of differences between TG and DRG nociceptors described in the literature but also suggests preferential utilization of key signaling pathways. Most prominently, we provide evidence that translational efficiency in mechanistic target of rapamycin (mTOR)-related genes is higher in the TG compared with DRG, whereas several genes associated with the negative regulator of mTOR, AMP-activated protein kinase, have higher translational efficiency in DRG nociceptors. Using capsaicin as a sensitizing stimulus, we show that behavioral responses are greater in the TG region and this effect is completely reversible with mTOR inhibition. These findings have implications for the relative capacity of these nociceptors to be sensitized upon injury. Together, our data provide a comprehensive, comparative view of transcriptome and translatome activity in TG and DRG nociceptors that enhances our understanding of nociceptor biology.SIGNIFICANCE STATEMENT The DRG and trigeminal ganglion (TG) provide sensory information from the body and head, respectively. Nociceptors in these tissues are critical first neurons in the pain pathway. Injury to peripheral neurons in these tissues can cause chronic pain. Interestingly, clinical and preclinical findings support the conclusion that injury to TG neurons is more likely to cause chronic pain and chronic pain in the TG area is more intense and more difficult to treat. We used translating ribosome affinity purification technology to gain new insight into potential differences in the translatomes of DRG and TG neurons. Our findings demonstrate previously unrecognized differences between TG and DRG nociceptors that provide new insight into how injury may differentially drive plasticity states in nociceptors in these two tissues.

NormExpression: An R Package to Normalize Gene Expression Data Using Evaluated Methods

Data normalization is a crucial step in the gene expression analysis as it ensures the validity of its downstream analyses. Although many metrics have been designed to evaluate the existing normalization methods, different metrics or different datasets by the same metric yield inconsistent results, particularly for the single-cell RNA sequencing (scRNA-seq) data. The worst situations could be that one method evaluated as the best by one metric is evaluated as the poorest by another metric, or one method evaluated as the best using one dataset is evaluated as the poorest using another dataset. Here raises an open question: principles need to be established to guide the evaluation of normalization methods. In this study, we propose a principle that one normalization method evaluated as the best by one metric should also be evaluated as the best by another metric (the consistency of metrics) and one method evaluated as the best using scRNA-seq data should also be evaluated as the best using bulk RNA-seq data or microarray data (the consistency of datasets). Then, we designed a new metric named Area Under normalized CV threshold Curve (AUCVC) and applied it with another metric mSCC to evaluate 14 commonly used normalization methods using both scRNA-seq data and bulk RNA-seq data, satisfying the consistency of metrics and the consistency of datasets. Our findings paved the way to guide future studies in the normalization of gene expression data with its evaluation. The raw gene expression data, normalization methods, and evaluation metrics used in this study have been included in an R package named NormExpression. NormExpression provides a framework and a fast and simple way for researchers to select the best method for the normalization of their gene expression data based on the evaluation of different methods (particularly some data-driven methods or their own methods) in the principle of the consistency of metrics and the consistency of datasets.

Comparative transcriptome profiling of the human and mouse dorsal root ganglia: an RNA-seq-based resource for pain and sensory neuroscience research

Molecular neurobiological insight into human nervous tissues is needed to generate next-generation therapeutics for neurological disorders such as chronic pain. We obtained human dorsal root ganglia (hDRG) samples from organ donors and performed RNA-sequencing (RNA-seq) to study the hDRG transcriptional landscape, systematically comparing it with publicly available data from a variety of human and orthologous mouse tissues, including mouse DRG (mDRG). We characterized the hDRG transcriptional profile in terms of tissue-restricted gene coexpression patterns and putative transcriptional regulators, and formulated an information-theoretic framework to quantify DRG enrichment. Relevant gene families and pathways were also analyzed, including transcription factors, G-protein-coupled receptors, and ion channels. Our analyses reveal an hDRG-enriched protein-coding gene set (∼140), some of which have not been described in the context of DRG or pain signaling. Most of these show conserved enrichment in mDRG and were mined for known drug-gene product interactions. Conserved enrichment of the vast majority of transcription factors suggests that the mDRG is a faithful model system for studying hDRG, because of evolutionarily conserved regulatory programs. Comparison of hDRG and tibial nerve transcriptomes suggests trafficking of neuronal mRNA to axons in adult hDRG, and are consistent with studies of axonal transport in rodent sensory neurons. We present our work as an online, searchable repository (https://www.utdallas.edu/bbs/painneurosciencelab/sensoryomics/drgtxome), creating a valuable resource for the community. Our analyses provide insight into DRG biology for guiding development of novel therapeutics and a blueprint for cross-species transcriptomic analyses.

Nociceptor Translational Profiling Reveals the Ragulator-Rag GTPase Complex as a Critical Generator of Neuropathic Pain

Nociceptors, sensory neurons in the DRG that detect damaging or potentially damaging stimuli, are key drivers of neuropathic pain. Injury to these neurons causes activation of translation regulation signaling, including the mechanistic target of rapamycin complex 1 (mTORC1) and mitogen-activated protein kinase interacting kinase (MNK) eukaryotic initiation factor (eIF) 4E pathways. This is a mechanism driving changes in excitability of nociceptors that is critical for the generation of chronic pain states; however, the mRNAs that are translated to lead to this plasticity have not been elucidated. To address this gap in knowledge, we used translating ribosome affinity purification in male and female mice to comprehensively characterize mRNA translation in Scn10a-positive nociceptors in chemotherapy-induced neuropathic pain (CIPN) caused by paclitaxel treatment. This unbiased method creates a new resource for the field, confirms many findings in the CIPN literature and also find extensive evidence for new target mechanisms that may cause CIPN. We provide evidence that an underlying mechanism of CIPN is sustained mTORC1 activation driven by MNK1-eIF4E signaling. RagA, a GTPase controlling mTORC1 activity, is identified as a novel target of MNK1-eIF4E signaling. This demonstrates a novel translation regulation signaling circuit wherein MNK1-eIF4E activity drives mTORC1 via control of RagA translation. CIPN and RagA translation are strongly attenuated by genetic ablation of eIF4E phosphorylation, MNK1 elimination or treatment with the MNK inhibitor eFT508. We identify a novel translational circuit for the genesis of neuropathic pain caused by chemotherapy with important implications for therapeutics.SIGNIFICANCE STATEMENT Neuropathic pain affects up to 10% of the population, but its underlying mechanisms are incompletely understood, leading to poor treatment outcomes. We used translating ribosome affinity purification technology to create a comprehensive translational profile of DRG nociceptors in naive mice and at the peak of neuropathic pain induced by paclitaxel treatment. We reveal new insight into how mechanistic target of rapamycin complex 1 is activated in neuropathic pain pointing to a key role of MNK1-eIF4E-mediated translation of a complex of mRNAs that control mechanistic target of rapamycin complex 1 signaling at the surface of the lysosome. We validate this finding using genetic and pharmacological techniques. Our work strongly suggests that MNK1-eIF4E signaling drives CIPN and that a drug in human clinical trials, eFT508, may be a new therapeutic for neuropathic pain.

Expression QTL (eQTLs) Analyses Reveal Candidate Genes Associated With Fruit Flesh Softening Rate in Peach [Prunus persica (L.) Batsch]

Significant differences in softening rate have been reported between melting flesh in peach and nectarine varieties. This trait seems to be controlled by several genes. We aimed to identify candidate genes involved in fruit softening rate by integrating quantitative trait loci (QTL) and expression QTL (eQTL) analyses, comparing siblings with contrasting softening rates. We used a segregating population derived from nectarine cv. 'Venus' selfing, which was phenotyped for softening rate during three seasons. Six siblings with high (HSR) and six with low softening rate (LSR) were sequenced using RNA-Seq. A group of 5,041 differentially expressed genes was identified. Also, we found a QTL with a LOD (logarithm of odds) score of 9.7 on LG4 in all analyzed seasons. Furthermore, we detected 1,062 eQTLs, of which 133 were found co-localizing with the identified QTL. Gene Ontology (GO) analysis showed 'Response to auxin' as one the main over-represented categories. Our findings suggest over-expression of auxin biosynthetic related genes in the HSR group, which implies a higher expression and/or accumulation of auxin, thereby triggering fast softening. Conversely, the LSR phenotype might be explained by an altered auxin-homeostasis associated with low auxin levels. This work will contribute to unraveling the genetic mechanisms responsible for the softening rate in peaches and nectarines and lead to the development of molecular markers.

A map of gene expression in neutrophil-like cell lines

BACKGROUND

Human neutrophils are central players in innate immunity, a major component of inflammatory responses, and a leading model for cell motility and chemotaxis. However, primary neutrophils are short-lived, limiting their experimental usefulness in the laboratory. Thus, human myeloid cell lines have been characterized for their ability to undergo neutrophil-like differentiation in vitro. The HL-60 cell line and its PLB-985 sub-line are commonly used to model human neutrophil behavior, but how closely gene expression in differentiated cells resembles that of primary neutrophils has remained unclear.

RESULTS

In this study, we compared the effectiveness of differentiation protocols and used RNA sequencing (RNA-seq) to compare the transcriptomes of HL-60 and PLB-985 cells with published data for human and mouse primary neutrophils. Among commonly used differentiation protocols for neutrophil-like cell lines, addition of dimethyl sulfoxide (DMSO) gave the best combination of cell viability and expression of markers for differentiation. However, combining DMSO with the serum-free-supplement Nutridoma resulted in increased chemotactic response, phagocytic activity, oxidative burst and cell surface expression of the neutrophil markers FPR1 and CD11b without a cost in viability. RNA-seq analysis of HL-60 and PLB-985 cells before and after differentiation showed that differentiation broadly increases the similarity in gene expression between the cell lines and primary neutrophils. Furthermore, the gene expression pattern of the differentiated cell lines correlated slightly better with that of human neutrophils than the mouse neutrophil pattern did. Finally, we created a publicly available gene expression database that is searchable by gene name and protein domain content, where users can compare gene expression in HL-60, PLB-985 and primary human and mouse neutrophils.

CONCLUSIONS

Our study verifies that a DMSO-based differentiation protocol for HL-60 and PLB-985 cell lines gives superior differentiation and cell viability relative to other common protocols, and indicates that addition of Nutridoma may be preferable for studies of chemotaxis, phagocytosis, or oxidative burst. Our neutrophil gene expression database will be a valuable tool to identify similarities and differences in gene expression between the cell lines and primary neutrophils, to compare expression levels for genes of interest, and to improve the design of tools for genetic perturbations.

Genome-wide characterization of mammalian promoters with distal enhancer functions

Gene expression in mammals is precisely regulated by the combination of promoters and gene-distal regulatory regions, known as enhancers. Several studies have suggested that some promoters might have enhancer functions. However, the extent of this type of promoters and whether they actually function to regulate the expression of distal genes have remained elusive. Here, by exploiting a high-throughput enhancer reporter assay, we unravel a set of mammalian promoters displaying enhancer activity. These promoters have distinct genomic and epigenomic features and frequently interact with other gene promoters. Extensive CRISPR-Cas9 genomic manipulation demonstrated the involvement of these promoters in the cis regulation of expression of distal genes in their natural loci. Our results have important implications for the understanding of complex gene regulation in normal development and disease.

cgCorrect: a method to correct for confounding cell-cell variation due to cell growth in single-cell transcriptomics

Accessing gene expression at a single-cell level has unraveled often large heterogeneity among seemingly homogeneous cells, which remains obscured when using traditional population-based approaches. The computational analysis of single-cell transcriptomics data, however, still imposes unresolved challenges with respect to normalization, visualization and modeling the data. One such issue is differences in cell size, which introduce additional variability into the data and for which appropriate normalization techniques are needed. Otherwise, these differences in cell size may obscure genuine heterogeneities among cell populations and lead to overdispersed steady-state distributions of mRNA transcript numbers. We present cgCorrect, a statistical framework to correct for differences in cell size that are due to cell growth in single-cell transcriptomics data. We derive the probability for the cell-growth-corrected mRNA transcript number given the measured, cell size-dependent mRNA transcript number, based on the assumption that the average number of transcripts in a cell increases proportionally to the cell's volume during the cell cycle. cgCorrect can be used for both data normalization and to analyze the steady-state distributions used to infer the gene expression mechanism. We demonstrate its applicability on both simulated data and single-cell quantitative real-time polymerase chain reaction (PCR) data from mouse blood stem and progenitor cells (and to quantitative single-cell RNA-sequencing data obtained from mouse embryonic stem cells). We show that correcting for differences in cell size affects the interpretation of the data obtained by typically performed computational analysis.

Neural transcriptome reveals molecular mechanisms for temporal control of vocalization across multiple timescales

Background

Vocalization is a prominent social behavior among vertebrates, including in the midshipman fish, an established model for elucidating the neural basis of acoustic communication. Courtship vocalizations produced by territorial males are essential for reproductive success, vary over daily and seasonal cycles, and last up to hours per call. Vocalizations rely upon extreme synchrony and millisecond precision in the firing of a homogeneous population of motoneurons, the vocal motor nucleus (VMN). Although studies have identified neural mechanisms driving rapid, precise, and stable neuronal firing over long periods of calling, little is known about underlying genetic/molecular mechanisms.

Results

We used RNA sequencing-based transcriptome analyses to compare patterns of gene expression in VMN to the surrounding hindbrain across three daily and seasonal time points of high and low sound production to identify candidate genes that underlie VMN’s intrinsic and network neuronal properties. Results from gene ontology enrichment, enzyme pathway mapping, and gene category-wide expression levels highlighted the importance of cellular respiration in VMN function, consistent with the high energetic demands of sustained vocal behavior. Functionally important candidate genes upregulated in the VMN, including at time points corresponding to high natural vocal activity, encode ion channels and neurotransmitter receptors, hormone receptors and biosynthetic enzymes, neuromodulators, aerobic respiration enzymes, and antioxidants. Quantitative PCR and RNA-seq expression levels for 28 genes were significantly correlated. Many candidate gene products regulate mechanisms of neuronal excitability, including those previously identified in VMN motoneurons, as well as novel ones that remain to be investigated. Supporting evidence from previous studies in midshipman strongly validate the value of transcriptomic analyses for linking genes to neural characters that drive behavior.

Conclusions

Transcriptome analyses highlighted a suite of molecular mechanisms that regulate vocalization over behaviorally relevant timescales, spanning milliseconds to hours and seasons. To our knowledge, this is the first comprehensive characterization of gene expression in a dedicated vocal motor nucleus. Candidate genes identified here may belong to a conserved genetic toolkit for vocal motoneurons facing similar energetic and neurophysiological demands.

Electronic supplementary material

The online version of this article (doi:10.1186/s12864-015-1577-2) contains supplementary material, which is available to authorized users.

Identification of copy number variants in whole-genome data using Reference Coverage Profiles

The identification of DNA copy numbers from short-read sequencing data remains a challenge for both technical and algorithmic reasons. The raw data for these analyses are measured in tens to hundreds of gigabytes per genome; transmitting, storing, and analyzing such large files is cumbersome, particularly for methods that analyze several samples simultaneously. We developed a very efficient representation of depth of coverage (150–1000× compression) that enables such analyses. Current methods for analyzing variants in whole-genome sequencing (WGS) data frequently miss copy number variants (CNVs), particularly hemizygous deletions in the 1–100 kb range. To fill this gap, we developed a method to identify CNVs in individual genomes, based on comparison to joint profiles pre-computed from a large set of genomes. We analyzed depth of coverage in over 6000 high quality (>40×) genomes. The depth of coverage has strong sequence-specific fluctuations only partially explained by global parameters like %GC. To account for these fluctuations, we constructed multi-genome profiles representing the observed or inferred diploid depth of coverage at each position along the genome. These Reference Coverage Profiles (RCPs) take into account the diverse technologies and pipeline versions used. Normalization of the scaled coverage to the RCP followed by hidden Markov model (HMM) segmentation enables efficient detection of CNVs and large deletions in individual genomes. Use of pre-computed multi-genome coverage profiles improves our ability to analyze each individual genome. We make available RCPs and tools for performing these analyses on personal genomes. We expect the increased sensitivity and specificity for individual genome analysis to be critical for achieving clinical-grade genome interpretation.