RNA-Seq Reveals Pathways Responsible for Meat Quality Characteristic Differences between Two Yunnan Indigenous Chicken Breeds and Commercial Broilers.

Poultry is a source of meat that is in great demand in the world. The quality of meat is an imperative point for shoppers. To explore the genes controlling meat quality characteristics, the growth and meat quality traits and muscle transcriptome of two indigenous Yunnan chicken breeds, Wuding chickens (WDs) and Daweishan mini chickens (MCs), were compared with Cobb broilers (CBs). The growth and meat quality characteristics of these two indigenous breeds were found to differ from CB. In particular, the crude fat (CF), inosine monophosphate content, amino acid (AA), and total fatty acid (TFA) content of WDs were significantly higher than those of CBs and MCs. In addition, it was found that MC pectoralis had 420 differentially expressed genes (DEGs) relative to CBs, and WDs had 217 DEGs relative to CBs. Among them, 105 DEGs were shared. The results of 10 selected genes were also confirmed by qPCR. The differentially expressed genes were six enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) biological pathways including lysosomes, phagosomes, PPAR signaling pathways, cell adhesion molecules, cytokine-cytokine receptor interaction, and phagosome sphingolipid metabolism. Interestingly, four genes (LPL, GK, SCD, and FABP7) in the PPAR signal pathway related to fatty acid (FA) metabolism were elevated in WD muscles, which may account for higher CF, inosine monophosphate content, and AA and FA contents, key factors affecting meat quality. This work laid the foundation for improving the meat quality of Yunnan indigenous chickens, especially WD. In future molecular breeding, the genes in this study can be used as molecular screening markers and applied to the molecular breeding of chicken quality characteristics.

Transcriptomic Analysis of the Spatiotemporal Axis of Oogenesis and Fertilization in C. elegans.

The oocyte germline of the C. elegans hermaphrodite presents a unique model to study the formation of oocytes. However, the size of the model animal and difficulties in retrieval of specific stages of the germline have obviated closer systematic studies of this process throughout the years. Here, we present a transcriptomic level analysis into the oogenesis of C. elegans hermaphrodites. We dissected a hermaphrodite gonad into seven sections corresponding to the mitotic distal region, the pachytene, the diplotene, the early diakinesis region and the 3 most proximal oocytes, and deeply sequenced the transcriptome of each of them along with that of the fertilized egg using a single-cell RNA-seq protocol. We identified specific gene expression events as well as gene splicing events in finer detail along the oocyte germline and provided novel insights into underlying mechanisms of the oogenesis process. Furthermore, through careful review of relevant research literature coupled with patterns observed in our analysis, we attempt to delineate transcripts that may serve functions in the interaction between the germline and cells of the somatic gonad. These results expand our knowledge of the transcriptomic space of the C. elegans germline and lay a foundation on which future studies of the germline can be based upon.

Validated Negative Regions (VNRs) in the VISTA Database might be Truncated Forms of Bona Fide Enhancers.

The VISTA enhancer database is a valuable resource for evaluating predicted enhancers in humans and mice. In addition to thousands of validated positive regions (VPRs) in the human and mouse genomes, the database also contains similar numbers of validated negative regions (VNRs). It is previously shown that the VPRs are on average half as long as predicted overlapping enhancers that are highly conserved and hypothesize that the VPRs may be truncated forms of long bona fide enhancers. Here, it is shown that like the VPRs, the VNRs also are under strong evolutionary constraints and overlap predicted enhancers in the genomes. The VNRs are also on average half as long as predicted overlapping enhancers that are highly conserved. Moreover, the VNRs and the VPRs display similar cell/tissue-specific modification patterns of key epigenetic marks of active enhancers. Furthermore, the VNRs and the VPRs show similar impact score spectra of in silico mutagenesis. These highly similar properties between the VPRs and the VNRs suggest that like the VPRs, the VNRs may also be truncated forms of long bona fide enhancers.

Annotations of four high-quality indigenous chicken genomes identify more than one thousand missing genes in subtelomeric regions and micro-chromosomes with high G/C contents.

