Livestock genome annotation: transcriptome and chromatin structure profiling in cattle, goat, chicken and pig

Background Functional annotation of livestock genomes is a critical step to decipher the genotype-to-phenotype relationship underlying complex traits. As part of the Functional Annotation of Animal Genomes (FAANG) action, the FR-AgENCODE project aims at profiling the landscape of transcription (RNA-seq) and chromatin accessibility and conformation (ATAC-seq and Hi-C) in four livestock species representing ruminants (cattle, goat), monogastrics (pig) and birds (chicken), using three target samples related to metabolism (liver) and immunity (CD4+ and CD8+ T cells). Results Standardized protocols were applied to produce transcriptome and chromatin datasets for the four species. RNA-seq assays considerably extended the available catalog of protein-coding and non-coding transcripts. Gene expression profiles were consistent with known metabolic/immune functions and revealed differentially expressed transcripts with unknown function, including new lncRNAs in syntenic regions. The majority of ATAC-seq peaks of chromatin accessibility mapped to putative regulatory regions, with an enrichment of predicted transcription factor binding sites in differentially accessible peaks. Hi-C provided the first set of genome-wide maps of three-dimensional interactions across livestock and showed consistency with results from gene expression and chromatin accessibility in topological compartments of the genomes. Conclusions We report the first multi-species and multi-assay genome annotation results obtained by a FAANG pilot project. The global consistency between gene expression and chromatin structure data in these four livestock species confirms previous findings in model animals. Overall, these results emphasize the value of FAANG for research on domesticated animals and strengthen the importance of future meta-analyses of the reference datasets being generated by this community on different species.


Background
The majority of complex trait-associated loci lie outside protein-coding regions, and comparative genomics studies have shown that the majority of mammalianconserved and recently adapted regions consisted of non-coding elements [1][2][3]. This evidence prompted the first large-scale efforts into genome annotation for human and model organisms [4][5][6]. The genome-wide annotation maps generated by these projects helped to shed light on the main features of genome activity. For example, chromatin conformation or transcription factor occupancy at regulatory elements can often be directly tied to the biology of the specific cell or tissue under study [3,7,8]. Moreover, although a subset of core regulatory systems are largely conserved across humans and mice, the underlying regulatory systems often diverge substantially [9][10][11], implying that understanding the phenotypes of interest requires organism-specific information for any specific physiological phase, tissue and cell.
The Functional Annotation of Animal Genomes (FAANG) initiative (www.faang. org) aims to support and coordinate the community in the endeavor of creating reference functional maps of the genomes of domesticated animals across different species, tissues and developmental stages, with an initial focus on farm and companion animals [12,13]. To achieve high quality functional annotation of these genomes, FAANG carries out activities to standardize core biochemical assay protocols, coordinate and facilitate data sharing. It also establishes an infrastructure for the analysis of genotype-to-phenotype data (www.faang.org, www.faang-europe.org). Substantial efforts are being dedicated to farm animal species, as deciphering the genotype-to-phenotype relationships underlying complex traits such as production efficiency and disease resistance is a prerequisite for exploiting the full potential of livestock [12].
Here we report the main results of a pilot project (FR-AgENCODE) launched at the beginning of the FAANG initiative. The broad aim was to generate standardized FAANG reference datasets from four livestock species (cattle, goat, chicken and pig) through the adaptation and optimization of molecular assays and analysis pipelines. We first collected a panel of samples from more than 40 tissues from two males and two females of four species: Bos taurus (cattle, Holstein breed), Capra hircus (goat, Alpine breed), Gallus gallus (chicken, White Leghorn breed) and Sus scrofa (pig, Large White breed), generating a total of 4,115 corresponding entries registered at the EMBL-EBI BioSamples database (see Methods). For molecular characterization, three tissues were chosen to represent a "hub" organ (liver) and two broad immune cell populations (CD4+ and CD8+ T cells). This allowed obtaining a partial representation of energy metabolisms and immunity functions, as well as the optimization of the protocols for experimental assays for both tissue-dissociated and primary sorted cells. In addition to the transcriptome, we analyzed chromatin accessibility by the assay for transposase-accessible chromatin using sequencing (ATACseq, [14]), and we characterized the three-dimensional (3D) genome architecture by coupling proximity-based ligation with massively parallel sequencing (Hi-C, [15]) ( Figure 1). Using this combination of tissues/assays, we assessed the expression of a large set of coding and non-coding transcripts in the four species, evaluated their patterns of differential expression in light of chromatin accessibility at promoter regions, and characterized active and inactive topological domains of these genomes.
The integrative analysis showed a global consistency of all data, emphasizing the value of a coordinated action to improve the genomic annotation of livestock species.

