# Genetic variability of the activity of bidirectional promoters: a pilot study in bovine muscle

Cédric Meerssemann, Rabia Letaief, Véronique Léjard, Emmanuelle Rebours, Gabriel Guillocheau, Diane Esquerre, Anis Djari, Amanda Chamberlain, Christy Vander Jagt, Christophe C. Klopp, et al.

► To cite this version:

Cédric Meerssemann, Rabia Letaief, Véronique Léjard, Emmanuelle Rebours, Gabriel Guillocheau, et al.. Genetic variability of the activity of bidirectional promoters: a pilot study in bovine muscle. DNA Research, 2017, 24 (3), pp.221-233. 10.1093/dnares/dsx004. hal-01606875

HAL Id: hal-01606875

https://hal.science/hal-01606875

Submitted on 26 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Full Paper

# Genetic variability of the activity of bidirectional promoters: a pilot study in bovine muscle

Cédric Meersseman¹,², Rabia Letaief¹, Véronique Léjard¹, Emmanuelle Rebours¹, Gabriel Guillocheau¹, Diane Esquerré³, Anis Djari⁴, Amanda Chamberlain⁵,⁶, Christy Vander Jagt⁵,⁶, Christophe Klopp⁴, Mekki Boussaha¹, Gilles Renand¹, Abderrahman Maftah², Daniel Petit², and Dominique Rocha¹,*

¹GABI, INRA, AgroParisTech, Université Paris-Saclay, 78350 Jouy-en-Josas, France, ²GMA, INRA, Université de Limoges, 87060 Limoges, France, ³GenPhySE, Université de Toulouse, INRA, INPT, ENVT, 31326 Castanet Tolosan, France, ⁴SIGENAE, UR 875, INRA, 31362 Castanet-Tolosan, France, ⁵Dairy Futures Cooperative Research Centre, AgriBio, Bundoora, Victoria, Australia, and ⁶AgriBio, Department of Economic Development, Jobs, Transport & Resources, Victoria, Australia

*To whom correspondence should be addressed. Tel. +33 134 652 432. Fax. +33 134 652 478.

Email: dominique.rocha@inra.fr

Edited by Dr. Takao Sekiya

Received 2 September 2016; Editorial decision 22 January 2017; Accepted 24 January 2017

## Abstract

Bidirectional promoters are regulatory regions co-regulating the expression of two neighbouring genes organized in a head-to-head orientation. In recent years, these regulatory regions have been studied in many organisms; however, no investigation to date has been done to analyse the genetic variation of the activity of this type of promoter regions. In our study, we conducted an investigation to first identify bidirectional promoters sharing genes expressed in bovine *Longissimus thoracis* and then to find genetic variants affecting the activity of some of these bidirectional promoters. Combining bovine gene information and expression data obtained using RNA-Seq, we identified 120 putative bidirectional promoters active in bovine muscle. We experimentally validated *in vitro* 16 of these bidirectional promoters. Finally, using gene expression and whole-genome genotyping data, we explored the variability of the activity in muscle of the identified bidirectional promoters and discovered genetic variants affecting their activity. We found that the expression level of 77 genes is correlated with the activity of 12 bidirectional promoters. We also identified 57 single nucleotide polymorphisms associated with the activity of 5 bidirectional promoters. To our knowledge, our study is the first analysis in any species of the genetic variability of the activity of bidirectional promoters.

**Key words:** bidirectional promoter, cattle genome, genetic variability, muscle
# 1. Introduction

Whole-genome analyses have identified a class of regulatory regions that contain elements that initiate transcription of two different genes positioned with a head-to-head (5'-5') orientation on opposite strands of the DNA. These regions termed 'bidirectional promoters'¹⁻³ (BiPs) often have fewer than 1,000 base pairs separating the transcription start sites (TSSs) of the two genes. Bidirectional promoters have been known for years;⁴ however, recent genome-scale studies have shown that the regulation of the expression of up to 10% of the genes is controlled by bidirectional promoters.³ Bidirectional promoters are a common feature within not only the human genome but are also present in many other genomes, including yeast, plants, Invertebrates and Vertebrates.¹⁻³,⁵⁻¹⁷ It has been shown that some bidirectional promoters, such as those of the histone genes, regulate the transcription of pair of genes that need coordinated expression to maintain stoichiometric relationships,¹⁸ or regulate the coexpression of genes that function in the same biological pathway.¹⁹ Detailed studies on bidirectional promoters have also shown that some are widely conserved among Vertebrates.²⁰⁻²³ The presence of such genomic architecture in many metazoans and the high inter-species conservation seen for some bidirectional promoters suggest that they are functionally important.

Despite substantial interest in bidirectional promoters, the functional mechanisms underlying the activation of bidirectional promoters are currently not well characterized and, to our knowledge, no analysis of the variability of the activity of bidirectional promoters has been done so far.

The aim of our study was to explore bidirectional promoter activity variation and to identify if genetic variants contribute to this variability. We performed this study in bovine *Longissimus thoracis* (LT). First, we identified bovine bidirectional promoters active in muscle, combining bovine gene information and RNA-Seq data. We then experimentally validated *in vitro* a subset of these bidirectional promoters. Finally, using gene expression and whole-genome genotyping data generated for 20 samples, we explore the variability of the activity in muscle of the identified bidirectional promoters and discovered genetic variants affecting this activity. We identified 120 bovine bidirectional promoters active in LT. A total of 16 bidirectional promoters were selected and all were validated experimentally *in vitro*. As we had RNA-Seq data for 20 different samples from the same tissue, it is possible to analyse the inter-individual variability of the activity of these 120 bidirectional promoters. As previously described, expression levels of gene pairs sharing a BiP are more correlated than those of neighbouring gene pairs in a unidirectional or convergent configuration.

In addition, we found that the expression level of 77 genes is highly correlated with the activity of 12 bidirectional promoters active in muscle. We also identified 57 SNPs associated significantly with the activity of 5 bidirectional promoters.

To our knowledge, our study is the first analysis in any species of the genetic variability of the activity of bidirectional promoters.

# 2. Materials and methods

## 2.1. Animals and tissue samples

The study was conducted with 11 Limousin bull calves from a large study on the genetic determinism of beef and meat quality traits.²⁴ These eleven bull calves were not closely related to one another (for at least four generations). They were fattened in a single feedlot and fed *ad libidum* with wet corn silage. They were humanely

slaughtered in an accredited commercial slaughterhouse when they reached 16 months. LT muscle samples were dissected immediately after death and tissue samples were snap frozen in liquid nitrogen and stored at -80 °C until analysis.

## 2.2. RNA sequencing and data analysis

RNA extraction and sequencing were performed as previously described.²⁵,²⁶ Briefly, after transfer to ice-cold RNeasy RLT lysis buffer (Qiagen), LT tissue samples were homogenized using a Precellys tissue homogeniser (Bertin Technologie). Total RNA was isolated using RNeasy Midi columns (Qiagen) and then treated with RNAse-free DNase I (Qiagen) for 15 min at room temperature according to the manufacturer's protocols. The concentration of total RNA was measured with a Nanodrop ND-100 instrument (Thermo Scientific) and the quality was assessed with an RNA 6000 Nano Labchip kit using an Agilent 2100 Bioanalyzer (Agilent Technologies). All 11 samples had an RNA integrity number (RIN) value greater than eight.

The mRNA-Seq libraries were prepared using the TruSeq RNA Sample Preparation Kit (Illumina) according to the manufacturer's instructions. Briefly, Poly-A containing mRNA molecules were purified from 4 µg total RNA of each sample using oligo(dT) magnetic beads and fragmented into 150-400 bp pieces using divalent cations at 94 °C for 8 min. The cleaved mRNA fragments were converted to double-stranded cDNA using SuperScript II reverse transcriptase (Life Technologies) and primed by random primers. The resulting cDNA was purified using Agencourt AMPure® XP beads (Beckman Coulter). Then, cDNA was subjected to end-repair and phosphorylation and subsequent purification was performed using Agencourt AMPure® XP beads (Beckman Coulter). These repaired cDNA fragments were 3'-adenylated producing cDNA fragments with a single 'A' base overhang at their 3'-ends for subsequent adapter-ligation. Illumina adapters containing indexing tags were ligated to the ends of these 3'-adenylated cDNA fragments followed by two purification steps using Agencourt AMPure® XP beads (Beckman Coulter). Ten rounds of PCR amplification were performed to enrich the adapter-modified cDNA library using primers complementary to the ends of the adapters. The PCR products were purified using Agencourt AMPure® XP beads (Beckman Coulter) and size-selected (200 ± 25 bp) on a 2% agarose Invitrogen E-Gel (Thermo Scientific). Libraries were then checked on an Agilent Technologies 2100 Bioanalyzer using the Agilent High Sensitivity DNA Kit and quantified by quantitative PCR with the QPCR NGS Library Quantification kit (Agilent Technologies). After quantification, three different tagged cDNA libraries were pooled in equal ratios and a final qPCR check was performed post-pooling. Each library pool was used for 2 × 100 bp paired-end sequencing on one lane of the Illumina HiSeq2000 with a TruSeq SBS v3-HS Kit (Illumina). After sequencing, the samples were demultiplexed and the indexed adapter sequences were trimmed using the CASAVA v1.8.2 software (Illumina).

The *Bos taurus* reference transcriptome was downloaded from Ensembl (version 63, Bos_taurus.Btau_4.0.63.cdna.all.fa). To align the reads back to the assembled reference transcriptome the BWA programme (version 0.5.9-r16) was used.²⁷ Reads were mapped for each sample separately to the assembled transcriptome. The BWA default values were used for mapping. Properly paired reads with a mapping quality of at least 30 (−q = 30) were extracted from the resulting BAM file using SAMtools²⁸ for further analyses. Properly paired is defined as both left and right reads mapped in opposite
directions on the same transcript at a distance compatible with the expected mean size of the fragments (<500 bp). Custom scripts were developed to identify paired-reads mapping to single locations and with the expected distance. Read pairs mapping to separate chromosomes were discarded for the present study. The number of paired-reads uniquely aligning to transcribed regions of each transcript was calculated for all genes in the annotated transcriptome. The transcript paired-read count was calculated as the number of unique paired-reads that aligned within the exons of each transcript, based on the coordinates of mapped reads.