BACKGROUND: Although multiple chicken genomes have been assembled and annotated, the numbers of protein-coding genes in chicken genomes and their variation among breeds are still uncertain due to the low quality of these genome assemblies and limited resources used in their gene annotations. To fill these gaps, we recently assembled genomes of four indigenous chicken breeds with distinct traits at chromosome-level. In this study, we annotated genes in each of these assembled genomes using a combination of RNA-seq- and homology-based approaches. RESULTS: We identified varying numbers (17,497-17,718) of protein-coding genes in the four indigenous chicken genomes, while recovering 51 of the 274 "missing" genes in birds in general, and 36 of the 174 "missing" genes in chickens in particular. Intriguingly, based on deeply sequenced RNA-seq data collected in multiple tissues in the four breeds, we found 571 ~ 627 protein-coding genes in each genome, which were missing in the annotations of the reference chicken genomes (GRCg6a and GRCg7b/w). After removing redundancy, we ended up with a total of 1,420 newly annotated genes (NAGs). The NAGs tend to be found in subtelomeric regions of macro-chromosomes (chr1 to chr5, plus chrZ) and middle chromosomes (chr6 to chr13, plus chrW), as well as in micro-chromosomes (chr14 to chr39) and unplaced contigs, where G/C contents are high. Moreover, the NAGs have elevated quadruplexes G frequencies, while both G/C contents and quadruplexes G frequencies in their surrounding regions are also high. The NAGs showed tissue-specific expression, and we were able to verify 39 (92.9%) of 42 randomly selected ones in various tissues of the four chicken breeds using RT-qPCR experiments. Most of the NAGs were also encoded in the reference chicken genomes, thus, these genomes might harbor more genes than previously thought. CONCLUSION: The NAGs are widely distributed in wild, indigenous and commercial chickens, and they might play critical roles in chicken physiology. Counting these new genes, chicken genomes harbor more genes than originally thought.

Artificial selection footprints in indigenous and commercial chicken genomes.

BACKGROUND: Although many studies have been done to reveal artificial selection signatures in commercial and indigenous chickens, a limited number of genes have been linked to specific traits. To identify more trait-related artificial selection signatures and genes, we re-sequenced a total of 85 individuals of five indigenous chicken breeds with distinct traits from Yunnan Province, China. RESULTS: We found 30 million non-redundant single nucleotide variants and small indels (< 50 bp) in the indigenous chickens, of which 10 million were not seen in 60 broilers, 56 layers and 35 red jungle fowls (RJFs) that we compared with. The variants in each breed are enriched in non-coding regions, while those in coding regions are largely tolerant, suggesting that most variants might affect cis-regulatory sequences. Based on 27 million bi-allelic single nucleotide polymorphisms identified in the chickens, we found numerous selective sweeps and affected genes in each indigenous chicken breed and substantially larger numbers of selective sweeps and affected genes in the broilers and layers than previously reported using a rigorous statistical model. Consistent with the locations of the variants, the vast majority (~ 98.3%) of the identified selective sweeps overlap known quantitative trait loci (QTLs). Meanwhile, 74.2% known QTLs overlap our identified selective sweeps. We confirmed most of previously identified trait-related genes and identified many novel ones, some of which might be related to body size and high egg production traits. Using RT-qPCR, we validated differential expression of eight genes (GHR, GHRHR, IGF2BP1, OVALX, ELF2, MGARP, NOCT, SLC25A15) that might be related to body size and high egg production traits in relevant tissues of relevant breeds. CONCLUSION: We identify 30 million single nucleotide variants and small indels in the five indigenous chicken breeds, 10 million of which are novel. We predict substantially more selective sweeps and affected genes than previously reported in both indigenous and commercial breeds. These variants and affected genes are good candidates for further experimental investigations of genotype-phenotype relationships and practical applications in chicken breeding programs.

High quality assemblies of four indigenous chicken genomes and related functional data resources.

Many lines of evidence indicate that red jungle fowl (RJF) is the primary ancestor of domestic chickens. Although multiple versions of RJF (galgal2-galgal5 and GRCg6a) and commercial chickens (GRCg7b/w and Huxu) genomes have been assembled since 2004, no high-quality indigenous chicken genomes have been assembled, hampering the understanding of chicken domestication and evolution. To fill the gap, we sequenced the genomes of four indigenous chickens with distinct morphological traits in southwest China, using a combination of short, long and Hi-C reads. We assembled each genome (~1.0 Gb) into 42 chromosomes with chromosome N50 90.5-90.9 Mb, amongst the highest quality of chicken genome assemblies. To provide resources for gene annotation and functional analysis, we also sequenced transcriptomes of 10 tissues for each of the four chickens. Moreover, we corrected many mis-assemblies and assembled missing micro-chromosomes 29 and 34-39 for GRCg6a. Our assemblies, sequencing data and the correction of GRCg6a can be valuable resources for studying chicken domestication and evolution.