Results and discussion
High-depth RNA-seq assays provide gene expression profiles in liver and immune cells from cattle, goat, chicken and pig For each animal (two males, two females) of the four species, we used RNA-seq to profile the transcriptome of liver, CD4+ and CD8+ cells (see Methods, Figure 1 and Table S1). Briefly, stranded libraries were prepared from polyA+ selected RNAs longer than 200bp, and sequencing was performed on an Illumina HiSeq3000 (see Methods). Between 250M (chicken) and 515M (goat) read pairs were obtained per tissue, of which 94% (chicken) to 98% (pig) mapped on their respective genome using STAR [16,17] (see Methods and Tables S2-S4). As an initial quality control step, we processed the mapped reads with RSEM [18] to estimate the expression levels of all genes and transcripts from the Ensembl reference annotation (hereafter called "reference" genes/transcripts/annotation) ( Figure S1, Supplementary data file 1). As expected, a large majority of the reads (from 62% in cattle to 72% in goat) fell in annotated exons of the reference genes ( Figure S2). In spite of the specialized scope of this initial study limited to liver and immune cells, a large share of all reference genes were detected (from 58% in chicken to 65% in goat) with an expression level higher than 0.1 TPM for the same transcript isoform in at least two samples (see Methods and Table 1).
For each species, we explored the similarity among samples using the expression profiles of the reference genes. Principal component analysis (PCA) revealed quite consistent patterns across species, where the first principal component (explaining 84 to 91% of the variability among samples) clearly separated samples according to their tissue of origin (liver vs. T cells). A more moderate yet systematic separation was observed between CD4+ and CD8+ T cells on the second principal component ( Figure S3). The consistency of these patterns across species supports the reliability of our RNA-seq data.
As expected, the gene expression ratio in males versus females was globally uniform genome-wide, with the exception of the sex chromosomes ( Figure 2). Dosage compensation (leading for instance to X chromosome inactivation in mammals) can indeed be observed in all species but chicken, which is consistent with previous reports on dosage compensation in mammals but not in avian species [19] (Figure 2).
To gain insights into the pattern of gene expression across species, we hierarchically clustered our samples using the expression of 9,461 genes found to be orthologous across the four species (see Methods and Supplementary data file 2). We observed that, regardless of the species, all liver samples clustered together while immune cell samples clustered first by species and then by type (i.e. CD4+ separately from CD8+) ( Figure 3). For most species and tissues, samples also clustered by sex ( Figure 3). This observation confirms the higher conservation of liver compared to blood gene expression across vertebrates, and highlights the fact that when RNAseq is performed on multiple tissues and species, samples could primarily cluster by either [20].
Most of the known reference genes are differentially expressed between liver and T cells To provide functional evidence supporting our RNA-seq data, we performed a differential gene expression analysis across tissues per species. The number of exonic reads was first counted for each gene of the reference annotation and was then TMM-normalized ([21], see Methods). Taking into account the specificities of our experimental design, in which samples from different tissues come from the same animal, we fitted generalized linear models (GLM) to identify genes with differential expression in: 1 All pairwise comparisons between liver, CD4+ and CD8+ 2 Liver vs. T cells globally 3 Males vs. females in liver and in T cells. As expected, liver and T cells were characterized by the largest differences in their respective RNA composition. In all species, the comparison between both tissue types yielded the majority of differentially expressed genes identified (from 7,000 genes in chicken to 10,500 genes in goat), while fewer differentially expressed genes were detected between CD4+ and CD8+ samples (Table S5 and Supplementary  data file 3). In addition, Gene Ontology (GO) analysis provided results in line with the role of liver in metabolism and of T cells in immunity for all species ( Figures  S4-7 and Methods).
Due in part to the limited number of male and female replicates, fewer genes could be revealed by the sex-specific differential analysis (Table S5). Among them, we identified 12 orthologous genes that showed consistent overexpression either in males or females across all the mammalian species, with various degrees of agreement considering the literature. We could indeed find previous reports of sex-specific effects or functions for 6 of these 12 orthologous genes (PLA2G7, L1CAM, PG-LYRP2, THBS2, MYBL1 and IGF2BP3 ) but not for five others (PATJ, CRYAB, NREP, PSAT1 and MCM2 ), [22][23][24][25][26][27]. For instance, the MYBL1 transcription factor, which is consistently overexpressed in female T-cells in our data, is described in mammals as a master regulator of meiosis in males [28]. Interestingly, the last of the 12 genes, CUX2, which was overexpressed in liver vs. T cells, was consistently overexpressed in males, as seen in human but not in mouse, where overexpression was observed in females [29,30].
As an additional control, we performed hierarchical clustering of RNA-seq samples by focusing on a manually curated list of T cell-or liver-related transcription factors (TF) (Supplementary data file 4). As expected, this also led to a clear separation between liver and T cells in all cases ( Figures S8-S9).
Few immune genes differentiate CD4+ and CD8+ cell populations In line with the results of the hierarchical clustering (Figure 3), most orthologous genes found to be differentially expressed in CD4+ versus CD8+ within species showed high variability of expression levels across species (not shown). This variability is likely caused by the natural heterogeneity in the relative proportions of cell subtypes among the different species, as already reported between mammals [31, 32]. Nevertheless, 39 orthologous genes could consistently differentiate CD4+ and CD8+ in the four species, including both mammals and chicken. Among those, 10 and 29 genes showed significant overexpression in CD4+ and in CD8+ cells respectively (Table S6). We searched for these genes in the baseline expression dataset of human CD4+ and CD8+ αβ T cell subsets generated by the Blueprint Epigenome Project [33, 34] and considered their relative enrichment in either cell subset (Table  S6). With one exception (ACVRL1 ), all genes were found to be expressed in human CD4+ and/or CD8+ αβ T cells and 25 of them show relative enrichment in CD4+ (or CD8+) human cells, which is consistent with our data across the four species. Out of these 25 genes, six and eight genes could be respectively associated with CD4+ and CD8+ T cell differentiation, activation and function (Table S6).
The limited number of genes composing the common signature of CD4+ and CD8+ T cell functions is primarily caused by the limited resolution of the sorting strategies (see Methods), leading to heterogeneity in the relative proportions of cell subtypes among the different species. As an example, the sorting of porcine CD8α+ cells included TCR-γδ T cells [31,35], an important subpopulation which is known to be present at higher frequency in pigs and ruminants (15-60%) compared to mice and humans (<5%) [32]. Despite these limits, and considering that only partial knowledge is available across the farm species used in this study, our results showed consistency with a previous study comparing human and mouse; although global transcription profiles were found conserved between corresponding immune cell lineages, several hundred genes show divergent expression [36].

Analysis of new transcripts improves and extends gene structure annotation
In order to test if our data could improve the reference gene annotation for each species, we used STAR and Cufflinks to identify all transcripts present in our samples and predict their exon-intron structures. We then quantified their expression in each sample using STAR and RSEM (Methods, Figure S1) and only retained the transcripts and corresponding genes expressed in at least two samples with TPM ≥ 0.1. We identified between 58,000 and 85,000 transcripts depending on the species (Table 1, Supplementary data file 5), hereafter called "FR-AgENCODE transcripts".
To characterize these FR-AgENCODE transcripts with respect to the reference ones, we grouped them into four positional classes (see Methods): (1) known: a FR-AgENCODE transcript whose exon-intron structure is strictly identical to a reference transcript (i.e. exact same exons and introns); (2) extension: same as (1), but the FR-AgENCODE transcript extends a reference transcript by at least 1 bp on at least one side; (3) alternative: a FR-AgENCODE transcript that shares at least one intron with a reference transcript but does not belong to the previous categories (only spliced transcripts can be in this class); and (4) novel : a FR-AgENCODE transcript that does not belong to any of the above categories. We found that the majority of the FR-AgENCODE transcripts were of the alternative class, therefore enriching the reference annotation with new splice variants of known genes. We note that the proportion of completely novel transcripts is relatively high for cattle, which may be due to the lack of RNA-seq data in the Ensembl gene annotation pipeline of the UMD3.1 version, leading to a less complete reference annotation (Figure S10A-C, Table 1 and Table S7).
In order to further validate the FR-AgENCODE gene annotation and identify interesting new transcripts involved in immunity and metabolism, we performed a gene expression analysis similar to the one performed on reference genes (see above and Methods). PCAs were highly similar to those previously found, with a clear separation between liver and T cell samples on the first component and a weaker separation on the second component between the two types of T cells ( Figure S11). Using the same approach as that used for the reference genes, we performed a differential expression (DE) analysis on the FR-AgENCODE genes (see above and Methods). As expected, the number of DE genes between T cells and liver was again higher than between CD4+ and CD8+ cells (Table S8 and Supplementary data file 6). The proportion of DE genes was smaller than for the reference genes (Table S5), which is expected given their globally lower gene expression levels. We also found a high level of consistency between the sets of DE genes found for the reference versus the FR-AgENCODE genes (more than 88% of the DE reference genes correspond to a DE FR-AgENCODE gene). Among the DE genes of the FR-AgENCODE set, between 202 (chicken) and 1,032 (goat) of them have at least one transcript that is classified as coding mRNA and does not overlap any reference gene on the same strand. This highlights the potential to identify novel interesting candidates for further functional characterization.