## 2.3. SNP identification

BWA was also used to map reads onto the bovine genome assembly UMD3.1 version.²⁹ Only reliable properly paired BWA mapped reads were considered for SNP calling. Indels were not considered because alternative splicing impedes reliable indel discovery. SNPs were called using the SAMtools software package. Genotype likelihoods were computed using the SAMtools utilities and variable positions in the aligned reads compared with the reference were called with the BCFtools utilities.³⁰ SNPs were called only for positions with a minimal mapping quality (−Q) of 30, a minimum coverage (−d) of 4 and a maximum read depth (−D) of 10,000,000.

## 2.4. Data and statistical analyses

Information on bovine genes were obtained from the Ensembl Genome Browser (version 84) using the BioMart tool.³¹ Functional annotation analysis of genes sharing putative bidirectional promoters was done using the FATIGO tool of the online software suite Babelomics.³² Genes were assigned their Ensembl identities as input for Babelomics. Only one copy of each gene was used. Default parameter settings were used for the analysis. Statistical assessment of annotation differences between the two sets of sequences (SNP-containing genes versus all the other bovine genes) was carried out for each FATIGO analysis, using the Fisher Exact Test with correction for multiple testing. Prediction of putative TFBSs was performed using the bioinformatics tools MATCH.³³ It uses a library of positional weight matrices from TRANSFAC 6.0 public version.³⁴ A cut-off selection of 0.7 and 0.75 was applied for matrix and core similarities, respectively. Spearman's rank correlation coefficients were calculated for the correlation studies using the statistical R package.

## 2.5. Bidirectional promoter validation

Each putative bidirectional promoter region was PCR amplified from bovine genomic DNA using a pair of primers. The genomic DNA sequences were retrieved from the UCSC Genome Database. Sequence repeats were masked using RepeatMasker³⁵ and Primer3³⁶ was then used to design primer pairs to amplify each putative bidirectional promoter region. PCR primers were synthesized by Eurofins MWG Operon. Restriction enzyme cutting sites for EcoRI (GAATTC) and BamHI (GGATCC) were artificially added into the PCR primers to facilitate directional cloning. Primer sequences are presented in Supplementary Table S1. Polymerase chain reactions were performed in 50 µl using 120 ng genomic DNA, 1 × PCR buffer, 1.5 mM MgCl₂, 0.2 mM of each dNTP, 0.3 µM of each primer and 1U AccuPrime GC-rich Taq DNA polymerase (Invitrogen). The following cycling protocol was used: 95°C for 15 min, followed by 35 cycles of 95°C for 1 min, 52°C for 1 min and 72°C for 1 min, and a final extension step at 72°C for 10 min. To check the quality of the amplification, 5 µl of PCR products were then analysed by gel electrophoresis with a 1% agarose gel.

The PCR products were purified using the Qiagen MinElute Gel Extraction kit (Qiagen), digested with EcoRI and BamHI (New England Biolabs), purified with Qiagen Reaction Cleanup kit, quantified and then ligated using T4 DNA ligase (New England Biolabs) to the pBiP0 vector at a 3:1 ratio. The vector was previously digested with EcoRI/BamHI and dephosphorylated with alkaline phosphatase (New England Biolabs). The litigation products were transformed into Escherichia coli DH5α competent cells (Invitrogen). Ten clones were then amplified and plasmids were purified with the Pureyield Plasmid Miniprep System DNA purification kit (Promega). Positive clones were identified by digesting plasmid DNA with EcoRI and BamHI. Digestion products were then visualized after gel electrophoresis on a 1% agarose gel. One clone carrying the right plasmid construct was then amplified and plasmids were purified with the Macherey-Nagel Midi Endotoxin-free plasmid DNA purification kit. The plasmids were then sequenced bidirectionally using BigDye Terminator v3.1 Cycle Sequencing kit (Applied Biosystems) and the primers used for the PCR amplification. After purification on Sephadex G50 superfine column (GE Healthcare), the sequencing reaction products were analyzed using 3130 Genetic Analyzer sequencer (Applied Biosystems). One clone containing for each cloning orientation the right plasmid construct was then chosen.

Murine C2C12 myoblastic cells (ATCC CRL-1772) were grown in Dulbecco's Modified Eagle's Medium with Glutamax-I (4.5 g/l glucose, Invitrogen) supplemented with 1% penicillin/streptomycin and 20% heat-inactivated fetal bovine serum. Cells were cultured at 37°C in a humidified atmosphere of 5% CO₂ in air. C2C12 (~25,000 cells/well) were seeded into 24-well plates 1 day before transfection. Transfections were performed on 40% confluent C2C12 cells with TurboFect transfection reagent (Fermentas) using 2 µl of reagent per µg of DNA plasmid, during 4 h according to the manufacturer's protocol. The cells were transfected with 0.8 µg of the vectors carrying unidirectional pCMV1/pCMV2 promoters or putative bidirectional promoters or with pBI-CMV2, pBiP-DsRed, pBiP0 control vectors.³⁷ All transient transfection experiments were done in triplicate and repeated three times. Thirty-six hours after transfection cells were washed two times with PBS 1×, and then viewed with an Axio Observer Z1 microscope (Zeiss) and images were acquired using a CoolSNAP HQ2 camera (Photometrics) driven by the Axiovision imaging system software. Analyses of micrographs were performed using the AxioVision 4.7.2 software (Zeiss). Filter sets 38HE [excitation BP470/40(HE) and emission BP525/50(HE)] and 43HE [excitation BP545/30(HE) and emission BP620/60(HE)] were used to visualize AcGFP1 and DsRed monomer signals, respectively. Each reporter gene assay was performed in triplicate and with three independent transfection experiments.

## 2.6. Whole-genome SNP genotyping

Muscle samples were obtained and genomic DNA extracted using the DNA Midi kit (Qiagen). Quality of DNA was checked using a Nanodrop ND-100 spectrophotometer (Thermo Scientific) and quantity was estimated with Quant-iT Picogreen dsDNA kit (Life Technologies) on an ABI 7900HT (Life Technologies). All DNA samples were standardized to 50 ng/µL.

DNA samples were genotyped on the BovineSNP50 Genotyping Beadchips (Illumina) at LABOGENA (Jouy-en-Josas, France) using the standard operating procedures recommended by the manufacturer.
# 3. Results and discussion

## 3.1. Bovine muscle whole-transcriptome analysis

To identify bidirectional promoters active in bovine muscle, we first obtained a global view of the bovine LT transcriptome using paired-end RNA sequencing (RNA-Seq) from twenty Limousin bull calves. We used already-published data from nine Limousin animals²⁵,²⁶ and poly(A)-enriched mRNA from 11 new Limousin bull calves were retrotranscribed and subjected to high-throughput sequencing. The 11 RNA-Seq libraries were barcode-tagged and sequenced on five lanes (3 libraries per lane) of an Illumina HiSeq2000 sequencer. Sequencing of cDNA libraries generated a total of 1,303 million raw paired-end reads with a length of 100 bases, resulting in a total of 130 gigabases. The reads were then de-multiplexed to assign reads to each sequenced sample according to its barcode index. Approximately 17–48 million paired-end reads were obtained for each library. Reads from each sample were then mapped back to the bovine reference transcriptome. We used the set of *Bos taurus* Ensembl transcripts v63 RefSeq genes as the reference transcriptome. This set contains transcripts for 22,915 known or novel genes but also pseudogenes. Based on mappings done using the Burrows-Wheeler Aligner (BWA) programme, 67–81% of the mapped reads were aligned properly paired (Supplementary Table S2). A total of 22,025 transcripts (18,055 genes) were identified, with at least one paired-end read in all samples analysed. Similar RNA-Seq read mapping rate and the number of genes identified were obtained in other RNA-Seq bovine studies. For example, Chitwood et al.³⁸ found that ~69% of the RNA-Seq reads they generated while sequencing the transcriptomes of single bovine blastocysts mapped uniquely onto the bovine genome. Raw gene expression levels were estimated by measuring the normalized count number for each transcript (number of reads per transcript divided by the total number of mapped reads, for each sample). The five most frequent transcripts are shown in Table 1. These five genes represented nearly 17% of all sequencing reads mapped to the bovine genome and are associated with muscle cell metabolism or structure. These results were consistent with the physiological role of genes expected in the surveyed tissue. Not all genes were expressed among all samples. A total of 14,700 transcripts corresponding to 12,314 different genes were detected in all 20 samples, while 840 different genes (~4%) were only expressed in one sample at a very low level (with an average less than 2 reads).

## 3.2. Identification of putative bidirectional promoters shared by genes co-expressed in bovine LT

In order to identify putative bidirectional promoters active in bovine muscle, the chromosomal locations of all the bovine genes were retrieved from the Ensembl database (version 84).³¹ The Ensembl gene start sites were used as the TSS. Pair of genes in a head-to-head arrangement and separated by less than 1 kb between their TSSs were identified and the region between paired-genes were considered

as putative bidirectional promoters. Among the 24,616 bovine genes annotated on the UMD3.1 genome assembly, we found 563 putative bidirectional promoters (Supplementary Table S3). The number of bovine bidirectional promoters we predicted is very similar to the one calculated by Xu et al.,³⁹ who found 574 putative bidirectional promoters in cattle.