High-quality genome assembly of a C. crossoptilon and related functional and genetics data resources.

There are four species in the Crossoptilon genus inhibiting at from very low to very high altitudes across China, and they are in varying levels of danger of extinction. To better understand the genetic basis of adaptation to high altitudes and genetic changes due to bottleneck, we assembled the genome (~1.02 Gb) of a white eared pheasant (WT) (Crossoptilon crossoptilon) inhibiting at high altitudes (3,000~7,000 m) in northwest of Yunnan province, China, using a combination of Illumina short reads, PacBio long reads and Hi-C reads, with a contig N50 of 19.63 Mb and only six gaps. To further provide resources for gene annotation as well as functional and population genetics analyses, we sequenced transcriptomes of 20 major tissues of the WT individual and re-sequenced another 10 WT individuals and a blue eared pheasant (Crossoptilon auritum) individual inhabiting at intermediate altitudes (1,500~3,000 m). Our assembled WT genome, transcriptome data, and DNA sequencing data can be valuable resources for studying the biology, evolution and developing conservation strategies of these endangered species.

Underlying causes for prevalent false positives and false negatives in STARR-seq data.

Self-transcribing active regulatory region sequencing (STARR-seq) and its variants have been widely used to characterize enhancers. However, it has been reported that up to 87% of STARR-seq peaks are located in repressive chromatin and are not functional in the tested cells. While some of the STARR-seq peaks in repressive chromatin might be active in other cell/tissue types, some others might be false positives. Meanwhile, many active enhancers may not be identified by the current STARR-seq methods. Although methods have been proposed to mitigate systematic errors caused by the use of plasmid vectors, the artifacts due to the intrinsic limitations of current STARR-seq methods are still prevalent and the underlying causes are not fully understood. Based on predicted cis-regulatory modules (CRMs) and non-CRMs in the human genome as well as predicted active CRMs and non-active CRMs in a few human cell lines/tissues with STARR-seq data available, we reveal prevalent false positives and false negatives in STARR-seq peaks generated by major variants of STARR-seq methods and possible underlying causes. Our results will help design strategies to improve STARR-seq methods and interpret the results.

Antisense transcription and its roles in adaption to environmental stress in E. coli.

It has been reported that a highly varying proportion (1% ∼ 93%) of genes in various prokaryotes have antisense RNA (asRNA) transcription. However, the extent of the pervasiveness of asRNA transcription in the well-studied E. coli K12 strain has thus far been an issue of debate. Furthermore, very little is known about the expression patterns and functions of asRNAs under various conditions. To fill these gaps, we determined the transcriptomes and proteomes of E. coli K12 at multiple time points in five culture conditions using strand-specific RNA-seq, differential RNA-seq, and quantitative mass spectrometry methods. To reduce artifacts of possible transcriptional noise, we identified asRNA using stringent criteria with biological replicate verification and transcription start sites (TSSs) information included. We identified a total of 660 asRNAs, which were generally short and largely condition-dependently transcribed. We found that the proportions of the genes which had asRNA transcription highly depended on the culture conditions and time points. We classified the transcriptional activities of the genes in six transcriptional modes according to their relative levels of asRNA to mRNA. Many genes changed their transcriptional modes at different time points of the culture conditions, and such transitions can be described in a well-defined manner. Intriguingly, the protein levels and mRNA levels of genes in the sense-only/sense-dominant mode were moderately correlated, but the same was not true for genes in the balanced/antisense-dominant mode, in which asRNAs were at a comparable or higher level to mRNAs. These observations were further validated by western blot on candidate genes, where an increase in asRNA transcription diminished gene expression in one case and enhanced it in another. These results suggest that asRNAs may directly or indirectly regulate translation by forming duplexes with cognate mRNAs. Thus, asRNAs may play an important role in the bacterium's responses to environmental changes during growth and adaption to different environments. IMPORTANCE: The cis -antisense RNA (asRNA) is a type of understudied RNA molecules in prokaryotes, which is believed to be important in regulating gene expression. Our current understanding of asRNA is constrained by inconsistent reports about its identification and properties. These discrepancies are partially caused by a lack of sufficient samples, biological replicates, and culture conditions. This study aimed to overcome these disadvantages and identified 660 putative asRNAs using integrated information from strand-specific RNA-seq, differential RNA-seq, and mass spectrometry methods. In addition, we explored the relative expression between asRNAs and sense RNAs and investigated asRNA regulated transcriptional activity changes over different culture conditions and time points. Our work strongly suggests that asRNAs may play a crucial role in bacterium's responses to environmental changes during growth and adaption to different environments.