Identification, classification and comparative analysis of lncRNAs
Using the FEELnc method [37] trained with mRNAs from the reference annotation as the coding set and shuffled mRNA sequences as the non-coding training set (Methods), we identified 22,724, 13,864, 7,502 and 12,587 putative lncRNA transcripts for cattle, goat, chicken and pig, respectively (Tables S9-10 and Supplementary data file 7), many more than compared to the reference (Ensembl and NCBI) annotations (0, 4,483, 4,643 and 875 lncRNAs respectively). The smaller number of lncRNAs found in chicken was due to the larger number of non-assembled small scaffolds in that genome, and can be extended to include the 2,718 additional unclassified lncRNAs located on them (217, 707 and 83 for the other species). Consistent with previous reports in several species including human [38], dog [37], and chicken [39], the predicted lncRNA genes were much less expressed than the reference protein-coding genes ( Figure S12). The structural features of these predicted lncRNA transcripts were consistent between the four species: lncRNAs are spliced but with fewer exons (1.5 vs. 10) and higher median exon length (660 vs. 130bp) compared to mRNAs ( Figure S12). lncRNAs are also smaller than mRNAs (1800 vs. 3600bp). Notably, the lower number of exons and consequent smaller size of lncRNAs compared to mRNAs could also be due to the weaker expression of lncR-NAs, which makes it more difficult for computational methods to properly identify their structure [40].
In addition to the coding/non-coding classification, FEELnc can also categorize lncRNAs as intergenic or intragenic based on their genomic positions with respect to a provided set of reference genes (usually protein coding genes), and considering their transcription orientation with respect to these reference genes. Table S10 shows this distribution for all species. Similarities were observed between species with, on average, 6% to 12% intragenic lncRNA genes and more than 88% intergenic lncRNA genes, as observed in human [38] and chicken [39], respectively.
We and others previously showed a sharp decrease in sequence conservation of lncRNAs when the phylogenetic distance increases [38,39,41], in particular between chicken and human that diverged 300M years ago. We therefore analyzed lncRNA conservation between the four livestock species using a synteny approach based on the orthology of protein-coding genes surrounding the lncRNA and not on the lncRNA sequence conservation itself ([39], see Methods). We found 73 such conserved lncRNAs across cattle, goat and pig, 19 across cattle, chicken and pig, and 6 across all four species (Supplementary data file 7). All were expressed in these species and located in the same orientation with respect to the flanking orthologous genes. An example of such a conserved lncRNA, hitherto unknown in our four species, is provided in Figure 4. In human, this lncRNA was named CREMos for "CREM opposite sense" since it is in divergent position with respect to the neighboring CREM protein-coding gene. The synteny is conserved across species from fishes to mammals. Interestingly, the CREMos lncRNA is overexpressed in T cells while the CREM protein-coding gene is overexpressed in liver in goat, cattle and chicken ( Figure 4).
Finally, we performed a comparative transcriptomic approach, using the 73 and 19 lncRNAs conserved between three species "cattle-goat-pig" and "cattle-chicken-pig" respectively. The clustering obtained with these two gene sets shows a clear tissue separation between liver and T cells but does not separate species according to their evolutionary distance, as was observed with the 9,461 orthologous protein-coding genes from the four species ( Figure S13). Such results, in line with studies using multiple tissues from 11 tetrapod species [42], probably reflects a faster evolution of lncRNA sequences with less constraints than protein coding sequences and a specific pattern of expression.