To overcome confusion during the subsequent analyses, we then used four filters (Supplementary Fig. S1). First, we removed all gene pairs for which at least one gene had paralogous sequences identified in the Ensembl database and predicted 192 bidirectional promoters. Second, we retained only gene pairs associated with these putative BiPs and for which we could detect the expression for both paired-genes in our RNA-Seq dataset, leaving 171 cases. Third, to avoid quantification problems relative to transcript isoforms, genes with known alternative transcripts were excluded from the study. Finally, for validation purposes, we only selected bidirectional promoters between protein-coding gene pairs. In total, 120 BiPs remained, each with both genes expressed in muscle (Supplementary Table S4).

## 3.3. Characteristics of genes shared by bidirectional promoters active in bovine LT

To characterise functions that are significantly enriched in genes sharing bidirectional promoters active in bovine muscle, a gene ontology (GO) analysis was performed. All the 120 co-expressed gene pairs could be analysed for functional classifications and were assigned to one or more GO annotations. GO term analysis showed a significant enrichment of specific GO terms when comparing the annotations of gene pairs associated with bidirectional promoters active in muscle against all transcripts from the bovine reference transcriptome. A summary of the classification of these genes into major biological process, cell component and molecular function categories is presented in Supplementary Table S5. Genes encoding proteins found in mitochondria, ribonucleoprotein complexes or organelle inner membrane or involved in chaperone binding and cellular response to DNA damage stimulus are significantly over-represented. No significant enrichment in KEGG terms/pathways was found. These results are in agreement with previous studies.¹⁻³,⁵⁻⁷,¹²,⁴⁰

The chromosomal location of all bovine genes, of all the predicted bidirectional promoters and of the 120 BiPs active in muscle is presented in Fig. 1. Interestingly, the chromosomal distribution of these 120 bidirectional promoters does not reflect the gene content of the chromosomes (Mann-Whitney *U*-test, *P*-value <0.05). For example, we found only one (0.83%) predicted BiP active in muscle on BTA4, whereas this chromosome harbours 855 (3.48%) bovine genes. BTA11 and BTA19 have the highest number of bidirectional promoters active in muscle (both 11.67%) whereas these chromosomes have only ~4 and ~5%, respectively, of all bovine genes. The chromosomal location of the 120 BiPs active in muscle rather follows the chromosomal distribution of all bidirectional promoters

<table><caption>**Table 1.** Top five transcripts with the most assigned reads</caption><thead><tr><th>Ensembl gene ID</th><th>Ensembl transcript ID</th><th>Description</th><th>Gene symbol</th><th>% total number of reads</th></tr></thead><tbody><tr><td>ENSBTAG00000026986</td><td>ENSBTAT00000061449</td><td>Titin</td><td>TTN</td><td>4.46</td></tr><tr><td>ENSBTAG00000018204</td><td>ENSBTAT00000009327</td><td>Myosin heavy chain 1</td><td>MYH1</td><td>4.30</td></tr><tr><td>ENSBTAG00000043561</td><td>ENSBTAT00000060569</td><td>Cytochrome c oxidase subunit I</td><td>COX1</td><td>3.61</td></tr><tr><td>ENSBTAG00000007090</td><td>ENSBTAT00000012797</td><td>Myosin heavy chain 2</td><td>MYH2</td><td>2.36</td></tr><tr><td>ENSBTAG00000004965</td><td>ENSBTAT00000006534</td><td>Nucleoporin 133 kDa</td><td>NUP133</td><td>1.99</td></tr></tbody></table>
<figure><img src="image_2.png" /><figcaption>Figure 1. Chromosomal distribution of the 120 putative bovine bidirectional promoters active in muscle.</figcaption></figure>

predicted in cattle (Pearson correlation coefficient *ρ* = 0.802). For example, we found 0.89 and 0.83% predicted bidirectional promoters and BiPs active in muscle, respectively, on chromosome 4.

## 3.4. Experimental validation of some predicted bidirectional promoters

To assess the promoter activity of the predicted bidirectional promoters, some bidirectional promoters were selected and cloned into the promoterless pBiP0 reporter vector we recently constructed.³⁷ This plasmid is a dual reporter vector with the *Green Fluorescent Protein* gene from the *Aequorea coerulescens* jellyfish (*AcGFP1*) inserted in a head-to-head configuration relative to a gene expressing a monomer mutant form of the *Red Fluorescent Protein* gene from the *Discosoma striata* reef coral (*DsRed monomer*). Both fluorescent protein genes are devoid of promoter and separated by a short multiple cloning region.

The sequences of predicted bidirectional promoters were downloaded from the UCSC database.⁴¹ To ensure that all regulatory elements were included each sequence was extended until the first ATG (ATG not included) of both paired genes. Among the 120 promoter sequences 18 contain cutting sites for EcoRI or BamHI and were eliminated, since these two restriction enzymes are used for cloning into the dual reporter vector. We previously validated one of these bidirectional promoters.²³ Sixteen new bidirectional promoters were then randomly selected for validation. Information relative to these 16 predicted bidirectional promoters is indicated in Table 2. The selected promoter regions were amplified by PCR and then cloned into the pBiP0 vector. Inserts were checked by sequencing.

The promoter activities of these sequences were then assayed after transient transfections in murine myoblastic cells (C2C12) as no bovine muscle cell line was available. The pBI-CMV2 and pBiP-DsRed vectors were used as positive controls, whereas the empty pBiP0 vector was used as a negative control to verify that no reporter gene expression could be observed without promoter. The expression of the reporter genes was assessed by fluorescence microscopy 36 h after

transfection to evaluate transfection efficiencies and check reporter activities. No significant reporter gene expression was observed with the empty pBiP0 transfected cells and positive signals were therefore attributed to the inserted sequences. Fluorescence microscopy analyses showed that all the 16 selected predicted bidirectional promoters seem to drive bidirectional expression in C2C12 cells, although variability in the reporter signals was observed (Fig. 2; Supplementary Fig. S2). These good results suggest that the procedure we developed to detect bovine bidirectional promoters active in muscle is accurate and has a low false-positive rate.

## 3.5. Co-expression analysis of genes shared by bidirectional promoters active in bovine muscle