ARBic: an all-round biclustering algorithm for analyzing gene expression data.

Identifying significant biclusters of genes with specific expression patterns is an effective approach to reveal functionally correlated genes in gene expression data. However, none of existing algorithms can simultaneously identify both broader and narrower biclusters due to their failure of balancing between effectiveness and efficiency. We introduced ARBic, an algorithm which is capable of accurately identifying any significant biclusters of any shape, including broader, narrower and square, in any large scale gene expression dataset. ARBic was designed by integrating column-based and row-based strategies into a single biclustering procedure. The column-based strategy borrowed from RecBic, a recently published biclustering tool, extracts narrower biclusters, while the row-based strategy that iteratively finds the longest path in a specific directed graph, extracts broader ones. Being tested and compared to other seven salient biclustering algorithms on simulated datasets, ARBic achieves at least an average of 29% higher recovery, relevance and[Formula: see text] scores than the best existing tool. In addition, ARBic substantially outperforms all tools on real datasets and is more robust to noises, bicluster shapes and dataset types.

Transcriptomic changes in single yeast cells under various stress conditions.

BACKGROUND: The stress response of Saccharomyces cerevisiae has been extensively studied in the past decade. However, with the advent of recent technology in single-cell transcriptome profiling, there is a new opportunity to expand and further understanding of the yeast stress response with greater resolution on a system level. To understand transcriptomic changes in baker's yeast S. cerevisiae cells under stress conditions, we sequenced 117 yeast cells under three stress treatments (hypotonic condition, glucose starvation and amino acid starvation) using a full-length single-cell RNA-Seq method. RESULTS: We found that though single cells from the same treatment showed varying degrees of uniformity, technical noise and batch effects can confound results significantly. However, upon careful selection of samples to reduce technical artifacts and account for batch-effects, we were able to capture distinct transcriptomic signatures for different stress conditions as well as putative regulatory relationships between transcription factors and target genes. CONCLUSION: Our results show that a full-length single-cell based transcriptomic analysis of the yeast may help paint a clearer picture of how the model organism responds to stress than do bulk cell population-based methods.

Accurate prediction of functional states of cis-regulatory modules reveals common epigenetic rules in humans and mice.

BACKGROUND: Predicting cis-regulatory modules (CRMs) in a genome and their functional states in various cell/tissue types of the organism are two related challenging computational tasks. Most current methods attempt to simultaneously achieve both using data of multiple epigenetic marks in a cell/tissue type. Though conceptually attractive, they suffer high false discovery rates and limited applications. To fill the gaps, we proposed a two-step strategy to first predict a map of CRMs in the genome, and then predict functional states of all the CRMs in various cell/tissue types of the organism. We have recently developed an algorithm for the first step that was able to more accurately and completely predict CRMs in a genome than existing methods by integrating numerous transcription factor ChIP-seq datasets in the organism. Here, we presented machine-learning methods for the second step. RESULTS: We showed that functional states in a cell/tissue type of all the CRMs in the genome could be accurately predicted using data of only 1~4 epigenetic marks by a variety of machine-learning classifiers. Our predictions are substantially more accurate than the best achieved so far. Interestingly, a model trained on a cell/tissue type in humans can accurately predict functional states of CRMs in different cell/tissue types of humans as well as of mice, and vice versa. Therefore, epigenetic code that defines functional states of CRMs in various cell/tissue types is universal at least in humans and mice. Moreover, we found that from tens to hundreds of thousands of CRMs were active in a human and mouse cell/tissue type, and up to 99.98% of them were reutilized in different cell/tissue types, while as small as 0.02% of them were unique to a cell/tissue type that might define the cell/tissue type. CONCLUSIONS: Our two-step approach can accurately predict functional states in any cell/tissue type of all the CRMs in the genome using data of only 1~4 epigenetic marks. Our approach is also more cost-effective than existing methods that typically use data of more epigenetic marks. Our results suggest common epigenetic rules for defining functional states of CRMs in various cell/tissue types in humans and mice.

A map of cis-regulatory modules and constituent transcription factor binding sites in 80% of the mouse genome.