RNA-seq mapped reads can confirm transcript 3' end positions
In order to confirm the 3' end positions of the annotated transcripts, we looked for non genomic stretches of As in RNA-seq mapped reads that could originate from unmapped polyA tails (see Methods). After filtering for internal priming and other artefacts, we obtained from 150,000 (chicken) to 250,000 (pig) distinct positions of potential polyA sites (Table S11 and Supplementary data file 8). Neighboring sites in the same orientation were further merged into 58,000 (chicken) to 108,000 (pig) polyA clusters, among which 81% contain a single polyA position. Despite their weak coverage that involved less than 0.02% of the genomic space in all species, these clusters tended to fall close to transcription termination sites (TTS). Approximately 11-14% of the reference and FR-AgENCODE transcripts had at least one polyA cluster near their TTS on the same strand (Table S11). The distribution of polyA clusters around reference and novel genes showed a clear peak at the 3' end of genes in all species ( Figure 5 and Figure S14), indicating that predicted polyAs might help support TTS positions. A quantitative analysis correlating the number of polyA-supporting reads 1Kb around the TTS of reference genes and the total number of RNA-seq reads in their exons showed only a moderate correlation (around 0.5, data not shown). Notably, most of the largest polyA clusters were found in the mitochondrial chromosome, where the cluster density was several orders of magnitude higher than genome-wide (from 11 to 30/Kb vs. 0.04 to 0.06/Kb). This suggests a high abundance of internally polyadenylated truncated mRNAs that may represent prematurely terminated transcripts or stalled degradation intermediates, as previously suggested from observations in human [43,44].
Mitochondrial polyadenylation thus appears to be prevalent across a wide range of vertebrates, from mammals to birds.
ATAC-seq peaks allow the characterization of accessible regions of chromatin in the four species We used ATAC-seq to profile the accessible chromatin of liver and of CD4+ and CD8+ immune cells in animals from the four species. Between 480M (chicken) and 950M (pig) ATAC-seq fragments were sequenced per species, and were processed by a standard pipeline: reads were trimmed, mapped to the genome (requiring a MAPQ higher than 10) and filtered for mitochondrial and PCR duplicate reads (see Methods). Most reads were filtered out by the MAPQ 10 and PCR duplicate filters ( Figure S15). Peaks were then called in each tissue separately (see Methods), resulting in between 26,000 (pig, liver) and 111,000 (pig, cd8) peaks per tissue (Table S12). Those peaks were further merged into a unique set of peaks per species, resulting in between 75,000 (goat) and 149,000 (pig) peaks (Table S12 and Supplementary data file 9), and covering between 1 and 5% of the genome. The average peak size was around 600 bp for all species, except for chicken where it was less than 500 bp. Merging tissue peaks did not result in much wider peaks ( Figure S16).
Similarly to GWAS variants associated to human diseases [3], the vast majority of our ATAC-seq peaks were either intronic or intergenic (Figure S17, Methods and Supplementary data file 9). In particular, from 38% (goat) to 55% (cattle) of the peaks lie at least 5Kb away from any reference gene ( Figure S17). These peaks potentially represent many new candidate regulatory regions that could be further correlated with genetic variants known from previous studies to be associated with phenotypes. Focusing on the promoter regions of the genes, we noted that around 15% of the ATAC-seq peaks lie at 5Kb or less from a reference TSS (Transcription Start Site). More precisely, the ATAC-seq peak distribution within and around reference genes showed a clear signal at the TSS for all species, thereby supporting the quality of both the reference gene annotation, and of our ATAC-seq data and processing ( Figure 5). Repeating the distribution analysis of the ATAC-seq peaks within and around the FR-AgENCODE transcripts that were not in the reference annotation (i.e. not from the known class; Figure S14E-H) resulted in a similar enhanced 5' signal (Figure S14A-D). This supports the quality of the structural annotations we propose, and the fact that ATAC-seq peaks can be used to support TSS prediction.
The ATAC-seq peaks of each species were quantified in each sample and resulting read counts were normalized using a loess correction (see Methods). Differential analyses similar to those used for RNA-seq genes were then performed on normalized counts (see Methods). The number of Differentially Accessible (DA) peaks between T cells and liver per species was between 4,800 (goat) and 13,600 (chicken), and, with the exception of pig, there were roughly as many more accessible peaks in T cells than in liver (Table S13 and Supplementary data file 10).
As accessible regions of the chromatin can contain regulatory sites for DNA binding proteins, we computed the density of transcription factor binding sites (TFBS) in ATAC-seq peaks genome-wide. Interestingly, the TFBS density was significantly higher in ATAC-seq peaks with a differential accessibility (p-value < 7.10e-4 for goat and < 2.2e-16 for chicken and pig, Wilcoxon test, see Methods). This was also observed for distal ATAC-seq peaks, at least 5kb away from promoters (data not shown). This suggests that differentially accessible peaks are more likely to have a regulatory role than the other ones.
Integrative analysis of RNA-seq and ATAC-seq data suggests complex regulatory mechanisms of gene expression We then investigated the correlation between chromatin accessibility and gene expression by comparing RNA-seq and ATAC-seq data. Given the specific distribution of the ATAC-seq signal, we initially focused on the gene promoter regions and considered proximal chromatin peaks (at a maximum distance of 1Kb from the TSS) to assign a promoter accessibility value to each gene.
Within each sample, genes with highly accessible promoters showed higher expression values globally ( Figure S18), as already reported in mouse and human for instance [45]. In pig and goat, the number of available samples allowed us to compute for each gene the correlation between promoter accessibility and gene expression across all samples ( Figure S19 and Methods). Interestingly, while the distribution of the correlation values appeared to be unimodal for non-differentially expressed genes, a bimodal distribution was obtained for differentially expressed genes, with an accumulation of both positive and negative correlation values ( Figure 6 Hi-C reveals topological aspects of the genome in the nucleus In order to profile the structural organization of the genome in the nucleus, we performed in situ Hi-C on liver cells from the two male and the two female samples of pig, goat and chicken. The in situ Hi-C protocol was applied as previously described [46] with slight modifications (see FAANG protocols online and Methods). Reads were processed using a bioinformatics pipeline based on HiC-Pro [47] (Methods). From 83 to 91% of the reads could be genomically mapped depending on the sample, and after filtering out all inconsistent mapping configurations we obtained a total of 185, 266 and 260M valid read pairs in goat, chicken and pig respectively (Table S14 and Figure S21). These sequencing depths allowed us to build interaction matrices (or Hi-C contact heatmaps) at the resolution of 40 and 500Kb in order to detect Topologically Associating Domains (TADs) and A/B compartments ( Figure S23).
We identified from ≈ 5,400 (chicken) to 11,000 (pig) TADs of variable sizes (150 to 220Kb on average, Table S15 and Supplementary data file 11), with a 79-86% genome-wide coverage. In order to validate the relevance of the predicted domains, obtained by maximizing the interaction densities across different resolutions [48] (see Methods), we used two approaches. First, we used the original TAD-finding metric called the Directionality Index (DI) to quantify the degree of upstream or downstream interaction bias for any genomic region [49]. As expected, the average DI along TADs showed a bias for downstream interactions at the beginning of the domains and for upstream interactions at the ends of the domains for all species ( Figure 7). Secondly, since TADs are known to be flanked by the CTCF insulating factor, at least in model organisms [49], we investigated if this could be detected by a prevalence of CTCF binding sites at the boundaries of our TADs. We observed that the density of in silico predicted CTCF sites (see Methods) consistently peaked at TAD boundaries (Figure 7), in agreement with previous reports on experimental data in human [46,49]. Altogether, these findings show a global consistency of our topological annotation with previous results and support the generality of the 3D structural organization across animals, ranging from mammals to birds.
Considering a higher-order topological organization of the chromosomes in the nucleus, we then investigated the partition of the genome into "A" and "B" compartments [15]. By applying the original method described in [15] as described in Figure S23 (see also Methods) we divided the genomes into "active" (A) and "inactive" (B) compartments and obtained from ≈ 580 to 700 compartments per genome with a mean size between 1.6Mb (chicken) and 3.4Mb (goat) (see Table S15 and Supplementary data file 11). To assess the robustness of the A/B compartment calling, we performed the same analysis on the interaction matrix of each replicate separately and compared the four resulting partitions within each species. About 80% of the informative loci (genomic bins of the matrix) were assigned to the same compartment in each animal of the three species ( Figure S23), emphasizing the stability of the detected signal in liver cells.
Considering functional aspects of genome topology, it has been proposed that A compartments represent genomic regions enriched for open chromatin and transcription compared to B compartments [49]. A higher gene expression and chromatin accessibility should therefore be expected in A versus B compartments. In all species for which we have Hi-C data, such a pattern was indeed observed, as both the average RNA-seq gene expression values and ATAC-seq chromatin accessibility values in liver samples were significantly higher in A than in B compartments ( Figure 8, p-value < 2.2e-16 for each comparison, Wilcoxon test). This result highlights the high consistency and coherence among the three types of data we have generated and validates the general relevance of the annotation we propose.

Conclusion
We report the first multi-species and multi-assay genome annotation results obtained by a pilot FAANG project. The main outcomes were consistent with our expectations: • Despite that only three tissue samples were used, a majority of the reference transcripts could be detected. Moreover, the new identified "FR-AgENCODE transcripts" considerably enrich the reference annotations of these species: the transcript repertoire is increased by about 100% in all species but cattle, where it is more than tripled. • Differential analyses of gene expression in liver and T cells yielded results consistent with known metabolism and immunity functions. • As expected, ATAC-seq data allowed to map an abundance of potential regulatory regions, and, upon integration with RNA-seq data, suggest complex regulatory mechanisms of gene expression.
• Hi-C experiments provided the first set of genome-wide 3D interaction maps of the same tissue from three livestock species. Beyond the chromosome topology annotation, the analysis shows high consistency with the gene expression and chromatin accessibility results.
As such, the FR-AgENCODE group has delivered a strong proof of concept of a successful collaborative approach at a national scale to apply FAANG guidelines to various experimental procedures and animal models. This notably includes the set up of a combination of sequencing assays on primary cells and tissue-dissociated cells, as well as a large collection of documented tissue samples available for further characterization. It also confirmed, in line with several studies in model species [4][5][6]8] the value of combining molecular assays on the same samples to simultaneously identify the transcriptomes and investigate the underlying regulatory mechanisms.
In the context of the global domesticated animal genome annotation effort, lessons learned from this pilot project confirmed conclusions drawn by the FAANG community regarding the challenges to be addressed in the future [12]. As all possible combinations of experimental assays, tissues, developmental stages and phenotypes cannot be undertaken by a single project on all species, standardized procedures for both data production and analysis need to be widely adopted [50]. Furthermore, the mosaic nature of a global annotation effort that gathers contributions from various partners worldwide emphasizes the challenge of translating recent advances from the field of data science into efficient methods for the integrative analysis of 'omics data and the importance of future metanalysis of several datasets.
Altogether, these annotation results will be useful for future studies aiming to determine which subsets of putative regulatory elements are conserved, or diverge, across animal genomes representing different phylogenetic taxa. This will be beneficial for devising efficient annotation strategies for the genomes of emerging domesticated species.