We have in our study the expression levels of genes expressed in LT for twenty samples. It is therefore possible to analyse the co-expression of each pair of genes sharing a bidirectional promoter. Using the normalized expression levels, we calculated the Spearman's rank correlation coefficient for each gene pair. We found that 34% (41/120) of the genes shared by bidirectional promoters are co-expressed (*P*-value < 0.05%) (Supplementary Table S6). In order to evaluate whether gene pairs sharing bidirectional promoters were significantly enriched for co-expressed genes, we then performed a similar co-expression analysis for pairs of non-overlapping neighbouring genes in a tail-tail orientation (convergent configuration) or in the same orientation (unidirectional configuration). We identified these two different types of gene pairs using the chromosomal locations of all the bovine genes from the Ensembl database and using the same filters as for the identification of gene pairs sharing BiPs (Supplementary Fig. S1). We found 218 and 163 gene pairs in the convergent and unidirectional configurations (Supplementary Tables S7 and S8), respectively. The Spearman's rank correlation coefficient for each gene pair was then calculated using the normalized expressed levels measured in the 20 samples. Thirty-five percent (76/218) and 30.1% (49/163) of the gene pairs in the convergent and unidirectional configurations, respectively, were co-expressed
<table><caption>**Table 2.** Information on the 16 tested predicted bovine bidirectional promoters</caption><thead><tr><th>BiP #</th><th>Ensembl gene ID</th><th>BTA</th><th>Gene start</th><th>Gene end</th><th>Strand</th><th>Description</th><th>Gene name</th><th>Distance</th></tr></thead><tbody><tr><td rowspan="2">2</td><td>ENSBTAG00000016337</td><td>23</td><td>32835596</td><td>32849025</td><td>-1</td><td>Acyl-CoA thioesterase 13</td><td>THEM2</td><td>293</td></tr><tr><td>ENSBTAG00000000365</td><td>23</td><td>32849318</td><td>32861445</td><td>1</td><td>tyrosyl-DNA phosphodiesterase 2</td><td>TDP2</td><td></td></tr><tr><td rowspan="2">12</td><td>ENSBTAG00000014646</td><td>17</td><td>74392831</td><td>74397072</td><td>-1</td><td>Mitotic spindle organizing protein 2B</td><td>MZT2</td><td>143</td></tr><tr><td>ENSBTAG00000002130</td><td>17</td><td>74397215</td><td>74412807</td><td>1</td><td>sphingomyelin phosphodiesterase 4, neutral membrane (neutral sphingomyelinase-3)</td><td>SMPD4</td><td></td></tr><tr><td rowspan="2">19</td><td>ENSBTAG00000003550</td><td>29</td><td>1032168</td><td>1056464</td><td>-1</td><td>Chromosome 11 open reading frame 54</td><td>C11orf54</td><td>163</td></tr><tr><td>ENSBTAG00000003545</td><td>29</td><td>1056627</td><td>1062580</td><td>1</td><td>TATA box binding protein (TBP)-associated factor, RNA polymerase I, D, 41kDa</td><td>TAF1D</td><td></td></tr><tr><td>24</td><td>ENSBTAG00000004295</td><td>11</td><td>93011815</td><td>93029730</td><td>-1</td><td>NADH dehydrogenase (ubiquinone) 1 alpha subcomplex, 8, 19kDa</td><td>NDUFA8</td><td>144</td></tr><tr><td rowspan="2">27</td><td>ENSBTAG00000004296</td><td>11</td><td>93029874</td><td>93064357</td><td>1</td><td>MORN repeat containing 5</td><td>MORN5</td><td></td></tr><tr><td>ENSBTAG00000004991</td><td>8</td><td>100163708</td><td>100212829</td><td>-1</td><td>Inhibitor of kappa light polypeptide gene enhancer in B-cells, kinase complex-associated protein family with sequence similarity 206, member A</td><td>IKBKAP</td><td>324</td></tr><tr><td rowspan="2">31</td><td>ENSBTAG00000006195</td><td>29</td><td>44691682</td><td>44693901</td><td>-1</td><td>Chromosome 29 open reading frame, human C11orf68</td><td>C29H11orf68</td><td>270</td></tr><tr><td>ENSBTAG00000006199</td><td>29</td><td>44694171</td><td>44696295</td><td>1</td><td>DR1-associated protein 1 (negative cofactor 2 alpha)</td><td>DRAP1</td><td></td></tr><tr><td rowspan="2">61</td><td>ENSBTAG00000023018</td><td>28</td><td>44283041</td><td>44397693</td><td>-1</td><td>Poly(ADP-ribose) glycohydrolase</td><td>PARG</td><td>121</td></tr><tr><td>ENSBTAG00000011694</td><td>28</td><td>44397814</td><td>44418582</td><td>1</td><td>translocase of inner mitochondrial membrane 23 homolog (yeast)</td><td>TIMM23</td><td></td></tr><tr><td rowspan="2">62</td><td>ENSBTAG00000011930</td><td>17</td><td>63448841</td><td>63466323</td><td>-1</td><td>DEAD (Asp-Glu-Ala-Asp) box polypeptide 54</td><td>DDX54</td><td>327</td></tr><tr><td>ENSBTAG00000011937</td><td>17</td><td>63466650</td><td>63472165</td><td>1</td><td>chromosome 17 open reading frame, human C12orf52 (C17H12orf52)</td><td>RITA</td><td></td></tr><tr><td rowspan="2">66</td><td>ENSBTAG00000012586</td><td>2</td><td>86438979</td><td>86449372</td><td>-1</td><td>Heat shock 60kDa protein 1 (chaperonin)</td><td>HSPD1</td><td>153</td></tr><tr><td>ENSBTAG00000012589</td><td>2</td><td>86449525</td><td>86451564</td><td>1</td><td>heat shock 10kDa protein 1 (chaperonin 10)</td><td>HSPE1</td><td></td></tr><tr><td rowspan="2">81</td><td>ENSBTAG00000015667</td><td>11</td><td>2324293</td><td>2335178</td><td>-1</td><td>Transmembrane protein 127</td><td>TMEM127</td><td>204</td></tr><tr><td>ENSBTAG00000015659</td><td>11</td><td>2335382</td><td>2342308</td><td>1</td><td>cytosolic iron-sulfur protein assembly 1</td><td>CIAO1</td><td></td></tr><tr><td rowspan="2">85</td><td>ENSBTAG00000016559</td><td>29</td><td>43146856</td><td>43149460</td><td>-1</td><td>Nudix (nucleoside diphosphate linked moiety X)-type motif 22</td><td>NUDT22</td><td>150</td></tr><tr><td>ENSBTAG00000016555</td><td>29</td><td>43149610</td><td>43152143</td><td>1</td><td>tRNA phosphotransferase 1</td><td>TRPT1</td><td></td></tr><tr><td rowspan="2">86</td><td>ENSBTAG00000016558</td><td>25</td><td>1354553</td><td>1360122</td><td>-1</td><td>splA/ryanodine receptor domain and SOCS box containing 3 nucleotide binding protein 2</td><td>SPSB3</td><td>457</td></tr><tr><td>ENSBTAG00000016561</td><td>25</td><td>1360579</td><td>1365802</td><td>1</td><td>Mitotic spindle organizing protein 1</td><td>NUBP2</td><td></td></tr><tr><td rowspan="2">100</td><td>ENSBTAG00000025028</td><td>12</td><td>47768049</td><td>47781910</td><td>-1</td><td>bora, aurora kinase A activator</td><td>MZT1</td><td>205</td></tr><tr><td>ENSBTAG00000019886</td><td>12</td><td>47782115</td><td>47809879</td><td>1</td><td>WD repeat domain 65</td><td>BORA</td><td></td></tr><tr><td rowspan="2">104</td><td>ENSBTAG00000021547</td><td>3</td><td>103353700</td><td>103422206</td><td>-1</td><td>EBNA1 binding protein 2</td><td>WDR65</td><td>399</td></tr><tr><td>ENSBTAG00000021544</td><td>3</td><td>103422605</td><td>103429520</td><td>1</td><td>Electron-transferring-flavoprotein dehydrogenase</td><td>EBNA1BP2</td><td></td></tr><tr><td rowspan="2">115</td><td>ENSBTAG00000016679</td><td>17</td><td>41262754</td><td>41307126</td><td>-1</td><td>chromosome 17 open reading frame, human C4orf46</td><td>ETFDH</td><td>252</td></tr><tr><td>ENSBTAG00000033486</td><td>17</td><td>41307378</td><td>41310601</td><td>1</td><td>WD repeat domain 54</td><td>WDR54</td><td>90</td></tr><tr><td rowspan="2">118</td><td>ENSBTAG00000000489</td><td>11</td><td>10186102</td><td>10189962</td><td>-1</td><td>chromosome 11 open reading frame, human C2orf81</td><td>C11H2orf81</td><td></td></tr><tr><td>ENSBTAG000000040226</td><td>11</td><td>10190052</td><td>10195932</td><td>1</td><td></td><td></td><td></td></tr></tbody></table>

Distance, distance between transcription start sites, in bp.

(Supplementary Tables S9 and S10). The percentage of co-expressed genes is comparable for the bidirectional, convergent and unidirectional configurations, suggesting that co-expression is most likely mediated by shared chromatin environment rather than only by specific regulatory regions. Similar findings were previously found in other species, including yeast.⁴² In addition, on average the expression levels in muscle of gene pairs sharing a BiP were more correlated

compared with pairs of neighbouring genes in a unidirectional or convergent configuration (Z-score test, P-value <0.05%). Using normalized gene expression counts previously measured in 18 bovine tissues and for three different cows,⁴³,⁴⁴ we calculated the Pearson correlation coefficient between the expression levels of the 120 gene pairs sharing a BiP, the 163 gene pairs in a unidirectional configuration and the 218 gene pairs in a convergent configuration. We then
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Fluorescence microscopy analysis of bidirectional promoter BiP100 in C2C12 cells. Images were taken at x20 magnification, 36h after transfection.</figcaption></figure>
compared the average expression correlation between these three types of gene pairs and found no significant differences between gene pairs in unidirectional and convergent configurations, but identified significant differences between gene pairs sharing a BiP and the other two types of gene pairs (Z-score test, P-value <0.05%), The higher correlation between expression levels of genes sharing BiPs has previously been described and is a common feature of bidirectional promoters.³⁹

### 3.6. The activity of the bidirectional promoters is highly variable

Not all the 240 genes sharing the 120 active bidirectional promoters were expressed in all the twenty samples. Among these genes, 227 genes (corresponding to 114 bidirectional promoters) were expressed in all 20 samples; whereas 13 genes (from 13 different bidirectional promoters) were expressed in all but one sample. In order to evaluate inter-individual variability of the activity of the bidirectional promoters active in muscle, we calculated for each sample the expression ratio of each gene pair (using normalized gene expression levels). We considered that the value of this ratio is a measure of the activity of the bidirectional promoter. We then calculated the coefficient of variation for the activity of each bidirectional promoter. We could calculate the activity (paired-gene ratio) only for 107 bidirectional promoters, as some genes were not expressed in some of the 20 samples. Interestingly there were significant differences when we compared the activity variability of these 107 bidirectional promoters. Thirteen bidirectional promoters had a slightly controled activity (CV <50%) whereas 41 had a very variable activity (CV >100%; Supplementary Table S11). These results suggest a strong or a loosened gene expression regulation of the gene pairs depending on the shared bidirectional promoter. This highly variable activity could be due to polymorphisms within the regulatory elements located within the bidirectional promoter regions.

### 3.7. The expression level of some genes is highly correlated with the activity of some bidirectional promoters active in muscle

Thanks to the availability of the expression levels of genes expressed in LT for 20 samples, it is also possible to analyse the effect of these genes on the activity of the identified bidirectional promoters. Using the normalized expression levels, we calculated the Spearman's rank correlation coefficient for each bidirectional promoter ratio with each gene expressed within the 20 muscle samples. We found after correction for multiple testing 77 significant correlations ($P < 1.9 \times 10^{-8}$), between 77 different genes and the activity of 12 different bidirectional promoters (Supplementary Table S12). We found 74 positive correlations and only 3 negative correlations (anti-correlations). Among these 77 correlations, we could find 9 significant correlations between the activity of a bidirectional promoter and the expression level of one of the two paired-genes sharing the bidirectional promoter. To characterize the functions of the genes showing expression levels highly correlated with the activity of these 12 bidirectional promoters, we performed a GO annotation. No significant enrichment in GO terms and KEGG terms/pathways was found.

Interestingly, with the exception of the paired-genes sharing a bidirectional promoter, no highly correlated genes were in close vicinity of the bidirectional promoter, indicating that the relationship between these genes and the bidirectional promoters are trans-effects. In addition, no genes are associated with more than a single

bidirectional promoter, suggesting that the regulation of bidirectional promoters is highly variable and no master regulator genes may be at play.

We have to be cautious with our results, as a significant correlation does not relate necessarily to an effect of the gene showing a high correlation on the activity of the bidirectional promoter. Indeed, the paired-genes sharing the bidirectional promoter could themselves have an impact on the regulation of the expression of the gene, for which the expression level shows a strong correlation. However, some of the genes showing a strong correlation with the activity of some of the bidirectional promoters encode for proteins involved in gene regulation, such as histones H2a.1 and H3.1 and transcription factors Hnf1a and Zfn133 and therefore might have direct effect on the activity of these bidirectional promoters. However, no binding sites for Hnf1a or Zfn133 were predicted within the sequence of bidirectional promoter #118 or #81, respectively. More work is needed to establish the directionality of the relationships detected with the correlation and the true involvement of identified genes with regulatory functions.