BACKGROUND: Mouse is probably the most important model organism to study mammal biology and human diseases. A better understanding of the mouse genome will help understand the human genome, biology and diseases. However, despite the recent progress, the characterization of the regulatory sequences in the mouse genome is still far from complete, limiting its use to understand the regulatory sequences in the human genome. RESULTS: Here, by integrating binding peaks in ~ 9,000 transcription factor (TF) ChIP-seq datasets that cover 79.9% of the mouse mappable genome using an efficient pipeline, we were able to partition these binding peak-covered genome regions into a cis-regulatory module (CRM) candidate (CRMC) set and a non-CRMC set. The CRMCs contain 912,197 putative CRMs and 38,554,729 TF binding sites (TFBSs) islands, covering 55.5% and 24.4% of the mappable genome, respectively. The CRMCs tend to be under strong evolutionary constraints, indicating that they are likely cis-regulatory; while the non-CRMCs are largely selectively neutral, indicating that they are unlikely cis-regulatory. Based on evolutionary profiles of the genome positions, we further estimated that 63.8% and 27.4% of the mouse genome might code for CRMs and TFBSs, respectively. CONCLUSIONS: Validation using experimental data suggests that at least most of the CRMCs are authentic. Thus, this unprecedentedly comprehensive map of CRMs and TFBSs can be a good resource to guide experimental studies of regulatory genomes in mice and humans.

The intestinal microbiome associated with lipid metabolism and obesity in humans and animals.

Intestinal microbiota is considered to play an integral role in maintaining health of host by modulating several physiological functions including nutrition, metabolism and immunity. Accumulated data from human and animal studies indicate that intestinal microbes can affect lipid metabolism in host through various direct and indirect biological mechanisms. These mechanisms include the production of various signalling molecules by the intestinal microbiome, which exert a strong effect on lipid metabolism, bile secretion in the liver, reverse transport of cholesterol and energy expenditure and insulin sensitivity in peripheral tissues. This review discusses the findings of recent studies suggesting an emerging role of intestinal microbiota and its metabolites in regulating lipid metabolism and the association of intestinal microbiota with obesity. Additionally, we discuss the controversies and challenges in this research area. However, intestinal micro-organisms are also affected by some external factors, which in turn influence the regulation of microbial lipid metabolism. Therefore, we also discuss the effects of probiotics, prebiotics, diet structure, exercise and other factors on intestinal microbiological changes and lipid metabolism regulation.

Dynamic changes in fecal bacterial microbiota of dairy cattle across the production line.

BACKGROUND: Microbiota play important roles in the gastrointestinal tract (GIT) of dairy cattle as the communities are responsible for host health, growth, and production performance. However, a systematic characterization and comparison of microbial communities in the GIT of cattle housed in different management units on a modern dairy farm are still lacking. We used 16S rRNA gene sequencing to evaluate the fecal bacterial communities of 90 dairy cattle housed in 12 distinctly defined management units on a modern dairy farm. RESULTS: We found that cattle from management units 5, 6, 8, and 9 had similar bacterial communities while the other units showed varying levels of differences. Hutch calves had a dramatically different bacterial community than adult cattle, with at least 10 genera exclusively detected in their samples but not in non-neonatal cattle. Moreover, we compared fecal bacteria of cattle from every pair of the management units and detailed the number and relative abundance of the significantly differential genera. Lastly, we identified 181 pairs of strongly correlated taxa in the community, showing possible synergistic or antagonistic relationships. CONCLUSIONS: This study assesses the fecal microbiota of cattle from 12 distinctly defined management units along the production line on a California dairy farm. The results highlight the similarities and differences of fecal microbiota between cattle from each pair of the management units. Especially, the data indicate that the newborn calves host very different gut bacterial communities than non-neonatal cattle, while non-neonatal cattle adopt one of the two distinct types of gut bacterial communities with subtle differences among the management units. The gut microbial communities of dairy cattle change dramatically in bacterial abundances at different taxonomic levels along the production line. The findings provide a reference for research and practice in modern dairy farm management.

PCRMS: a database of predicted cis-regulatory modules and constituent transcription factor binding sites in genomes.