Animals and breeds
Well-characterized breeds were chosen in order to obtain well-documented samples. Holstein is the most widely used breed for dairy cattle. For goat, the Alpine breed is one of the two most commonly used dairy breeds, and for pigs, the Large white breed is widely used as a dam line. For chickens, the White Leghorn breed was chosen as it provides the genetic basis for numerous experimental lines and is widely used for egg production.
Four animals were sampled for each species, two males and two females. They all had a known pedigree. Animals were sampled at an adult stage, so that they were sexually mature and had performance records, obtained in known environmental conditions. Females were either lactating or laying eggs.
All animals were fasted at least 12 hours before slaughter. No chemicals were injected before slaughtering, animals were stunned and bled in a licensed slaughter facility at the INRA research center in Nouzilly.

Samples
Liver samples of 0.5 cm 3 were taken from the edge of the organ, avoiding proximity with the gallbladder and avoiding blood vessels. Time from slaughter to sampling varied from 5 minutes for chickens to 30 minutes for goats and pigs and 45 minutes for cattle. For the purpose of RNA-seq, samples were immediately snapfrozen in liquid nitrogen, stored in 2ml cryotubes and temporarily kept in dry ice until final storage at -80°C.
For mammals, whole blood was sampled into EDTA tubes before slaughter; at least one sampling took place well before slaughter (at least one month) and another just before slaughter, in order to obtain at least 50 ml of whole blood for separation of lymphocytes (PBMC). PBMC were re-suspended in a medium containing 10% FCS, counted, conditioned with 10% DMSO and stored in liquid nitrogen prior to the sorting of specific cell types: CD3+CD4+ ("CD4") and CD3+CD8+ ("CD8").
For chicken, spleen was sampled after exsanguination. Spleen leucocytes were purified by density-gradient separation to get rid of nucleated erythrocytes contamination and stored in liquid nitrogen prior to CD4+ and CD8+ T cell sorting.
All protocols for liver sampling, PBMC separation, splenocyte purification and T cell sorting can be found at ftp://ftp.faang.ebi.ac.uk/ftp/protocols/ samples/

Experimental assays and protocols
All assays were performed according to FAANG guidelines and recommendations, available at http://www.faang.org. All detailed protocols used for RNA extraction and libraries production for RNA-seq, ATAC-seq and Hi-C are available at http: //ftp.faang.ebi.ac.uk/ftp/protocols/assays/.

RNA extraction
Cells and tissues were homogenized in TRIzol reagent (Thermo) using an ULTRA-TURRAX (IKA-Werke) and total RNAs were extracted from the aqueous phase. They were then treated with TURBO DNase (Ambion) to remove remaining genomic DNA and then process to separate long and small RNAs using the mirVana miRNA Isolation kit. Small and long RNAs quality was assessed using an Agilent 2100 Bioanalyzer and RNA 6000 nano kits (Agilent) and quantified on a Nanodrop spectrophotometer.

RNA-seq
Stranded mRNA libraries were prepared using the TruSeq Stranded mRNA Sample Prep Kit -V2 (Illumina) on 200 ng to 1µg of total long RNA with a RNA Integrity Number (RIN) over 8 following the manufacturer's instructions. Libraries were PCR amplified for 11 cycles and libraries quality was assessed using the High Sensitivity NGS Fragment Analysis Kit DNF-474 and the Fragment Analyser system (AATI). Libraries were loaded onto a High-seq 3000 (Illumina) to reach a minimum read numbers of 100M paired reads for each library.
Hi-C In situ Hi-C libraries were performed according to [46] with a few modifications. For all species, fresh liver biopsies were dissociated using Accutase, and each resulting cell suspension was filtered using a 70 µm cell strainer. Cells were then fixed with 1% formaldehyde for 10 minutes at 37°C and fixation was stopped by adding Glycine to a final concentration of 0.125M. After two washes with PBS, cells were pelleted and kept at -80°C for long term storage. Subsequently, cells were thawed on ice and 5 million cells were processed for each Hi-C library. Cell membranes were disrupted using a potter-Elvehjem PTFE pestle and nuclei were then permeabilized using 0.5% SDS with digestion overnight with HindIII endonuclease. HindIII-cut restriction sites were then end-filled in the presence of biotin-dCTP using the Klenow large fragment and were religated overnight at 4°C. Nuclei integrity was checked using DAPI labelling and fluorescence microscopy. Nuclei were then lysed and DNA was precipitated and purified using Agencourt AMPure XP beads (Beckman Coulter) and quantified using the Qubit fluorimetric quantification system (Thermo). Hi-C efficiency was controlled by PCR using specific primers for each species and, if this step was successful, DNA was used for library production. DNA was first treated with T4 DNA polymerase to remove unligated biotinylated ends and sheared by sonication using a M220 Covaris ultra-sonicator with the DNA 550pb SnapCap microtube program (Program length: 45s; Picpower 50; DutyF 20; Cycle 200; Temperature 20°C).
Sheared DNA was then sized using magnetic beads, and biotinylated fragments were purified using M280 Streptavidin Dynabeads (Thermo) and reagents from the Nextera Mate Pair Sample preparation kit (Illumina). Purified biotinylated DNA was then processed using the TrueSeq nano DNA kit (Illumina) following the manufacturer's instructions. Libraries were amplified for 10 cycles and then purified using Agencourt AMPure XP beads. Library quality was assessed on a Fragment Analyser (AATI) and by endonuclease digestion using NheI endonuclease. Once validated, each library was sequenced on an Illumina Hi-Seq 3000 to reach a minimum number of 150M paired reads per library.
ATAC-seq ATAC-seq libraries were prepared according to Buenrostro et al. (2013) with a few modifications. For liver, cells were dissociated from the fresh tissue to get a single cell suspension. Cells were counted and 50,000 cells were processed for each assay. Transposition reaction was performed using the Tn5 Transposase and TD reaction buffer from the Nextera DNA library preparation kit (Illumina) for 30 minutes at 37°C. DNA was then purified using the Qiagen MinElute PCR purification kit. Libraries were first amplified for 5 cycles using custom-synthesized index primers (see supplementary methods) and then a second amplification was performed. The appropriate number of additional PCR cycles was determined using real-time PCR, permitting the cessation of amplification prior to saturation. The additional number of cycles needed was determined by plotting the Rn versus Cycle and then selecting the cycle number corresponding to one-third of the maximum fluorescent intensity. After PCR amplification, libraries were purified using a Qiagen MinElute PCR purification kit followed by an additional clean-up and sizing step using AMPure XP beads (160 µl of bead stock solution was added to 100 µl of DNA in EB buffer) following the manufacturer's instructions. Library quality was assessed on a BioAnalyser (Agilent) using Agilent High Sensitivity DNA kit (Agilent), and libraries were quantified using a Qubit Fluorometer (Thermo).