### 3.8. Some genetic variants might alter the activity of the bidirectional promoters active in muscle

After assessing if some genes expressed in muscle influence the activity of the identified bidirectional promoters, we performed a similar analysis to this time to identify DNA polymorphisms altering the activity of these bidirectional promoters.

DNA samples from the twenty animals were genotyped with the Illumina BovineSNP50 Genotyping Beadchips. Among the 54,001 SNPs included in the chip, only polymorphic single nucleotide polymorphisms (SNPs) mapping to bovine autosomes on the bovine genome assembly UMD3.1 version,²⁹ genotyped on ≥90% of the individuals, without a low minor allele frequency (MAF ≥0.05) and in Hardy-Weinberg Equilibrium (P-value <0.001) were retained for further analysis. 21,470 informative SNPs were left after all the filtering steps. We calculated the Spearman's rank correlation coefficient for each bidirectional promoter ratio with the genotypes of each SNP. We found after correction for multiple testing 27 significant correlations (11 and 16 positive and negative correlations respectively, P-value <1.94 × 10⁻⁸), between 27 different SNPs and the activity of 4 different bidirectional promoters. Annotation of these variants indicates that 20 SNPs are intergenic while the remaining 7 SNPs are within genes but none introduce missense changes in the corresponding proteins or are located in putative regulatory regions (Table 3). GO annotation of these seven genes did not reveal any significant enrichment. None of the intergenic or intronic SNPs were located within known genes encoding bovine long non-coding RNAs.²⁶,⁴⁵ The 27 SNPs were neither in close vicinity of the bidirectional promoter and they are all associated with only one bidirectional promoter. This suggest, as seen with the genes whose expression were significantly correlated with the activity of some bidirectional promoters that the effect of these SNPs are trans-effects and there might not be master regulator genes. In addition, annotation of these variants indicates that most SNPs are intergenic and it is unlikely that these 27 SNPs have any direct effect on the activity of the bidirectional promoter.