More accurate and more complete predictions of cis-regulatory modules (CRMs) and constituent transcription factor (TF) binding sites (TFBSs) in genomes can facilitate characterizing functions of regulatory sequences. Here, we developed a database predicted cis-regulatory modules (PCRMS) (https://cci-bioinfo.uncc.edu) that stores highly accurate and unprecedentedly complete maps of predicted CRMs and TFBSs in the human and mouse genomes. The web interface allows the user to browse CRMs and TFBSs in an organism, find the closest CRMs to a gene, search CRMs around a gene and find all TFBSs of a TF. PCRMS can be a useful resource for the research community to characterize regulatory genomes. Database URL: https://cci-bioinfo.uncc.edu/.

Integrative analysis of transcriptomics and metabolomics to reveal the melanogenesis pathway of muscle and related meat characters in Wuliangshan black-boned chickens.

BACKGROUND: Melanin is an important antioxidant in food and has been used in medicine and cosmetology. Chicken meat with high melanin content from black-boned chickens have been considered a high nutritious food with potential medicinal properties. The molecular mechanism of melanogenesis of skeletal muscle in black-boned chickens remain poorly understood. This study investigated the biological gene-metabolite associations regulating the muscle melanogenesis pathways in Wuliangshan black-boned chickens with two normal boned chicken breeds as control. RESULTS: We identified 25 differentially expressed genes and 11 transcription factors in the melanogenesis pathways. High levels of the meat flavor compounds inosine monophosphate, hypoxanthine, lysophospholipid, hydroxyoctadecadienoic acid, and nicotinamide mononucleotide were found in Wuliangshan black-boned chickens. CONCLUSION: Integrative analysis of transcriptomics and metabolomics revealed the dual physiological functions of the PDZK1 gene, involved in pigmentation and/or melanogenesis and regulating the phospholipid signaling processes in muscle of black boned chickens.

Antibiotic Resistance Genes and Associated Phenotypes in Escherichia coli and Enterococcus from Cattle at Different Production Stages on a Dairy Farm in Central California.

The objectives of this study were to characterize overall genomic antibiotic resistance profiles of fecal Escherichia coli and Enterococcus spp. from dairy cattle at different production stages using whole-genome sequencing and to determine the association between antimicrobial resistance (AMR) phenotypes and their corresponding genotypes. The Comprehensive Antibiotic Resistance Database (CARD) and ResFinder, two publicly available databases of antimicrobial resistance genes, were used to annotate isolates. Based on the ResFinder database, 27.5% and 20.0% of tested E. coli isolates (n = 40) harbored single and ≥3 antimicrobial resistance genes, respectively; for Enterococcus spp., we observed 87.8% and 8.2%, respectively. The highest prevalence of AMR genes in E. coli was for resistance to tetracycline (27.5%), followed by sulphonamide (22.5%) and aminoglycoside (20.0%); the predominant antimicrobial resistance genes in Enterococcus spp. targeted macrolide drugs (77.6%). Based on the CARD database, resistance to ≥3 antimicrobial classes was observed in all E. coli and 77.6% in Enterococcus spp. isolates. A high degree of agreement existed between the resistance phenotype and the presence of resistance genes for various antimicrobial classes for E. coli but much less so for isolates of Enterococcus. Consistent with prior work, fecal E. coli and Enterococcus spp. isolates from calves harbored a wide spectrum of resistance genes, compared to those from cattle at other production stages, based on the cross-sectional samples from the studied farm.

Accurate prediction of cis-regulatory modules reveals a prevalent regulatory genome of humans.

cis-regulatory modules(CRMs) formed by clusters of transcription factor (TF) binding sites (TFBSs) are as important as coding sequences in specifying phenotypes of humans. It is essential to categorize all CRMs and constituent TFBSs in the genome. In contrast to most existing methods that predict CRMs in specific cell types using epigenetic marks, we predict a largely cell type agonistic but more comprehensive map of CRMs and constituent TFBSs in the gnome by integrating all available TF ChIP-seq datasets. Our method is able to partition 77.47% of genome regions covered by available 6092 datasets into a CRM candidate (CRMC) set (56.84%) and a non-CRMC set (43.16%). Intriguingly, the predicted CRMCs are under strong evolutionary constraints, while the non-CRMCs are largely selectively neutral, strongly suggesting that the CRMCs are likely cis-regulatory, while the non-CRMCs are not. Our predicted CRMs are under stronger evolutionary constraints than three state-of-the-art predictions (GeneHancer, EnhancerAtlas and ENCODE phase 3) and substantially outperform them for recalling VISTA enhancers and non-coding ClinVar variants. We estimated that the human genome might encode about 1.47M CRMs and 68M TFBSs, comprising about 55% and 22% of the genome, respectively; for both of which, we predicted 80%. Therefore, the cis-regulatory genome appears to be more prevalent than originally thought.