Bioinformatics and Data Analysis
All software used in this project along with the corresponding versions are listed in Table S3. The reference gene annotation was obtained from the Ensembl v90 release. Since Capra hircus was not part of the Ensembl release, we used the NCBI CHIR 1.0.102 annotation (see Table S2).

RNA-seq pipeline
Prior to any processing, all RNA-seq reads were trimmed using cutadapt version 1.8.3. Reads were then mapped twice using STAR v2.5.1.b [16,17]: first on the genome indexed with the reference gene annotation to quantify expression of reference transcripts, and secondly on the same genome indexed with the newly generated gene annotation (FR-AgENCODE transcripts) (see below and Figure  S1) [51]. The STAR --quantMode TranscriptomeSAM option was used in both cases in order to additionally generate a transcriptome alignment (bam) file. After read mapping and CIGAR-based softclip removal, each sample alignment file (bam file) was processed with Cufflinks 2.2.1 [52,53] with the max-intron-length (100000) and overlap-radius (5) options, guided by the reference gene annotation (--GTF-guide option) ( [51], Figure S1). All cufflinks models were then merged into a single gene annotation using Cuffmerge 2.2.1 [52,53] with the --ref-gtf option. The transcript and gene expressions on both reference and newly generated gene annotation were quantified as TPM (transcripts per million) using RSEM 1.3.0 [18] on the corresponding transcriptome alignment files ( [51], Figure S1). The newly generated transcripts were then processed with FEELnc version 0.1.0 [37] in order to classify them into "lncRNA", "mRNA" and "otherRNA" ( Figure S1, Tables S9-10, Figure S10). The newly generated transcripts with a TPM value of at least 0.1 in at least 2 samples were called FR-AgENCODE transcripts and kept as part of the new annotation. The 0.1 threshold was chosen knowing that the expression values of polyadenylated transcripts usually go from 0.01 to 10,000 [54] and that we wanted to simultaneously capture long non coding RNAs that are generally lowly expressed and reduce the risk of calling artefactual transcripts.

PCA based on gene expression
Principal Component Analysis (PCA) was performed using the mixOmics R package [55] on the RNA-seq sample quantifications of each species. This was done using the expression (TPM) of two different sets of genes: reference genes with TPM 0.1 in at least two samples ( Figure S3) and FR-AgENCODE genes with TPM 0.1 in at least two samples ( Figure S15).

Annotated gene orthologs
We used Ensembl Biomart [56] to define the set of orthologous genes across cattle, chicken and pig. Only "1 to 1" orthology relationships were kept (11,001 genes, Supplementary data file 1). Since goat was not part of the Ensembl annotation, goat gene IDs were added to this list using gene name as a correspondence term. The resulting 4-species orthologous set contained 9,461 genes (Supplementary data file 1).

RNA-seq sample hierarchical clustering
Based on the expression of the 9,461 orthologous genes in the 39 RNA-seq samples from the four species, the sample-by-sample correlation matrix was computed using the Pearson correlation of the log 10 gene TPM values. We then represented this sample by sample correlation matrix as a heatmap where the samples were also clustered using a complete linkage hierarchical clustering (Figure 3).

RNA-seq normalization and differential analysis
To perform the differential analysis of gene expression, we used the expected read counts provided by RSEM [18]. RNA-seq library size normalization factors were calculated using the weighted Trimmed Mean of M-values (TMM) approach of [57] as implemented in the R/Bioconductor package edgeR [21]. The same package was used to fit three different per-gene negative binomial (NB) generalized log-linear models [58].
• In Model 1, the expression of each gene was explained by a tissue effect; because all three tissues (liver, CD4, CD8) were collected from each animal, an animal effect was also included to account for these repeated measures: log µ gi s i = β g,tissue(i) + γ g,animal(i) , where µ gi represents the mean expression of gene g in sample i, s i the TMM normalization factor for sample i, tissue(i) ∈ {liver, CD4, CD8} and animal(i) ∈ {1, 2, 3, 4} the tissue and animal corresponding to sample i, and β g,tissue(i) and γ g,animal(i) the fixed tissue and animal effects, respectively, of gene g in sample i. Hypothesis tests were performed to identify significantly differentially expressed genes among each pair of tissues, e.g.
• Model 2 is identical to the previous model, where gene expression was modeled using both a tissue and an animal effect, with the exception that the CD4 and CD8 tissues were collapsed into a single group. In this model, the only hypothesis of interest is thus between the liver and global CD cell group: H 0g : β g,liver = β g,CD .
• In Model 3, the expression of each gene was explained using a nested factorial model, which included a sex effect as well as tissue and animal effects nested within sex: where µ gi and s i are as in Model 1, sex(i) ∈ {M, F}, tissue(i) ∈ {liver, CD4, CD8}, and animal(i) ∈ {1, 2}, respectively represent the sex, tissue, and animal nested within sex corresponding to sample i, β g,sex(i) is the sex effect for gene g in sample i, and γ g,sex(i):tissue(i) and δ g,sex(i):animal(i) are the nested interaction effects of tissue and animal with sex. Note that this model enables comparisons between sexes to be made between animals, while comparisons among tissues are performed within animals. Specifically, we focused on two sets of hypothesis tests: 1) pairwise comparisons of tissues, after correction of the sex effect; and 2) comparisons of differences between males and females within each tissue. For the former, in any of the pairwise tissue comparisons or any of the four species, 50% to 85% of reference genes (56% to 90% of the FR-AgENCODE genes) identified as DE in Model 1 or 3, were commonly identified as DE by both Models 1 and 3. All hypothesis tests were performed using likelihood-ratio tests and were corrected for multiple testing with the Benjamini-Hochberg (FDR, [59]) procedure. Genes with an FDR smaller than 5% and an absolute log-fold change larger than 1 were declared differentially expressed. Due to the small number of available samples and the fact that Model 3 was only possible for a subset of species due to incomplete designs, we primarily focused on the simpler but more robust Models 1 and 2 by default, and only used Model 3 to investigate the sex effect. For chicken, since RNAseq libraries were missing for T cells, the sex effect was estimated using liver samples only.

GO analysis of differentially expressed genes
For each species, GO term enrichment analysis was performed on the genes found to be over-expressed in liver versus T cells, and reciprocally. This analysis was done separately for each species (Figures S4-S6) and then on genes present in all species ( Figure S7), using the three following ontologies: biological process (BP), molecular function (MF) and cell compartment (CC), and using the GOstat R/Bioconductor package [60]). This analysis tests the differential genes against all the genes, and here we always restricted analyses to genes with a human ortholog.