In order to identify SNPs altering directly the activity of some of the bidirectional promoters active in muscle, we performed another correlation analysis between the activity of the bidirectional promoters and this time the genotypes from coding SNPs identified with the RNA-Seq data. Mapped paired-end reads from each sample were
<table><caption>Table 3. List of SNPs from the Illumina BovineSNP50 Beadchip with genotypes highly correlated to the activity of some bidirectional promoters</caption><thead><tr><th>SNP ID</th><th>BTA</th><th>Position</th><th>Alleles</th><th>Consequence</th><th>Ensembl gene ID</th><th>Ensembl transcript ID</th><th>Gene symbol</th><th>BiP #</th><th>Rho</th><th>P-value</th></tr></thead><tbody><tr><td>Hapmap42177-BTA-31679</td><td>1</td><td>61,125,554</td><td>A/G</td><td>intergenicleder</td><td></td><td></td><td></td><td>118</td><td>0.99</td><td>1.78E-24</td></tr><tr><td>Hapmap42933-BTA-40454</td><td>1</td><td>85,731,453</td><td>A/G</td><td>intergenicleder</td><td></td><td></td><td></td><td>101</td><td>-0.99</td><td>1.95E-20</td></tr><tr><td>ARS-BFGL-NGS-108514</td><td>2</td><td>105,521,127</td><td>A/G</td><td>intergenicleder</td><td></td><td></td><td></td><td>118</td><td>-0.99</td><td>1.78E-24</td></tr><tr><td>Hapmap54229-rs29017613</td><td>4</td><td>30,200,987</td><td>A/G</td><td>intergenicleder</td><td></td><td></td><td></td><td>118</td><td>0.99</td><td>1.78E-24</td></tr><tr><td>ARS-BFGL-NGS-35869</td><td>5</td><td>21,707,809</td><td>A/G</td><td>intergenicleder</td><td></td><td></td><td></td><td>118</td><td>-0.99</td><td>1.78E-24</td></tr><tr><td>ARS-BFGL-NGS-93953</td><td>5</td><td>22,737,219</td><td>A/G</td><td>intergenicleder</td><td></td><td></td><td></td><td>118</td><td>0.99</td><td>1.78E-24</td></tr><tr><td>Hapmap49852-BTA-107572</td><td>6</td><td>29,790,005</td><td>A/G</td><td>intergenicleder</td><td></td><td></td><td></td><td>104</td><td>-0.99</td><td>1.78E-24</td></tr><tr><td>Hapmap25168-BTC-033275</td><td>6</td><td>33,713,818</td><td>A/G</td><td>intergenicleder</td><td></td><td></td><td></td><td>104</td><td>0.99</td><td>1.78E-24</td></tr><tr><td>Hapmap36813-</td><td>6</td><td>33,768,128</td><td>A/G</td><td>intergenicleder</td><td></td><td></td><td></td><td>104</td><td>0.99</td><td>1.78E-24</td></tr><tr><td>SCAFFOLD50174_9004</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Hapmap23923-BTC-066021</td><td>6</td><td>39,721,727</td><td>A/C</td><td>intergenicleder</td><td></td><td></td><td></td><td>104</td><td>-0.99</td><td>1.78E-24</td></tr><tr><td>Hapmap36286-</td><td>7</td><td>111,161,115</td><td>A/C</td><td>intergenicleder</td><td></td><td></td><td></td><td>104</td><td>-0.99</td><td>1.78E-24</td></tr><tr><td>SCAFFOLD260285_24265</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>ARS-BFGL-NGS-99031</td><td>8</td><td>101,251,865</td><td>A/C</td><td>intergenicleder</td><td></td><td></td><td></td><td>118</td><td>0.99</td><td>1.78E-24</td></tr><tr><td>ARS-BFGL-NGS-60678</td><td>9</td><td>85,454,475</td><td>A/G</td><td>intergenicleder</td><td></td><td></td><td></td><td>104</td><td>0.99</td><td>1.78E-24</td></tr><tr><td>ARS-BFGL-NGS-4488</td><td>10</td><td>6,476,252</td><td>A/G</td><td>intronleder</td><td></td><td></td><td></td><td>104</td><td>-1</td><td>0</td></tr><tr><td>ARS-BFGL-NGS-5819</td><td>10</td><td>10,428,184</td><td>A/G</td><td>intronleder</td><td></td><td></td><td></td><td>118</td><td>0.99</td><td>1.78E-24</td></tr><tr><td>ARS-BFGL-NGS-27341</td><td>11</td><td>91,438,914</td><td>A/G</td><td>intergenicleder</td><td></td><td></td><td></td><td>118</td><td>-0.99</td><td>1.78E-24</td></tr><tr><td>ARS-BFGL-NGS-115889</td><td>12</td><td>79,814,959</td><td>A/C</td><td>intronleder</td><td></td><td></td><td></td><td>64</td><td>0.99</td><td>1.95E-20</td></tr><tr><td>ARS-BFGL-NGS-23509</td><td>14</td><td>70,636,087</td><td>A/G</td><td>intergenicleder</td><td></td><td></td><td></td><td>118</td><td>-0.99</td><td>1.78E-24</td></tr><tr><td>ARS-BFGL-NGS-859</td><td>14</td><td>73,919,098</td><td>A/G</td><td>intergenicleder</td><td></td><td></td><td></td><td>104</td><td>-0.99</td><td>1.78E-24</td></tr><tr><td>ARS-BFGL-NGS-99802</td><td>16</td><td>74,999,809</td><td>A/C</td><td>5-primeleder</td><td></td><td></td><td></td><td>104</td><td>0.99</td><td>1.78E-24</td></tr><tr><td>Hapmap52466-rs29015577</td><td>19</td><td>3,617,183</td><td>A/C</td><td>intergenicleder</td><td></td><td></td><td></td><td>118</td><td>-0.99</td><td>1.78E-24</td></tr><tr><td>Hapmap47625-BTA-44726</td><td>19</td><td>21,878,635</td><td>A/G</td><td>intergenicleder</td><td></td><td></td><td></td><td>104</td><td>-0.99</td><td>1.78E-24</td></tr><tr><td>ARS-BFGL-NGS-14187</td><td>19</td><td>25,165,920</td><td>A/G</td><td>intronleder</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>ARS-BFGL-NGS-109291</td><td>23</td><td>13,517,193</td><td>T/A</td><td>intronleder</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>ARS-BFGL-NGS-57958</td><td>26</td><td>23,000,155</td><td>A/G</td><td>intergenicleder</td><td></td><td></td><td></td><td>104</td><td>-0.99</td><td>1.78E-24</td></tr><tr><td>ARS-BFGL-NGS-42033</td><td>28</td><td>14,243,381</td><td>C/G</td><td>intergenicleder</td><td></td><td></td><td></td><td>118</td><td>-0.99</td><td>1.78E-24</td></tr><tr><td>ARS-BFGL-NGS-77028</td><td>28</td><td>26,689,199</td><td>A/G</td><td>intronleder</td><td></td><td></td><td></td><td>118</td><td>0.99</td><td>1.78E-24</td></tr></tbody></table>
<table><caption>**Table 4.** List of coding SNPs with genotypes highly correlated to the activity of some bidirectional promoters</caption><thead><tr><th>SNP ID</th><th>BTA</th><th>Position</th><th>dbSNP ID</th><th>Alleles</th><th>Consequence</th><th>Ensembl gene ID</th><th>Ensembl transcript ID</th><th>Gene symbol</th><th>BiP #</th><th>Rho</th><th>P-value</th></tr></thead><tbody><tr><td>ENSBTAT0000002157_515</td><td>9</td><td>61,291,956</td><td></td><td>G/A</td><td>intronled variant</td><td>ENSBTAG00000001644</td><td>ENSBTAT00000002157</td><td>MDN1</td><td>104</td><td>-0.99</td><td>1.78E-24</td></tr><tr><td>ENSBTAT0000003078_106</td><td>15</td><td>82,337,162</td><td></td><td>G/A</td><td>intronled variant</td><td>ENSBTAG00000002381</td><td>ENSBTAT00000003078</td><td>ZDHHC5</td><td>118</td><td>0.99</td><td>1.78E-24</td></tr><tr><td>ENSBTAT0000003826_479</td><td>21</td><td>22,218,555</td><td></td><td>C/T</td><td>upstreamled geneled variant</td><td>ENSBTAG00000002939</td><td>ENSBTAT00000003826</td><td>FURIN</td><td>104</td><td>-0.99</td><td>1.78E-24</td></tr><tr><td>ENSBTAT0000005923_585</td><td>2.5</td><td>2,977,744</td><td></td><td>T/C</td><td>upstreamled geneled variant</td><td>ENSBTAG00000004509</td><td>ENSBTAT00000005923</td><td>SLX4</td><td>118</td><td>-0.99</td><td>1.78E-24</td></tr><tr><td>ENSBTAT0000008129_240</td><td>16</td><td>19,500,611</td><td></td><td>G/C</td><td>intronled variant</td><td>ENSBTAG00000006186</td><td>ENSBTAT00000008129</td><td>KCTD3</td><td>118</td><td>0.99</td><td>1.78E-24</td></tr><tr><td>ENSBTAT0000008507_123</td><td>15</td><td>82,265,194</td><td></td><td>T/C</td><td>intronled variant</td><td>ENSBTAG00000006493</td><td>ENSBTAT00000008507</td><td>CLP1</td><td>118</td><td>-0.99</td><td>1.78E-24</td></tr><tr><td>ENSBTAT0000008728_747</td><td>12</td><td>18,197,212</td><td></td><td>C/T</td><td>intronled variant</td><td>ENSBTAG00000006640</td><td>ENSBTAT00000008728</td><td>RBI</td><td>118</td><td>0.99</td><td>1.78E-24</td></tr><tr><td>ENSBTAT00000018163_665</td><td>21</td><td>66,866,347</td><td></td><td>T/G</td><td>upstreamled geneled variant</td><td>ENSBTAG00000013666</td><td>ENSBTAT00000018163</td><td>SLC25A29</td><td>104</td><td>-0.99</td><td>1.78E-24</td></tr><tr><td>ENSBTAT00000020493_130</td><td>7</td><td>62,925,278</td><td></td><td>A/G</td><td>intronled variant</td><td>ENSBTAG00000015419</td><td>ENSBTAT00000020493</td><td>ARHGEF37</td><td>64</td><td>-0.99</td><td>1.95E-20</td></tr><tr><td>ENSBTAT00000025492_402</td><td>19</td><td>55,675,666</td><td></td><td>T/C</td><td>intronled variant</td><td>ENSBTAG00000019153</td><td>ENSBTAT00000025492</td><td>JMJD6</td><td>24</td><td>-1</td><td>6.38E-137</td></tr><tr><td>ENSBTAT00000028060_520</td><td>21</td><td>16,720,681</td><td></td><td>C/T</td><td>intronled variant</td><td>ENSBTAG00000037383</td><td>ENSBTAT00000028060</td><td>AKAP13</td><td>118</td><td>-0.99</td><td>1.78E-24</td></tr><tr><td>ENSBTAT00000028277_365</td><td>15</td><td>76,757,181</td><td></td><td>T/G</td><td>intronled variant</td><td>ENSBTA00000021223</td><td>ENSBTAT00000028277</td><td>CRY2</td><td>118</td><td>-0.99</td><td>1.78E-24</td></tr><tr><td>ENSBTAT00000028656_569</td><td>14</td><td>75,708,914</td><td></td><td>T/G</td><td>intronled variant</td><td>ENSBTAG00000039968</td><td>ENSBTAT00000028656</td><td>TMEM55A</td><td>104</td><td>0.99</td><td>1.78E-24</td></tr><tr><td>ENSBTAT00000028865_103</td><td>2.5</td><td>39,432,034</td><td></td><td>G/A</td><td>intronled variant</td><td>ENSBTAG00000037400</td><td>ENSBTAT00000028865</td><td>TNRC18</td><td>118</td><td>0.99</td><td>1.78E-24</td></tr><tr><td>ENSBTAT00000029400_195</td><td>17</td><td>73,658,214</td><td></td><td>C/T</td><td>upstreamled geneled variant</td><td>ENSBTAG00000021656</td><td>ENSBTAT00000029400</td><td>SPECC1L</td><td>104</td><td>0.99</td><td>1.78E-24</td></tr><tr><td>ENSBTAT00000029403_203</td><td>18</td><td>62,789,937</td><td></td><td>A/G</td><td>intronled variant</td><td>ENSBTAG00000030393</td><td>ENSBTAT00000029403</td><td>RDH13</td><td>104</td><td>-0.99</td><td>1.78E-24</td></tr><tr><td>ENSBTAT00000033704_233</td><td>6</td><td>113,701,647</td><td>rs442770236</td><td>G/T</td><td>upstreamled geneled variant</td><td>ENSBTAG00000004316</td><td>ENSBTAT00000033704</td><td>BOD1L</td><td>104</td><td>0.99</td><td>1.78E-24</td></tr><tr><td>ENSBTAT00000035362_466</td><td>20</td><td>10,305,279</td><td></td><td>G/A</td><td>downstreamled geneled variant</td><td>ENSBTAG000000027980</td><td>ENSBTAT00000035362</td><td>TAF9</td><td>118</td><td>-0.99</td><td>1.78E-24</td></tr><tr><td>ENSBTAT00000037465_138</td><td>25</td><td>10,053,612</td><td></td><td>T/C</td><td>intronled variant</td><td>ENSBTAG00000026375</td><td>ENSBTAT00000037465</td><td>RM12</td><td>118</td><td>-0.99</td><td>1.78E-24</td></tr><tr><td>ENSBTAT00000043122_702</td><td>10</td><td>86,384,804</td><td></td><td>A/G</td><td>upstreamled geneled variant</td><td>ENSBTAG00000020379</td><td>ENSBTAT00000043122</td><td>AREL1</td><td>118</td><td>0.99</td><td>1.78E-24</td></tr><tr><td>ENSBTAT00000043260_900</td><td>10</td><td>80,246,761</td><td></td><td>A/G</td><td>upstreamled geneled variant</td><td>ENSBTAG00000014334</td><td>ENSBTAT00000043260</td><td>ZFYVE26</td><td>118</td><td>0.99</td><td>1.78E-24</td></tr><tr><td>ENSBTAT00000043778_165</td><td>15</td><td>82,285,586</td><td></td><td>T/C</td><td>intronled variant</td><td>ENSBTAG00000002411</td><td>ENSBTAT00000043778</td><td>CTNNND1</td><td>118</td><td>0.99</td><td>1.78E-24</td></tr><tr><td>ENSBTAT00000043778_492</td><td>15</td><td>82,284,426</td><td></td><td>G/A</td><td>intronled variant</td><td>ENSBTAG00000002411</td><td>ENSBTAT00000043778</td><td>CTNNND1</td><td>118</td><td>-0.99</td><td>1.78E-24</td></tr><tr><td>ENSBTAT00000043778_760</td><td>15</td><td>82,284,694</td><td></td><td>C/T</td><td>intronled variant</td><td>ENSBTAG00000002411</td><td>ENSBTAT00000043778</td><td>CTNNND1</td><td>118</td><td>-0.99</td><td>1.78E-24</td></tr><tr><td>ENSBTAT00000054096_509</td><td>8</td><td>103,426,705</td><td></td><td>G/A</td><td>upstreamled geneled variant</td><td>ENSBTAG00000038335</td><td>ENSBTAT00000054096</td><td>IGBP1</td><td>118</td><td>0.99</td><td>1.78E-24</td></tr><tr><td>ENSBTAT00000056520_352</td><td>2.3</td><td>27,378,739</td><td></td><td>C/T</td><td>upstreamled geneled variant</td><td>ENSBTAG00000039620</td><td>ENSBTAT00000056520</td><td>MGC151586</td><td>104</td><td>0.99</td><td>1.78E-24</td></tr><tr><td>ENSBTAT00000061199_185</td><td>14</td><td>20,989,316</td><td></td><td>A/G</td><td>upstreamled geneled variant</td><td>ENSBTAG00000044106</td><td>ENSBTAT00000061199</td><td>SPIDR</td><td>118</td><td>-0.99</td><td>1.78E-24</td></tr><tr><td>ENSBTAT00000061451_269</td><td>2.3</td><td>3,434,098</td><td></td><td>G/A</td><td>intronled variant</td><td>ENSBTAG00000021237</td><td>ENSBTAT00000061451</td><td>DST</td><td>118</td><td>-0.99</td><td>1.78E-24</td></tr><tr><td>ENSBTAT00000061451_704</td><td>23</td><td>3,432,111</td><td></td><td>G/T</td><td>intronled variant</td><td>ENSBTAG00000021237</td><td>ENSBTAT00000061451</td><td>DST</td><td>118</td><td>0.99</td><td>1.78E-24</td></tr><tr><td>ENSBTAT00000034195_104</td><td>8</td><td>53,736,890</td><td></td><td>G/A</td><td>intronled variant</td><td>ENSBTAG00000017734</td><td>ENSBTAT00000034195</td><td>VP513A</td><td>118</td><td>0.99</td><td>1.78E-24</td></tr></tbody></table>
used with the SAMtools package for variant calling. Using stringent parameters (e.g. minimum coverage of 8 reads and mapping quality of 20) we detected 8,121 different billalelic SNPs. Among the SNPs identified, similarly to the SNPs from the Illumina BovineSNP50 Genotyping Beadchip, we used only polymorphisms mapping to bovine autosomes, genotyped on ≥90% of the individuals, without a low minor allele frequency (MAF ≥0.05) and in Hardy-Weinberg Equilibrium (P-value <0.001). A total of 8,116 informative SNPs were left after these filtering steps. We calculated the Spearman's rank correlation coefficient for each bidirectional promoter ratio with the genotypes of each SNP. We found after correction for multiple testing 30 significant correlations (P-value <5.13 × 10⁻⁸), between 30 different SNPs (in 27 different genes) and the activity of 4 different bidirectional promoters (Table 4). Annotation of these variants indicates that all these SNPs are within protein-coding genes. GO annotation of the 27 genes containing the 30 coding SNPs (cSNPs) did not reveal any significant enrichment. None of the 30 cSNPs introduce missense changes in the corresponding proteins; however, 11 cSNPs are located in the upstream or downstream regions.