Sample and gene hierarchical clustering based on T cell and metabolism Transcription Factors
Manually curated lists of T cell and metabolism related Transcription Factors (TF) were established (Supplementary data file 3). As a quality control we performed in each species and for both T cell and metabolism related TFs, a double hierarchical clustering of RNA-seq samples and TF genes, using one minus the squared Pearson correlation on the base 10 logarithm of the TPM as a distance and the complete linkage aggregation (Figures S8-9).

FR-AgENCODE transcript positional classification
The FR-AgENCODE transcript models were first classified according to their position with respect to reference transcripts: known the FR-AgENCODE transcript is strictly identical to a reference transcript (same intron chain and same initial and terminal exons) extension the FR-AgENCODE transcript extends a reference transcript (same intron chain but one of its two most extreme exons extends the reference transcript by at least one base pair) alternative the FR-AgENCODE transcript corresponds to a new isoform (or variant) of a reference gene, i.e. the FR-AgENCODE transcript shares at least one intron with a reference transcript but does not belong to the above categories novel the FR-AgENCODE transcript is in none of the above classes

FR-AgENCODE transcript coding classification
The FR-AgENCODE transcript models were also classified according to their coding potential. For this, the FEELnc program (release v0.1.0) was used to discriminate long non-coding RNAs from protein-coding RNAs. FEELnc includes three consecutive modules: FEELnc filter , FEELnc codpot and FEELnc classifier . The first module, FEELnc filter , filters out non-lncRNA transcripts from the assembled models, such as transcripts smaller than 200 nucleotides or those with exons strandedly overlapping exons from the reference annotation. This module was used with default parameters except -b transcript_biotype=protein_coding, pseudogene to remove novel transcripts overlapping protein-coding and pseudogene exons from the reference. The FEELnc codpot module then calculates a coding potential score (CPS) for the remaining transcripts based on several predictors (such as multi k-mer frequencies and ORF coverage) incorporated into a random forest algorithm [61]. In order to increase the robustness of the final set of novel lncRNAs and mRNAs, the options --mode=shuffle and --spethres=0.98,0.98 were set. Finally, the FEELnc classifier classifies the resulting lncRNAs according to their positions and transcriptional orientations with respect to the closest annotated reference transcripts (sense or antisense, genic or intergenic) in a 1Mb window (--maxwindow=1000000). It is worth noting that between 83 and 2,718 lncRNA transcripts were not classified because of their localization on the numerous unassembled contigs in livestock species with no annotated genes.

Syntenic conservation of lncRNAs
Briefly, a lncRNA was considered as "syntenically" conserved between two species if (1) the lncRNA was located between two orthologous protein-coding genes, (2) the lncRNA was the only one in each species between the two protein-coding genes and (3) the relative gene order and orientation of the resulting triplet was identical between species. Using these criteria, we found six triplets shared between the four species, 73 triplets shared between cattle, goat and pig, and 19 triplets shared between cattle, chicken and pig.

Annotation of 3' polyA sites
The identification of potential polyA sites was performed from the already mapped RNA-seq data in three steps: prediction, cleaning and clustering. Prediction The samToPolyA utility (Table S3) was used on RNA-seq mapped reads to find traces of sequenced yet unmapped A stretches in reads that overlap polyA sites. Taking a bam file as input, sam-ToPolyA finds stretches of As in the terminal part of the reads that had to be trimmed to allow the read to be mapped (softclip in the CIGAR string). The bam file of each sample was processed with the options --minClipped=20 --minAcontent=0.9 --discardInternallyPrimed, which outputs one bed record per potential site. For each site the number of supporting reads was computed per tissue and in total.

Cleaning
Extensive filtering was done to discard false positive sites that could have originated from internal priming during the library construction or from excessive read trimming during the mapping. The first can be detected by the presence of an A-rich region upstream of the candidate polyA site (12 As or more within the upstream 20bp), the second by a similar region downstream. Any site with such a region upstream or downstream was discarded. Unknown bases "N" were considered as "A"s for this filtering step. Clustering Consecutive positions on the same strand that were not separated by more than 10bp were then grouped into clusters of polyA sites. PolyA clusters are provided for each species as additional files in BED format, where the positions of the distinct polyA sites within the clusters and the corresponding numbers of supporting reads at each of them are indicated in fields #7 and #8 respectively. The total number of supporting reads for each cluster is reported in the score field #5.
ATAC-seq data analysis pipeline ATAC-seq reads were trimmed with trimgalore 0.4.0 using the --stringency 3, -q 20, --paired and --nextera options (Table S3). The trimmed reads were then mapped to the genome using bowtie 2 2.3.3.1 with the -X 2000 and -S options [62]. The resulting sam file was then converted to a bam file with samtools 1.3.1, and this bam file was sorted and indexed with samtools 1.3.1 [63]. The reads for which the mate was also mapped and with a MAPQ ≥ 10 were retained using samtools 1.3.1 (-F 12 and -q 10 options, [63]), and finally only the fragments where both reads had a MAPQ ≥ 10 and which were on the same chromosome were retained. Mitochondrial reads were then filtered out, as well as duplicate reads (with picard tools, MarkDuplicates subtool, see Table S14). The peaks were called using MACS2 2.1.1.20160309 [64] in each tissue separately using all the mapped reads from the given tissue (-t option) and with the --nomodel, -f BAMPE and --keep-dup all options. To get a single set of peaks per species, the tissue peaks were merged using mergeBed version 2.26.0 [65]. These peaks were also quantified in each sample by simply counting the number of mapped reads overlapping the peak.
ATAC-seq peaks were also classified with respect to the reference gene annotation using these eight genomic domains and allowing a peak to be in several genomic domains: exonic the ATAC-seq peak overlaps an annotated exon by at least one bp intronic the ATAC-seq peak is totally included in an annotated intron tss the ATAC-seq peak includes an annotated TSS tss1Kb the ATAC-seq peak overlaps an annotated TSS extended 1Kb both 5' and 3' tss5Kb the ATAC-seq peak overlaps an annotated TSS extended 5Kb both 5' and 3' tts the ATAC-seq peak includes an annotated TTS tts1Kb the ATAC-seq peak overlaps an annotated TTS extended 1Kb both 5' and 3' tts5Kb the ATAC-seq peak overlaps an annotated TTS extended 5Kb both 5' and 3' intergenic the ATAC-seq peak does not overlap any gene extended by 5KB on each side Differential analysis: normalization and model. Contrary to RNA-seq counts, ATAC-seq counts exhibited trended biases visible in log ratio-mean average (MA) plots between pairwise samples after normalization using the TMM approach, suggesting that an alternative normalization strategy was needed. In particular, trended biases are problematic as they can potentially inflate variance estimates or log foldchanges for some peaks. To address this issue, a fast loess approach [66] implemented in the normOffsets function of the R/Bioconductor package csaw [67] was used to correct differences in log-counts vs log-average counts observed between pairs of samples.
As for for RNA-seq, we used three different differential models: Model 1 for tissue pair comparisons, Model 2 for T cell versus liver comparisons and Model 3 for tissue pair and sex comparisons taking the sex effect into account (see corresponding RNAseq section for more details). In any of the pairwise tissue comparisons or any of the four species, 26% to 56% of ATAC-seq peaks identified as differentially accessible (DA) in Model 1 or 3 were commonly identified as DA by both Models 1 and 3.
Although sex-specific effects appeared to be stronger in the ATAC-seq data than in the RNA-seq data, for simplicity and consistency, we have focused our discussions on results from Models 1 and 2 in the manuscript, unless specifically discussing sex effects.
ATAC-seq peak TFBS density. In order to identify Transcription Factor Binding Sites (TFBS) genome-wide, we used the FIMO [68] software (Table S3) to look for genomic occurrences of the 519 TFs catalogued and defined in the Vertebrate 2016 JASPAR database [69] . We then intersected these occurrences with the ATAC-seq peaks of each species and computed the TFBS density in differential vs non differential ATAC-seq peaks ( Figure S18). Among the predicted TFBSs, those obtained from the CTCF motif were used to profile the resulting density with respect to Topological Associating Domains from Hi-C data (Figure 7).

Hi-C data analysis pipeline
Our Hi-C analysis pipeline includes HiC-Pro v2.9.0 [70] (Table S3) for the read cleaning, trimming, mapping (this part is internally delegated to Bowtie 2 v2.3.3.1), matrix construction, and matrix balancing ICE normalization [71]. TAD finding was done using Armatus v2.1 [48] (Table S3). Graphical visualization of the matrices was produced with the HiTC R/Bioconductor package v1.18.1 [47] (Table S3). Export to JuiceBox [72] was done through the Juicer Tools v0.7.5 (Table S3). These tools were called through a pipeline implemented in Python. Because of the high number of unassembled scaffolds (e.g. for goat) and/or micro-chromosomes (e.g. for chicken) in our reference genomes, only the longest 25 chromosomes were considered for TAD and A/B compartment calling. For these processes, each chromosome was considered separately.
The Directionality Index (DI) was computed using the original definition introduced by [49] to indicate the upstream vs. downstream interaction bias of each genomic region. Interaction matrices of each chromosome were merged across replicates and the score was computed for each bin of 40Kb. CTCF sites were predicted along the genomes by running FIMO with the JASPAR TFBS catalogue (see section "ATAC-seq peak TFBS density").
A/B compartments were obtained using the method described in [15] as illustrated in Figure S23: first, ICE-normalized counts, K ij , were corrected for a distance effect with: in which K ij is the distance-corrected count for the bins i and j, K d is the average count over all pairs of bins at distance d = d(i, j) and σ d is the standard deviation of the counts over all pairs of bins at distance d. Within-chromosome Pearson correlation matrices were then computed for all pairs of bins based on their distancecorrected counts and a PCA was performed on this matrix. The overall process was performed similarly to the method implemented in the R/Bioconductor package HiTC [47]. Boundaries between A and B compartments were identified according to the sign of the first PC (eigenvector). Since PCAs had to be performed on each chromosome separately, the average counts on the diagonal of the normalized matrix were used to identify which PC sign (+/-) should be assigned to A and B compartments for each chromosome. This allowed to automatically obtain an homogeneous assignment across chromosomes. In line with what was originally observed in human, where the first PC was the best criterion for separating A from B compartments (apart from few exceptions like the chromosome 14 for instance [15]), we also observed a good agreement between the plaid patterns of the normalized correlation matrices and the sign of the first PC ( Figure S23). To estimate the robustness of A/B compartment calling, the method was tested on each replicate separately (four animals). Since the HiTC filtering method can discard a few bins in some matrices, resulting in missing A/B labels, the proportion of bins with no conflicting labels across replicates was computed among the bins that had at least two informative replicates ( Figure S24).
Integrative analysis ATAC-seq vs. RNA-seq correlation: intra-and inter-sample analysis. For each ATAC-seq peak that overlapped a promoter region (1Kb upstream of the TSS, as suggested in Figure 5) its loess-normalized read count value (see differential analysis) was associated with the TMM-normalized expression of the corresponding gene from the reference annotation. Intra-and inter-sample correlations were then investigated: within each sample, genes were ranked according to their expressions and the distribution of the corresponding ATAC-seq values was computed for each quartile ( Figure S19). Across samples, the Pearson correlation coefficient was computed for each gene using the ATAC-seq and RNA-seq normalized values from all the available samples with both values (10 for pig, for instance), as illustrated in Figure S20.
Chromatin accessibility and gene expression in A/B compartments. To compute the general chromatin accessibility in A and B compartments, we first computed the average of the normalized read count values across all liver samples for each ATAC-seq peak. For each compartment, the mean value of all contained peaks was then reported and the resulting distributions for all A and B compartments were reported ( Figure 8).
The same approach was used to assess the general expression of genes in A and B compartments, using the average of the normalized expression values from the liver samples.
Difference between A and B distributions was tested for statistical significance using a Wilcoxon test.

Competing interests
The authors declare that they have no competing interests.  Figure 1 Experimental design overview. For each species, samples from liver and T cells of two males and two females were processed by RNA-seq, ATAC-seq and Hi-C assays. See Table S1 for a complete list of performed experiments and available datasets.    Phylogenetic representation based on the NCBI taxonomy of the 22 annotated species from fishes to mammals using iTOL [73] (A) Three gene syntenic region centered on CREMos lncRNA with NCBI IDs for lncRNAs already annotated in reference databases and distance between entities in nucleotides ; The cases where CREM and CREMos genes are overlapping are indicated by the "0*" distance. (B) Expression of the 3 genes in cattle, goat, chicken and pig: CREMos is generally less expressed in liver than in CD cells (in cattle, chicken and pig) whereas CREM is generally more expressed in liver than in CD cells (in cattle, chicken and goat). Mean coverage values of ATAC-seq peaks (y-axis, left) and polyA clusters (y-axis, right) were computed genome-wide relatively to the position of the genes (x-axis). Dotted lines indicate gene boundaries. Since genes have variable sizes, an absolute scale is used for the upstream and downstream flanking regions (from -5Kb before to +5Kb after the genes) while a relative scale is used within the gene regions (from 0 to 100% of the gene length). As expected, the ATAC-seq and the polyA signals tend to accumulate at the 5' and 3' ends of the genes respectively.  "active" compartments. B: "repressed" compartments. As Hi-C data was only available from liver, only RNA-seq and ATAC-seq values from the same samples were considered. The significant and systematic difference of gene expression and chromatin accessibility values between A and B compartments (p-values < 2.2e-16, Wilcoxon test) confirm a general consistency between RNA-seq, ATAC-seq and Hi-C data across species.