We found among these correlations a significant negative correlation between a cSNP within *Taf9*, a gene encoding the RNA polymerase II, TATA box binding protein-associated factor (TAF) 32kD (also known as *TafII32* or *Taf2g*), and the activity of bidirectional promoter #118. TAF9 is involved in the initiation of the transcription by RNA polymerase II.⁴⁶ Interestingly, Taf1, another TATA box binding factor has motifs over-represented in human bidirectional promoters⁴⁷ and two-hybrid experiments performed in yeast have shown that Taf1 is able to bind to DNA but also to different TATA box binding factors, including TAF9.⁴⁸ Many bidirectional promoters lack a TATA box;³,¹¹ however, it has been shown that several TATA-less promoters require TAFs for transcription.⁴⁷⁻⁵³ The cSNP found in *Taf9* is therefore an interesting variant and might point out an effect of Taf9 on the activity of bidirectional promoter #118.

The 11 cSNPs might have an impact on the regulation of the corresponding gene, which in turn might affect the activity of some bidirectional promoters. Spearman's rank correlation coefficients were calculated for the genotypes of these 11 cSNPs and the normalized expression level of the corresponding gene (Table 5). Significant correlations were found at P-value <0.05 for 4 cSNPs, suggesting that these cSNPs have an effect on the regulation of the expression of their corresponding gene. The seven other cSNPs, including the variant found in *Taf9*, might have most probably an effect on the post-transcriptional regulation of the corresponding gene, which consequently alter the amount of proteins which then affect the activity of the bidirectional promoter.

### 3.9. Some bidirectional promoter regions are located within quantitative trait loci for meat-related phenotypes

The positions of the 120 bidirectional promoters sharing genes expressed in bovine LT were compared with the position of known quantitative trait loci (QTLs) deposited in the public database AnimalQTLdb.⁵⁴ Fifty-four bidirectional promoters were located in 244 different QTL regions (Supplementary Table S14). For example, 7 different bidirectional promoters are found in 8 QTL regions for intramuscular fat; whereas 18 bidirectional promoters are within 18 QTLs for shear force. Karim et al.⁵⁵ have shown that a genetic variant located within a bidirectional promoter shared by *PLAG1* and

<table><caption>**Table 5. Correlation scores between cSNPs and the expression level of the corresponding gene**</caption><thead><tr><th>cSNP ID</th><th>BTA</th><th>Position</th><th>dbSNP ID</th><th>Alleles</th><th>Consequence</th><th>Ensembl transcript ID</th><th>Gene symbol</th><th>Rho</th><th>P-value</th></tr></thead><tbody><tr><td>ENSBTAT0000003826_4794</td><td>21</td><td>22,218,555</td><td></td><td>C/T</td><td>upstream_geneleder</td><td>ENSBTAT00000003826</td><td>FURIN</td><td>-0.52</td><td>0.02</td></tr><tr><td>ENSBTAT00000018163_665</td><td>21</td><td>66,866,347</td><td></td><td>T/G</td><td>upstream_geneleder</td><td>ENSBTAT00000018163</td><td>SLC25A29</td><td>-0.49</td><td>0.03</td></tr><tr><td>ENSBTAT00000043122_702</td><td>10</td><td>86,384,804</td><td></td><td>A/G</td><td>upstream_geneleder</td><td>ENSBTAT00000043122</td><td>AREL1</td><td>0.49</td><td>0.03</td></tr><tr><td>ENSBTAT00000043260_900</td><td>10</td><td>80,246,761</td><td></td><td>A/G</td><td>upstream_geneleder</td><td>ENSBTAT00000043260</td><td>ZFYVE26</td><td>0.46</td><td>0.04</td></tr><tr><td>ENSBTAT0000005923_585</td><td>25</td><td>2,977,744</td><td></td><td>T/C</td><td>upstream_geneleder</td><td>ENSBTAT0000005923</td><td>SLX4</td><td>-0.43</td><td>0.06</td></tr><tr><td>ENSBTAT00000035362_4667</td><td>20</td><td>10,305,279</td><td></td><td>G/A</td><td>downstream_geneleder</td><td>ENSBTAT00000035362</td><td>TAF9</td><td>-0.43</td><td>0.06</td></tr><tr><td>ENSBTAT00000029400_195</td><td>17</td><td>73,658,214</td><td></td><td>C/T</td><td>upstream_geneleder</td><td>ENSBTAT00000029400</td><td>SPECC1L</td><td>0.35</td><td>0.14</td></tr><tr><td>ENSBTAT00000061199_1851</td><td>14</td><td>20,989,316</td><td></td><td>A/G</td><td>upstream_geneleder</td><td>ENSBTAT00000061199</td><td>SPIDR</td><td>-0.29</td><td>0.22</td></tr><tr><td>ENSBTAT00000033704_233</td><td>6</td><td>113,701,647</td><td>rs442770236</td><td>G/T</td><td>upstream_geneleder</td><td>ENSBTAT00000033704</td><td>BOD1L</td><td>0.17</td><td>0.46</td></tr><tr><td>ENSBTAT00000054096_509</td><td>8</td><td>103,426,705</td><td></td><td>G/A</td><td>upstream_geneleder</td><td>ENSBTAT00000054096</td><td>IGBP1</td><td>0.14</td><td>0.54</td></tr><tr><td>ENSBTAT00000056520_3528</td><td>23</td><td>27,378,739</td><td></td><td>C/T</td><td>upstream_geneleder</td><td>ENSBTAT00000056520</td><td>MGC151586</td><td>-0.12</td><td>0.63</td></tr></tbody></table>
CHCHD7 influences bovine stature. It will be interesting to investi-
gate if variants within some of the bidirectional promoter regions we
identified affect quantitative traits, including meat quality traits.

4. Conclusion

In this study, we identified 120 bidirectional promoters active in bo-
vine muscle and validated *in vitro* 16 of them. Using RNA-seq data
from 20 muscle samples, we found that the activity of these bidirec-
tional promoters is highly variable and that the expression level of
77 genes is highly correlated with the activity of 12 of the bidirec-
tional promoters. We also identified 57 SNPs associated significantly
with the activity of 5 bidirectional promoters. Interestingly, we found
that a coding SNP within *Taf9* has an effect on the activity of bidirec-
tional promoter #118. To our knowledge, our study is the first ana-
lysis in any species of the genetic variability of the activity of
bidirectional promoters.

Data availability

The sequencing data have been submitted to the European
Nucleotide Archive (accession numbers ERP002220, E-MTAB-2646
and E-MTAB-4625).

Acknowledgements

We are grateful to Sébastien Legardinier (UMR1061, Limoges) for the gift of the C2C12 cells and Jérôme Chapuis (UR892 INRA VIM, Jouy-en-Josas) for his help with the fluorescence microscope.

Conflict of interest

None declared.

**Supplementary data**

Supplementary data are available at www.dnaresearch.oxfordjournals.org.

Funding

The work was supported by the French National Research Agency (Regulomix project, contract ANR-09-GENM-011), French National Institute for Agricultural Research and the Limousin Regional Council.

References

1. Adachi, N. and Lieber, M. 2002, Bidirectional gene organization: a common architectural feature of the human genome. *Cell* **109**, 807–9.

2. Takai, D. and Jones, P.A. 2004, Origins of bidirectional promoters: com-
putational analyses of intergenic distance in the human genome. *Mol.
Biol. Evol.* **21**, 463–7.

3. Trinklein, N.D., Aldred, S.F., Hartman, S.J., Schroeder, D.L., Otillar, R.P. and Myers, R.M. 2004, An abundance of bidirectional promoters in the human genome. *Genome Res.* **14**, 62–6.

4. Williams, T.J. and Fried, M. 1986, The MES-1 murine enhancer element is closely associated with the heterogeneous 5' ends of two divergent transcription units. *Mol. Cell. Biol.* **6**, 4558–69.

5. Herr, D.R. and Harris, G.L. 2004, Close head-to-head juxtaposition of genes favors their coordinate regulation in *Drosophila melanogaster*. *FEBS Lett.* **572**, 147–53.

6. Koyanagi, K.O., Hagiwara, M., Itoh, T., Gojobori, T. and Imanishi, T. 2005, Comparative genomics of bidirectional gene pairs and its implications for the evolution of a transcriptional regulation system. *Gene* **353**, 169–76.

7. Engstrom, P.G., Suzuki, H., Ninomiya, N., et al. 2006, Complex loci in human and mouse genomes. *PLoS Genet.* **2**, e47.

8. Li, Y.Y., Yu, H., Guo, Z.M., Guo, T.Q., Tu, K. and Li, X.Y. 2006, Systematic analysis of head-to-head gene organization: evolutionary conservation and potential biological relevance. *PLoS Comput. Biol.* **2**, e74.

9. Krom, N. and Ramakrishna, W. 2008, Comparative analysis of divergent and convergent gene pairs and their expression patterns in rice, *Arabidopsis*, and populus. *Plant Physiol.* **147**, 1763–73.

10. Yang, M.Q. and Elnitski, L.L. 2008, Prediction-based approaches to char-
acterize bidirectional promoters in the mammalian genome. BMC
Genomics 9 (Suppl 1), S2.

11. Yang, M.Q. and Elnitski, L.L. 2008, Diversity of core promoter elements comprising human bidirectional promoters. *BMC Genomics* 9 (Suppl 2), S3.

12. Yang, M.Q., Taylor, J. and Elnitski, L. 2008, Comparative analyses of bi-directional promoters in vertebrates. *BMC Bioinformatics* 9 (Suppl 6), S9.

13. Piontkivska, H., Yang, M.Q., Larkin, D.M., Lewin, H.A., Reecy, J. and Elnitski, L. 2009, Cross-species mapping of bidirectional promoters enables prediction of unannotated 5' UTRs and identification of species-specific transcripts. *BMC Genomics* **10**, 189.

14. Wang, Q., Wan, L., Li, D., Zhu, L., Qian, M. and Deng, M. 2009, Searching for bidirectional promoters in *Arabidopsis thaliana*. *BMC Bioinformatics* **10** (Suppl 1), S29.

15. Krom, N. and Ramakrishna, W. 2010, Conservation, rearrangement, and deletion of gene pairs during the evolution of four grass genomes. *DNA Res.* **17**, 343–52.

16. Woo, Y.H. and Li, W.H. 2011, Gene clustering pattern, promoter architecture, and gene expression stability in eukaryotic genomes. *Proc. Natl. Acad. Sci. U.S.A.* **108**, 3306–11.

17. Behura, S.K. and Severson, D.W. 2015, Bidirectional promoters of insects: genome-wide comparison, evolutionary implication and influence on gene expression. *J. Mol. Biol.* **427**, 521–36.

18. Maxson, R., Cohn, R., Kedes, L. and Mohun, T. 1983, Expression and organization of histone genes. *Annu. Rev. Genet.* **17**, 239–77.

19. Momota, R., Sugimoto, M., Oohashi, T., Kigasawa, K., Yoshioka, H. and Ninomiya, Y. 1998, Two genes, COL4A3 and COL4A4 coding for the human z3(IV) and z4(IV) collagen chains are arranged head-to-head on chromosome 2q36. *FEBS Lett.* **424**, 11–6.

20. Harrison, L., Ascione, A.G., Takiguchi, Y., Wilson, D.M., Chen, D.J. and Temple, B. 1997, Comparison of the promoters of the mouse (APEX) and human (APE) apurinic endonuclease genes. *Mutat. Res.* **385**, 159–72.

21. Tan, M., Wei, C. and Price, C.M. 2003, The telomeric protein Rap1 is conserved in vertebrates and is expressed from a bidirectional promoter positioned between the Rap1 and KARS genes. *Gene* **323**, 1–10.

22. Uwanogho, D.A., Yasin, S.A., Starling, B. and Price, J. 2010, The inter-
genic region between the mouse Recql4 and Lrrc14 genes functions as an
evolutionary conserved bidirectional promoter. *Gene* **449**, 103–17.

23. Meersseman, C., Lejard, V., Rebours, E., et al. 2014, Bovine TWINKLE and mitochondrial ribosomal protein L43 genes are regulated by an evolutionary conserved bidirectional promoter. *Gene* **537**, 154–63.

24. Allais, S., Leveziel, H., Hocquette, J.F., et al. 2014, Fine mapping of quantitative trait loci underlying sensory meat quality traits in three French beef cattle breeds. *J. Anim. Sci.* **92**, 4329–41.

25. Djari, A., Esquerre, D., Weiss, B., et al. 2013, Gene-based single nucleo-
tide polymorphism discovery in bovine muscle using next-generation tran-
scriptomic sequencing. BMC Genomics **14**, 307.

26. Billerey, C., Boussaha, M., Esquerre, D., et al. 2014, Identification of large intergenic non-coding RNAs in bovine muscle using next-generation transcriptomic sequencing. BMC Genomics **15**, 499.

27. Li, H. and Durbin, R. 2009, Fast and accurate short read alignment with Burrows-Wheeler transform. *Bioinformatics* **25**, 175417–60.

28. Li, H., Handsaker, B., Wysoker, A., et al. 2009, The Sequence alignment/
map (SAM) format and SAMtools. *Bioinformatics* **25**, 2078–9.
29. Zimin, A.V., Delcher, A.L., Florea, L., et al. 2009, A whole-genome assembly of the domestic cow. *Bos taurus*. *Genome Biol.* **10**, R42.

30. Li, H; 2011. A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. *Bioinformatics* **27**, 2987-93.

31. Birney, E., Andrews, D., Caccamo, M., et al. 2006, Ensembl 2006. *Nucleic Acids Res.* **34** (Database issue), D556-D561.

32. Al-Shahrour, F., Minguez, P., Vaquerizas, J.M., Conde, L. and Dopazo, J. 2005, BABELOMICS: a suite of web tools for functional annotation and analysis of groups of genes in high-throughput experiments. *Nucleic Acids Res.* **33** (Web Server issue), W460-4.

33. Kel, A.E., Gössling, E., Reuter, I., Cheremushkin, E., Kel-Margoulis, O.V. and Wingender E. 2003, MATCH: a tool for searching transcription factor binding sites in DNA sequences. *Nucleic Acids Res.* **31**, 3576-9.

34. Matys, V., Fricke, E., Geffers, R., et al. 2003, TRANSFAC: transcriptional regulation, from patterns to profiles. *Nucleic Acids Res.* **31**, 374-378.

35. Homepage Repeat Masker Server, 2013. http://www.repeatmasker.org.

36. Rozen, S. and Skaletsky, H. 2000, Primer3 on the WWW for general users and for biologist programmers. *Methods Mol. Biol.* **132**, 365-86.

37. Lejard, V., Rebours, E., Meersseman, C. and Rocha, D. 2014, Construction and validation of a novel dual reporter vector for studying mammalian bidirectional promoters. *Plasmid* **74**, 1-8.

38. Chitwood, J.L., Rincon, G., Kaiser, G.G., Medrano, J.F. and Ross, P.J. 2003, RNA-seq analysis of single bovine blastocysts. *BMC Genomics* **14**, 350.

39. Xu, C., Chen, J. and Shen, B. 2012, The preservation of bidirectional promoter architecture in eukaryotes: what is the driving force? *BMC Syst. Biol.* **6** (Suppl 1), S21.

40. Liu, B., Chen, J. and Shen, B. 2011, Genome-Wide analysis of the transcription factor binding preference of human bi-directional promoters and functional annotation of related gene pairs. *BMC Syst. Biol.* **5** (Suppl 1), S2.

41. Karolchik, D., Baertsch, R., Diekhans, M., et al. 2003, The UCSC genome browser database. *Nucleic Acids Res.* **31**, 51-54.

42. Tsai, H.K., Su, C.P., Lu, M.Y., Shih, C.H and Wang, D. 2007, Co-expression of adjacent genes in yeast cannot be simply attributed to shared regulatory system. *BMC Genomics* **8**, 352.

43. Koufariotis, L.T., Chen, Y.P., Chamberlain, A., Vander Jagt, C. and Hayes, B.J. 2015. A catalogue of novel bovine long noncoding RNA across 18 tissues. *PLoS One* **10**, e0141225.

44. Chamberlain, A.J., Vander Jagt, C.J., Hayes, B.J., et al. 2015. Extensive variation between tissues in allele specific expression in an outbred mammal. *BMC Genomics* **16**, 993.

45. Liu, C., Bai, B., Skogerbø, G., et al. 2005, NONCODE: an integrated knowledge database of non-coding RNAs. *Nucleic Acids Res.* **33** (Database issue), D112-5.

46. Hisatake, K., Ohta, T., Takada, R., et al. 1995, Evolutionary conservation of human TATA-binding-polypeptide-associated factors TAFII31 and TAFII80 and interactions of TAFII80 with other TAFs and with general transcription factors. *Proc. Natl. Acad. Sci. U.S.A.* **92**, 8195-9.

47. Chen, Y., Li, Y., Wei, J. and Li, Y.Y. 2014, Transcriptional regulation and spatial interactions of head-to-head genes. *BMC Genomics* **15**, 519.

48. Yatherajam, G., Zhang, L., Kraemer, S.M. and Stargell, L.A. 2003, Protein-protein interaction map of yeast TFIID. *Nucleic Acids Res.* **31**, 1252-60.

49. Pugh, B.F. and Tjian, R. 1991, Transcription from a TATA-less promoter requires a multisubunit TFIID complex. *Genes Dev.* **5**, 1935-45.

50. Zhou, Q., Lieberman, P.M., Boyer, T.G. and Berk, A.J. 1992, Holo-TFIID supports transcriptional stimulation by diverse activators and from a TATA-less promoter. *Genes Dev.* **6**, 1964-74.

51. Sadowski, C.L., Henry, R.W., Lobo, S.M. and Hernandez, N. 1993, Targeting TBP to a non-TATA box cis-regulatory element: a TBP-containing complex activates transcription from snRNA promoters through the PSE. *Genes Dev.* **7**, 1535-48.

52. Zhou, Q., Boyer, T.G. and Berk, A.J. 1993, Factors (TAFs) required for activated transcription interact with TATA box-binding protein conserved core domain. *Genes Dev.* **7**, 180-7.

53. Martinez, E., Chiang, C.M., Ge, H. and Roeder, R.G. 1994, TATA-binding protein-associated factor(s) in TFIID function through the initiator to direct basal transcription from a TATA-less class II promoter. *EMBO J.* **13**, 3115-312.

54. Hu, Z.L., Fritz, E.R. and Reecy, J.M. 2007, AnimalQTLdb: a livestock QTL database tool set for positional QTL information mining and beyond. *Nucleic Acids Res.* **35** (Database issue), D604-9.

55. Karim, L., Takeda, H., Lin, L., et al. 2011, Variants modulating the expression of a chromosome domain encompassing PLAG1 influence bovine stature. *Nat. Genet.* **43**, 405-13.