# Immunome differences between porcine ileal and jejunal Peyer's patches revealed by global transcriptome sequencing of gut-associated lymphoid tissues

T. Maroilley, M. Berri, G. Lemonnier, D. Esquerré, C. Chevaleyre, S. Melo, F. Meurens, J. Coville, J. Leplat, Andrea Rau, et al.

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

## To cite this version:

T. Maroilley, M. Berri, G. Lemonnier, D. Esquerré, C. Chevaleyre, et al., Immunome differences between porcine ileal and jejunal Peyer's patches revealed by global transcriptome sequencing of gut-associated lymphoid tissues. Scientific Reports, 2018, 8 (1), pp.1-12. 10.1038/s41598-018-27019-7. hal-02356215

HAL Id: hal-02356215

https://hal.science/hal-02356215

Submitted on 26 May 2020

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN

Received: 17 November 2017

Accepted: 18 May 2018

Published online: 13 June 2018

# Immunome differences between porcine ileal and jejunal Peyer's patches revealed by global transcriptome sequencing of gut-associated lymphoid tissues

T. Maroilley<sup>D 1</sup>, M. Berri<sup>2</sup>, G. Lemonnier<sup>1</sup>, D. Esquerré<sup>3</sup>, C. Chevaleyre<sup>2</sup>, S. Mélo<sup>2</sup>,

F. Meurens<sup>2,8</sup>, J. L. Coville<sup>1</sup>, J. J. Leplat<sup>1,4</sup>, A. Rau<sup>1</sup>, B. Bed'hom<sup>1</sup>, S. Vincent-Naulleau<sup>1,4</sup>,

M. J. Mercat<sup>5</sup>, Y. Billon<sup>6</sup>, P. Lepage<sup>7</sup>, C. Rogel-Gaillard<sup>1</sup> & J. Estellé<sup>D 1</sup>

The epithelium of the intestinal mucosa and the gut-associated lymphoid tissues (GALT) constitute an essential physical and immunological barrier against pathogens. In order to study the specificities of the GALT transcriptome in pigs, we compared the transcriptome profiles of jejunal and ileal Peyer's patches (PPs), mesenteric lymph nodes (MLNs) and peripheral blood (PB) of four male piglets by RNA-Seq. We identified 1,103 differentially expressed (DE) genes between ileal PPs (IPPs) and jejunal PPs (JPPs), and six times more DE genes between PPs and MLNs. The master regulator genes *FOXP3*, *GATA3*, *STAT4*, *TBX21* and *RORC* were less expressed in IPPs compared to JPPs, whereas the transcription factor *BCL6* was found more expressed in IPPs. In comparison between IPPs and JPPs, our analyses revealed predominant differential expression related to the differentiation of T cells into Th1, Th2, Th17 and iTreg in JPPs. Our results were consistent with previous reports regarding a higher T/B cells ratio in JPPs compared to IPPs. We found antisense transcription for respectively 24%, 22% and 14% of the transcripts detected in MLNs, PPs and PB, and significant positive correlations between PB and GALT transcriptomes. Allele-specific expression analyses revealed both shared and tissue-specific cis-genetic control of gene expression.

The digestive tract provides nutrients and energy from food to the organism, and it is also a major potential entry site for pathogens. The epithelium of the intestinal mucosa, together with its associated immune system, are in direct contact with both symbiotic and pathogenic microorganisms and thus constitute an essential immunological and physical barrier<fnref n="1" />. The specialized intestine immune system comprises the gut-associated lymphoid tissues (GALT) that include the Peyer's patches (PPs), the mesenteric lymph nodes (MLNs) and isolated lymphoid follicles. PPs and MLNs both act as secondary lymphoid organs specifically associated with the gut tissues². MLNs are induction sites of immune responses, as they initiate peptide presentation to lymphocytes by antigen presenting cells arriving from the intestine². The PPs are distributed along the ileum and jejunum segments and they are thus classified as either ileal PPs (IPPs) or jejunal PPs (JPPs). The PPs are covered with the follicle-associated epithelium that harbors microfold cells (M cells). The sub-epithelial dome region contains abundant professional antigen presenting cells including dendritic cells (DCs) that take up antigens retrieved by M cells from the lumen of the digestive tract. The DCs present peptides to the T and B cells located in the PPs and trigger the initiation of

<footnote n="1">¹GABI, INRA, AgroParisTech, Université Paris-Saclay, 78350, Jouy-en-Josas, France. ²ISP, INRA, Université de Tours, 37380, Nouzilly, France. ³GenPhySE, INRA, INPT, ENVT, Université de Toulouse, 31326, Castenet-Tolosan, France. ⁴LREG, IRCM, DRF, CEA, Université Paris-Saclay, 78350, Jouy-en-Josas, France. ⁵BIOPORC and IFIP-Institut du porc, La Motte au Vicomte, BP 35104, 35651, Le Rheu, France. ⁶GENESI, INRA, 17700, Surgères, France. ⁷MICALIS Institute, INRA, AgroParisTech, Université Paris-Saclay, 78350, Jouy-en-Josas, France. ⁸Present address: BIOEPAR, INRA, Oniris, La Chantrerie, 44307, Nantes, France. Correspondence and requests for materials should be addressed to C.R.-G. (email: claire.rogel-gaillard@inra.fr)</footnote>
antigen-specific immune responses and adaptive immunity. The PPs are also the primary inductive sites where most IgA immune responses are initiated by the class-switching of B cells from IgM to IgA antibodies³.

In humans, JPPs and IPPs are multiple isolated follicles, but in other species such as ruminants and pigs, JPPs are isolated follicles while IPPs are instead continuous structures along the ileum⁴. The overall anatomical differences between JPPs and IPPs in pigs suggest functional differences and their specific roles still need to be clarified⁵,⁶. Indeed, several studies have described histological and functional differences between the two porcine PPs. For instance, the cellular composition of the GALT has been described as follows⁷: B cells represent 90% and 45% of IPP and JPP cells, respectively, and 35% of MLN cells. Other studies have explored the functional and transcriptomic disparities for limited gene sets. For instance, it has been reported that IPPs would not be necessary for the development of the B cell systemic pool⁸ or specific IgA production⁹. In contrast, in a targeted gene expression approach, Levast *et al.*¹⁰ showed differential expression of 29 immune genes between IPPs and JPPs and observed a more diverse IgA repertoire in IPPs. Similarly, Gourbeyre *et al.*¹¹ focused on pattern recognition receptors involved in innate immunity and observed differential expression for *TLR2*, *6*, *7*, *9*, *10* and *NLRP3* between IPPs and JPPs.

The existing published results strongly suggest shared but subtle differences between the functionalities of IPPs and JPPs. To further explore this question, we are presenting in this work a study focusing on transcriptome profiles in pigs with no clinical signs of infection, by stranded RNA-sequencing of IPPs, JPPs, MLNs and peripheral blood as a complementary non-GALT immune tissue. We report the differential gene expression between PPs and MLNs, and more specifically between IPPs and JPPs, together with the analysis of antisense transcription and allele-specific expression (ASE) in the four tissues.

## Materials and Methods

**Animals and sample collection.** All animals were Large White pigs bred in the INRA experimental farm at Le Magneraud (GENESI, UE 1372, France). All animals were weaned at 28 days of age and fed *ad libitum*. They belong to the previously reported SUS\_FLORA cohort¹² for which all animals were sampled for blood at 60 days of age. Among this cohort, 36 pigs were slaughtered at 70 days of age to sample IPPs, JPPs and MLNs. A subset of four uncastrated males was selected for the transcriptome study. The remaining group of 32 pigs that comprised 16 males and 16 females was used as a validation group for quantitative real-time PCR (qRT-PCR) studies. All animal experiments were carried out in accordance with European Guidelines for the Care and Use of Animals for Research Purposes. The animal protocol was approved by the local ethics committee in Poitou Charentes and assigned the approval number CE2013-2.

Peripheral blood (PB) was sampled in the jugular vein using PAXgene Blood RNA tubes (PreAnalytiX, Qiagen, Germany). Pieces of MLNs, IPPs and JPPs were collected, snap frozen in liquid nitrogen, and stored at -80 °C until use.

**RNA extraction and sequencing.** Tissue RNA extractions were performed as previously reported¹¹,¹². In brief, 60 mg of tissue were lysed in 1 ml of Trizol (Invitrogen, Cergy Pontoise, France) with ceramic beads (Bertin technologies, St Quentin en Yvelines, France), and total RNA was purified using RNeasy Mini Kit (Qiagen, Courtaboeuf, France) according to the manufacturer's recommendations. Residual genomic DNA was removed using DNase digestion with RNase-free DNase I Amplification Grade (Invitrogen, Cergy Pontoise, France) following the recommended protocol. For blood samples, no globin depletion procedure was applied, and RNA purification was performed as reported in Maroilley *et al*.¹³. One µg of total RNA was reverse transcribed for 90 min at 37 °C in a 20 µl volume containing 0.25 mM dNTP of OligodT, 25 U of MuMLV reverse transcriptase in 4 µl 5X MuMLV buffer (Eurogentec, Liège, Belgium). After heat-inactivation at 93 °C for 5 min, generated cDNA was stored at -80 °C until use. Finally, reverse-transcribed RNA samples of MLNs, IPPs, JPPs and PB from the four pigs chosen for transcriptome analysis were sequenced using three lanes in an Illumina HiSeq 2000 with the Illumina TruSeq Stranded mRNA Library Prep Kit. On average, 47 million reads (100 bp paired-end) were obtained per sample. The raw reads are available at NCBI's SRA repository (Bioproject PRLNA416857; accessions SAMN07966348 to SAMN07966363).

**Differential expression analyses of RNA-Seq data.** After quality control using the FastQC tool, sequencing reads were filtered with the Trimmomatic tool (v.0.32)¹⁴ by trimming leading and trailing bases with phred qualities less than 3 and dropping reads shorter than 36 bases long and those with average phred qualities per base less than 15. Reads were subsequently aligned on the reference genome Sscrofa10.2 using TopHat2 (v2.0.14), with the default parameters and the gene annotation available for *Sus scrofa* in the Ensembl 85 release database¹⁵. Based on the read mapping with TopHat2, gene expression was quantified by obtaining read counts with the HTSeq-count software (v.0.6.1p1)¹⁶ with a default parameter that discards all reads that map to multiple locations. Expression levels of antisense transcripts were also analyzed by using the strand read information. We quantified antisense transcription according to the number of reads mapping to the corresponding gene reference sequence on the opposite strand. Therefore, we did not use the annotation of antisense transcripts available in Ensembl 89 but only the annotation of reference genes with their genomic position. Read counts were further analyzed using the edgeR R/Bioconductor package (v.3.12.1)¹⁷. We retained genes as expressed in a tissue when the count per million (CPM) sense reads was greater than one for at least two animals. Similarly, a gene was found to have an antisense transcription if the antisense read CPM was greater than one for at least two animals. Read counts were then normalized according to the total number of reads of each sample using the Trimmed Mean of M-values normalization method (TMM)¹⁸ implemented in edgeR (Supplementary Tables S1 and S2). The sense and antisense expression data structures were explored with multi-dimensional scaling (MDS) plots. The correlation between sense and antisense transcription across tissues was estimated by Spearman correlation using the normalized read counts.
The differentially expressed (DE) genes between GALT tissues were detected by fitting a negative binomial generalized linear model (GLM) in edgeR. In order to take into account intra-individual variability, the model included covariates for both tissue and individual. A likelihood ratio test was performed to identify DE genes among each pair of tissues, and *P*-values were corrected for multiple testing using the Benjamini-Hochberg control of the False Discovery Rate (FDR < 0.05). Finally, smear plots of log-fold expression changes versus log-concentration were produced using the “plotSmear” function in the edgeR package.

**Joint genotyping and inference of allelic specific expression.** For ASE analysis, we realigned reads along the reference genome Sscrofa10.2 using the STAR 2-pass protocol (v2.4.0i)¹⁹. After marking duplicates, we split the reads, realigned them on indels and recalibrated the sequence data following the GATK (v3.7) Best Practices for calling variants in RNA-Seq data²⁰, and as described in Maroilley et al.¹³. SNPs were called with the HaplotypeCaller of the GATK tool (v3.7) using all samples and setting a base quality score higher than 10 and a coverage greater than 10 reads. We then removed SNPs with a low minor allele count (less than 3 reads that harbor the alternative allele) using VCFtools (v0.1.12a)²¹. In order to limit alignment biases of reads due to the presence of alternative SNP alleles, we created a masked genome reference sequence with an N at each SNP position using the Bedtools maskbast tool (v2-2.24.0) and re-aligned reads on the masked genome with the same STAR 2-pass protocol. Using the dataset from this second alignment, counts of reads covering each allele at selected SNPs were obtained using ASEReadCounter of GATK v3.7 with parameters ensuring an adequate coverage and quality (read depth greater than 10, mapping quality greater than 10, base quality greater than 2). We then adapted ASEReadCounter output for the format required by QuASAR v0.1.²² and used the associated statistical method to genotype individuals from multiple RNA-Seq samples (fitAseNullMulti) and determine for each animal the probability of each possible genotype according to the alleles found in each tissue (with a minimum read coverage equal to five). Finally, we performed binomial tests on allelic read counts to detect potential allelic imbalance for the SNPs found to be heterozygous in each pig. The calculated *P*-values were corrected for multiple testing using the Benjamini-Hochberg control of the False Discovery Rate (FDR < 0.05).

**Functional enrichment analysis.** The Sus scrofa Ensembl gene IDs (release 89 corresponding to Sscrofa10.2 assembly) annotated from DE and ASE analyses were submitted to Biomart in order to identify lists of human orthologous genes for subsequent functional enrichment analyses. The annotation was further updated using the Ensembl 91 release (Sscrofa11.1 assembly). For 869 genes still lacking GeneName HUGO Genome Nomenclature Committee (HGNC) annotation, we aligned the protein sequences corresponding to the pig gene Ensembl ID against the human RefSeq protein database using BLASTP²³. The GOrilla tool²⁴ was used for enrichment analyses of expressed genes in each tissue. For DE genes, we used the Ingenuity Pathway Analysis tool (IPA; http://www.ingenuity.com) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway mapper²⁵. SNPs and ASE SNPs were annotated with the Variant Effect Predictor tool (Ensembl).

**Quantitative real-time PCR assays.** In order to validate the differential expression of a set of immune response mediator genes, qRT-PCRs were carried out on IPPs, JPPs and MLNs from a group of 70-day-old pigs. Three genes, *β-2-microglobulin* (B2M), *hydroxymethylbilane synthase* (HMBS) and *glyceraldehyde-3-phosphate dehydrogenase* (GAPDH), were used as reference genes because of their stable expression between samples as evaluated using the geNorm algorithm²⁶. All gene primers (Supplementary Table S3) were designed using the Clone Manager 9 software (Scientific & Educational Software, Cary, NC, USA). For cDNA synthesis, one µg of total RNA was reverse transcribed for 90 min at 37°C in a 20 µl volume containing 0.25 mM dNTP of OligoT, 25 U of MuMLV reverse transcriptase in 4 µl 5X MuMLV buffer (Eurogentec, Liège, Belgium). After heat-inactivation at 93°C for 5 min, generated cDNA was stored at −80°C until use. qRT-PCRs were performed on a Bio-Rad Chromo 4 apparatus (Bio-Rad, Hercules, CA, USA) using the Mesa Green qPCR MasterMix Plus for SYBR assay (Eurogentec, Liège, Belgium) as previously described²⁷⁻²⁹. All qRT-PCRs displayed efficiency between 90 and 110%, and expression data were given as relative values using modified delta delta Ct method³⁰ after Genex macroanalysis (Bio-Rad, Hercules, CA, USA)²⁶.

**Pyrosequencing validation of selected ASE signals.** A subset of ASE SNPs annotated to six genes were selected for validation studies by pyrosequencing with a PyroMark Q24 instrument (Qiagen). Supplementary Table S4 shows primer sequences used for PCR amplification and subsequent sequencing. For each cDNA sample, PCRs were performed using the HotStart Taq polymerase following the conditions specified in Supplementary Table S4. The allele quantification was subsequently performed in the PyroMark Q24 system according to manufacturer's instructions.

## Results

**Global transcription in GALT and peripheral blood.** The transcriptome analyses were carried out by RNA-Seq on the same four male piglets for blood samples (60-day-old), and for MLNs, IPPs and JPPs (70-day-old). The global expression data structure was explored with a multi-dimensional scaling plot (MDS plot, Fig. 1a) based on the TMM normalized gene expressions for each sample. In this two-dimensional plot, samples are positioned according to the statistical distance of their expression profiles, highlighting here that samples are grouped rather by tissue than by individual. Similar results were obtained from normalized antisense read counts (Fig. 1b). For both sense and antisense transcription levels, we observed a clear separation between MLNs, PB and PPs, whereas IPPs and JPPs grouped together, thus revealing strong similarities between PP types at this scale of comparison.

Globally, we found over 13,000 expressed genes (according to the pig genome annotation) in IPPs, JPPs and MLNs. In PB, only ~9,200 genes were detected (Table 1), which is likely an under-estimation since no globin depletion procedure was applied on blood samples, thus limiting the detection of weakly expressed genes due
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Summary plots of global gene expression analyses performed by RNA-Seq on MLN, IPP, JPP and PB. Graphs show MDS plots for normalized gene expression levels (a) and normalized antisense transcription levels (b).</figcaption></figure>

<table><caption>**Table 1.** Number of genes with sense and/or antisense transcription in the four tissues.</caption><thead><tr><th></th><th>JPP</th><th>IPP</th><th>MLN</th><th>PB</th></tr></thead><tbody><tr><td>Expressed annotated genes</td><td>13,767</td><td>13,646</td><td>13,463</td><td>9,201</td></tr><tr><td>Number of genes that represented 50% of sense reads</td><td>648</td><td>619</td><td>719</td><td>2</td></tr><tr><td>Genes with antisense expression</td><td>3,155</td><td>3,138</td><td>3,616</td><td>1,345</td></tr><tr><td>Number of genes that represented 50% of antisense reads</td><td>38</td><td>34</td><td>43</td><td>1</td></tr><tr><td>Genes with both sense and antisense transcription</td><td>2,841 (21%)</td><td>2,836 (21%)</td><td>3,242 (24%)</td><td>1,059 (12%)</td></tr><tr><td>Genes harboring only antisense transcription</td><td>314</td><td>302</td><td>374</td><td>286</td></tr></tbody></table>

to an over-representation of globin genes. Indeed, HTSeq-count counted over 18 million reads for the HBB (ENSSSCG00000014725) gene and almost 12 million reads for the HBA (ENSSSCG00000007978) gene, corresponding to 43% and 27% of the reads, respectively. Additionally, an HBB-like gene (ENSSSCG00000014727) was represented by four million reads (4% of the reads).

We quantified the transcription level of elements transcribed from the opposite strand of annotated genes in Ensembl (release 85) and found over 3,100 genes with measurable antisense transcription levels in PPs and MLNs, and less than half of that number in PB (Table 1). Globally, 4,138 loci annotated as genes on one strand in Ensembl were detected as transcribed on the other strand in at least one tissue.

As presented in Table 1, the GALT and PB showed an antisense transcription for over 21% and 12% of expressed genes, respectively. The five genes with the largest number of antisense overlapping reads were the same in all tissues: GYS1 (ENSSSCG00000003154) with 48% of antisense reads in blood, NDUFB10 (ENSSSCG000000030328), RNF167 (ENSSSCG00000017906), ND6 (ENSSSCG00000018092) and FIB (ENSSSCG00000012969). Of note, NDUFB10 and ND6 are both members of the NADH: ubiquinone oxidoreductase enzymatic complex. Besides, we detected an antisense transcription level for 598 non-expressed genes in our dataset (see full list in Supplementary Table S5). Interestingly, GALT shared 230 common annotated genes harboring only antisense transcription, and 101 antisense expressed events were shared by the four tissues.

We mapped the sense and antisense reads along the genome (see Supplementary Fig. S1 in Supplementary Information file). We observed that they are not homogeneously distributed, and that antisense transcription was preferentially localized on genomic regions showing sense transcription. Our results highlighted similar genomic regions across the four tissues that appear to be actively transcribed on both DNA strands.

We estimated correlations between sense and antisense gene transcription within and between tissues for the 14,596 genes with a CPM greater than 1 in at least one tissue and for the 4,138 genes with an antisense transcription level (CPM greater than 1 with at least 10 reads). The correlation of sense and antisense gene expression levels was found to be higher between IPPs and JPPs than between PPs and MLNs (Table 2). The correlation was lower between GALT and PB, and MLNs appeared to be more correlated with PB. No significant correlations were found between the antisense and sense read counts for a same gene, whatever the tissue.

**Differential gene expression analysis between GALT.** We performed a differential expression analysis within the GALT to explore differences in the transcriptome between PPs and MLNs as well as between IPPs and JPPs, based on the 14,502 genes found to be expressed (sense transcription with CPM > 1) in at least one of the three tissues (Fig. 2a). JPPs, IPPs and MLNs shared the same five most highly expressed genes: cytochrome
<table><caption>**Table 2.** Spearman correlation coefficients of gene expression (below diagonal) and antisense transcription (above diagonal) between tissues; Spearman correlation coefficients between sense and antisense gene transcription in each tissue (diagonal).</caption><thead><tr><th colspan="5">Antisense transcription</th></tr><tr><th>Gene expression</th><th>MLN</th><th>IPP</th><th>JPP</th><th>PB</th></tr></thead><tbody><tr><th>MLN</th><td>-0.04</td><td>0.72</td><td>0.72</td><td>0.71</td></tr><tr><th>IPP</th><td>0.83</td><td>-0.06</td><td>0.84</td><td>0.56</td></tr><tr><th>JPP</th><td>0.84</td><td>0.97</td><td>-0.06</td><td>0.58</td></tr><tr><th>PB</th><td>0.79</td><td>0.68</td><td>0.69</td><td>-0.1</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Summary plots of gene expression analyses within GALT. Graphs show an MDS plot for normalized gene expression levels for GALT (a), a Venn diagram for the comparison of expressed genes in GALT (b), and smear plots of differential gene expression analyses by comparing IPP and JPP (c), MLN and IPP (d) and JPP and MLN (e).</figcaption></figure>

c oxidase subunit 1 (COX-1), cytochrome c oxidase subunit 2 (COX-2), cytochrome c oxidase subunit 3 (COX-3), ATP synthase F0 subunit 6 (ATP6) and eukaryotic translation elongation factor 1 alpha 1 (EEF1A1). The COX-1, COX-2 and COX-3 mitochondrial genes, coding for cytochrome c oxidases, were represented by 8% and 5% of the sequencing reads in PPs and MLNs, respectively. ATP6 is also a mitochondrial gene coding for an ATP synthase, and EEF1A1 codes for a translation elongation factor.

In a first exploratory approach, an MDS plot for the GALT (Fig. 2b) showed that IPP and JPP samples were separated by intestinal segment, revealing differences in gene expression between IPPs and JPPs. We subsequently tested whether gene expression was significantly different between the GALT tissues using the edgeR package (Table 3). The complete lists of DE genes is available in Supplementary Tables S6 to S8. There were six times more differences between MLNs and PPs than between IPPs and JPPs. However, since fewer DE genes were detected, JPPs seemed to be more similar to MLNs than IPPs. Smear plots showed that the bulk of genes was centered around a log-fold change of zero indicating that any composition bias between libraries was successfully removed (Fig. 2c-e). We observed that the scale of fold changes between MLNs and PPs was larger than between IPPs and JPPs.

DE genes between the PPs and MLNs were involved in cell surface receptor signaling pathway (GO:0007166) and signal transduction (GO:0007165). They were also involved in immune response (GO:0006955). They were implicated in cell adhesion (GO:0007155) and regulation of cell adhesion (GO:0030155), regulation of
<table><caption>**Table 3.** Differential gene expression among the GALT samples.</caption><thead><tr><th rowspan="2"></th><th colspan="3">FDR &lt; 0.05</th><th colspan="3">FDR &lt; 0.05 and |log<sub>2</sub>FC| &gt; = 1</th></tr><tr><th>DE</th><th>Over-expressed</th><th>Under-expressed</th><th>DE</th><th>Over-expressed</th><th>Under-expressed</th></tr></thead><tbody><tr><td>IPP VS MLN</td><td>6,852</td><td>3,405</td><td>3,447</td><td>3,922</td><td>2,143</td><td>1,779</td></tr><tr><td>JPP VS MLN</td><td>6,110</td><td>3,058</td><td>3,052</td><td>3,255</td><td>2,008</td><td>1,247</td></tr><tr><td>IPP VS JPP</td><td>1,103</td><td>444</td><td>659</td><td>549</td><td>217</td><td>332</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** The most significantly enriched biological processes associated with differentially expressed genes between JPP and MLN (red) and IPP and MLN (in blue), performed with the GOrilla tool on human ortholog GeneNames. The bars represent the number of genes found DE involved in a given biological process, without taking into account whether the gene is under- or over-expressed in the comparison. The full lists of under- and over-expressed DE genes are available in Supplementary Tables S6–8.</figcaption></figure>

locomotion (GO:0040012) and regulation of cell motility (GO:2000145). They also affected metabolic processes such as lipid (GO:0006629) or organic metabolism processes (GO:0006082) (Fig. 3). The DE genes between IPPs and JPPs were found to be involved in 60 pathways referenced in the KEGG database²⁵, where each of these pathways contained more than ten DE genes (Supplementary Table S9). The pathway that included the most DE genes (106) was “Metabolic pathway”, the second (44 DE genes) was “Pathways in cancer”, and the third (40 DE genes) was “Cytokine-cytokine receptor interaction”. The “Th1 and Th2 cell differentiation” pathway was also enriched and included 19 DE genes, of which 18 and one were under- and over-expressed in IPPs compared to JPPs, respectively (Fig. 4).

By comparing these results to our previous study that targeted transcriptome analysis along the small intestine (duodenum, jejunum, ileum)¹², we found an overlap of 431 genes between the set of DE genes previously identified between jejunum and ileum and the 1,105 DE genes between IPPs and JPPs. The fold change ranges were quite similar in both studies (respectively 1.7 versus 1.9 log₂FC), and 97% of the shared genes showed a similar trend of differential expression, with only 14 genes presenting contradictory differential expression patterns between ileum/jejunum and IPPs/JPPs. The 431 shared DE genes were involved in pathways referred to as “LPS/IL-1 Mediated Inhibition of RXR Function”, “FXR/RXR Activation”, and “LXR/RXR Activation”. The top ten DE genes in term of absolute log-fold change most DE genes were assembled by IPA in various networks related to lipid metabolism (Supplementary Table S10), as reported in Mach *et al.*¹². In addition, over 60% of the DE genes (672 genes) between both PP types were not found DE in the intestinal segments in Mach *et al.*'s results¹². In enrichment analyses, DE genes specific to lymphoid tissues were found to be involved in several immunity-related canonical pathways: “Th1 and Th2 Activation Pathway”, “Primary Immunodeficiency Signaling”, “Th1 Pathway”, “iCOS-iCOSL Signaling in T Helper Cells”, and “T Cell Receptor Signaling”. The most enriched functions were “Hematological System Development and Function”, “Tissue Morphology”, and “Lymphoid Tissue Structure and Development”. The top ten DE genes between IPPs and JPPs assembled in IPA networks were related to biological functions linked to lymphoid tissue structure, cellular development and cell-mediated immune response.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Th1 and Th2 cell differentiation pathways found enriched in DE genes between IPP and JPP by KEGG²⁵. Genes over-expressed in IPP are shown in yellow and those over-expressed in JPP in red. The log-fold changes of each DE gene have been added to the figure.</figcaption></figure>

**Allelic specific expression analyses in porcine GALT and blood transcriptomes.** Beyond the observed differential expression patterns between the GALT, we further explored the *cis*-genetic control of gene expression in each tissue using an ASE analysis. Overall, 432,198 SNPs were called across the 16 samples (MLN, IPP, JPP and PB from four individuals), after filtering. Each animal was genotyped according to the RNA-Seq data for the four tissues with the QuASAR R package²², and 157,178 SNPs were found to be heterozygous in at least one animal. While 46% of the SNPs were shared by at least two pigs and 5,600 SNPs were heterozygous for the four pigs, each animal still harbored a large number of specific heterozygous SNPs (17,224 to 25,160).

For each tissue, we tested SNPs heterozygous in at least one animal (Table 4) for possible ASE. In total, 157,178 heterozygous SNPs overlapping 12,755 genes were tested for ASE in at least one tissue (Supplementary Table S11). Results showed 16,797 significant ASE SNPs, 52% of which were found significant for every heterozygous animal (Table 4). Overall, 37% of the 12,755 genes tested had at least one SNP with ASE in at least one tissue (Table 4), and we found 2,509 genes affected by more than one heterozygous SNP ASE.

The proportion of genes showing significant ASE in at least one individual was close to 25% for GALT and 17% for blood. Of the genes showing significant ASE, 16% displayed ASE in the four tissues and were found to be involved in immune related functions such as antigen processing and presentation, T cell receptor signaling pathway of type I interferon signaling pathway by enrichment analyses. In PB, 141 genes were specifically under a *cis*-genetic control while they were also expressed in other tissues. In MLNs, 669 genes specifically displayed an ASE. By comparing MLNs and PPs, 113 genes expressed in GALT displayed specific ASE in PPs. 500 genes expressed in JPPs and 572 genes expressed in IPPs had ASE while expressed in the other tissues with no allelic imbalance (Supplementary Table S12). Nevertheless, the ASE genes specific to each of these tissues were not enriched in any functional category.

**Validation of differentially expressed genes and ASE signals.** qRT-PCRs were carried out for 24 genes (Supplementary Table S3) using tissues collected from the group of 32 pigs of 70 days of age belonging to the same cohort as the four animals included in the transcriptome analysis. The results confirmed that the *CCR10* and *CCR9* genes were under-expressed in IPPs versus JPPs (Fig. 5 and Supplementary Table S13). Furthermore, we showed that *CCR10* was less expressed in IPPs compared to MLNs, but more highly expressed in JPPs compared to MLNs and IPPs. We also validated that *SLA-7* was not differentially expressed between GALT samples. qRT-PCRs confirmed that *CCL28* and *IL-15* were more expressed in PPs than in MLNs, and that *SLA-6*, *SLA-8*, *MIC-2*, *MR1*, *IL-10*, *CCR9*, *MADCAM1*, *NFKB1*, *SMAD2* and *SMAD3* were down-regulated in PPs in comparison to MLNs.

As a first validation step of the ASE results, we compared the current ASE analysis on blood with a previous study that we conducted on 38 60-day-old Large-White pigs on peripheral blood¹³. Here, we revealed 1,005 ASE-SNPs (FDR < 0.05, at least 10 reads) found in at least two heterozygous pigs, of which 39% of these SNPs
<table><caption>**Table 4.** Summary of ASE results obtained from RNA-Seq datasets for the four tissues analyzed in this study (MLN, IPP, JPP and PB). ¹Genes with a heterozygous SNP (<5 kb distance from gene according to Ensembl's Variant Effect Predictor). ²For each animal heterozygote at this position, the SNP was detected with a significant ASE.</caption><thead><tr><th>TISSUE</th><th>HETEROZYGOUS SNPS</th><th>GENES<sup>1</sup></th><th>ASE SNPs (% tested)</th><th>Consistent ASE SNPs<sup>2</sup> (% ASE SNPs)</th><th>ASE Genes (% tested)</th><th>ASE Genes > 1 ASE SNP (% ASE GENES)</th></tr></thead><tbody><tr><td>MLN</td><td>129,122</td><td>11,667</td><td>8,295 (6.4%)</td><td>4,046 (49%)</td><td>2,945 (25%)</td><td>1,543 (52%)</td></tr><tr><td>IPP</td><td>111,378</td><td>11,576</td><td>8,774 (7.9%)</td><td>4,265 (49%)</td><td>2,956 (26%)</td><td>1,524 (52%)</td></tr><tr><td>JPP</td><td>113,475</td><td>11,570</td><td>8,583 (7.5%)</td><td>4,161 (48%)</td><td>2,877 (25%)</td><td>1,420 (49%)</td></tr><tr><td>PB</td><td>40,584</td><td>7,410</td><td>3,339 (8.2%)</td><td>1,787 (54%)</td><td>1,254 (17%)</td><td>614 (49%)</td></tr><tr><td>Total</td><td>157,178</td><td>12,755</td><td>16,787 (11%)</td><td>8,746 (52%)</td><td>4,769 (37%)</td><td>2,509 (53%)</td></tr></tbody></table>

(395) were already described as ASE-SNPs in our precedent study. Among them, 332 SNPs are described in dbSNP and include 55% missense variants, 43% synonymous variants and 1% stop gained variants. These 332 SNPs affected 93 genes.

In a second step, we used the pyrosequencing allele quantification approach to evaluate the consistency of the ASE effects detected in six genes and representing a range of allele frequencies and tissue specificities: *lysozyme* (LYZ), *immunoglobulin lambda variable 10 (IGVL10)*, *complement factor D (CFD)*, *C-X-C motif chemokine ligand 13 (CXCL13)*, 2,−5,−oligoadenylate synthetase 2 (OAS2) and *germinal center-associated, signaling and motility (GCSAM)*. Out of these six genes, we validated the RNA-Seq results in three of the cases representing a tissue-specific ASE (LYZ, with ASE confirmed only in MLNs and PB) and two strong ASE effects (IGLV10 and OAS2). However, ASE effects of smaller magnitude could not be confirmed for *CFD*, *CXCL13* and *GCSAM* (Supplementary Table S14).

Discussion

Improving individual health in pig production systems requires a good understanding of immune and intestine pig functions to foster new disease preventive approaches. The pig species is also considered to be an excellent model for biomedical research, particularly for intestinal diseases due to anatomical and physiological similarities with the human digestive system³¹. However, swine PPs are an exception to this rule because, as in ruminants and rabbits, they have the particularity of forming a continuous structure in the ileum, while in humans they constitute independent follicles in all intestinal segments as in the pig jejunum³¹,³⁰,³².

We report here a comparative analysis of the transcriptomic profiles of GALT in pigs including IPPs, JPPs and MLNs, and PB as non-GALT reference immune tissue. The results showed extensive expression differences across GALT that were likely to reflect tissue specificities. The mapping of the cis-genetic control of gene expression with ASE analyses uncovered both shared and specific ASE signals in each tissue. Overall, our results propose a global picture of the transcriptomic control and profiles of GALT tissues in 70-day-old pigs in good health as a reference study to pave the way for a better understanding of the molecular mechanisms of PP functions in pigs.

The porcine IPPs and JPPs have transcriptomic specificities linked with immunity. The anatomical differences between IPPs and JPPs in pigs have been previously described⁴,³³, and several studies have reported the transcriptome profiles of JPPs³⁴,³⁵. However, to our knowledge, comparative analyses of gene expression between JPPs and IPPs in pigs have only been investigated in limited sets of target immune response-related-genes¹⁰,¹¹. Thus, our study is the first whole-transcriptome approach to this question, revealing 1,103 DE genes between IPPs and JPPs in pigs. Since, to the best of our knowledge, there are not equivalent reports on ruminants like sheep and cattle, our results could also be relevant for these species. The differential gene expression profiles between the two PP types could be partly related to their anatomical position along the small intestine. We thus compared the IPP and JPP transcriptomes to those of the ileum and jejunum sections previously reported¹². We confirmed that genes found DE between IPPs and JPPs represented specific differences between the two PP types and were not simply reflecting the gene expression profiles of their corresponding intestinal segment.

We analyzed the gene sets found DE between IPPs and JPPs by performing IPA-based enrichment studies in order to explore which pathways were affected. We observed that canonical pathways related to T cell activity were down-regulated in IPPs while the B cell receptor-signaling pathway was up-regulated in IPPs compared to JPPs. These differences were consistent with the variations in lymphocyte composition described in pigs, since IPPs were reported to contain more B cells than JPPs (90% and 45% of B lymphocytes, respectively), and fewer T cells⁷. Several studies have described that the relative number of cells entering the IPPs was lower than that entering the JPPs in pigs³⁶,³⁷. Accordingly, we found that the receptor *CCR7* and its ligand *CCL21* that mediate the organization of the T zone³⁸ were less expressed in IPPs. The *CCR10* gene required for optimal T cell dependent IgA-antibody secreting cells accumulation in mice³⁹ was also found to be less expressed in the IPPs (logFC = −1.3). Furthermore, the receptor *CTLA4*, involved in the regulation of activated T-lymphocytes⁴⁰, was found to be less expressed in the IPPs, and *CCR9*, which promotes the homing of IgA-antibody secretion T cells towards the intestine⁴¹, was found more expressed in JPPs (logFC = 0.8) than in IPPs as described previously¹⁰. Thus, all together our results are consistent with a higher proportion of T cells in the JPPs compared to the IPPs.

Our results highlighted gene expression differences between IPPs and JPPs related to the differentiation and fate of T cells (see Supplementary Table S9). The master regulators of CD4+ T cell commitment and plasticity have been extensively studied⁴². A slight but significant decrease in the master regulators known to be involved in naïve CD4+ T cell differentiation was observed in IPPs compared to JPPs (FC between −1.8 and −2.6 and logFC between −0.87 and −1.42 see Fig. 4). The transcription factors found less expressed in IPPs were *FOXP3* (differentiation of induced
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Validation of the differential expression of a subset of genes between IPP, JPP and MLN by qRT-PCR.</figcaption></figure>

T regulator cells or iTreg), GATA3 (Th2 differentiation), STAT4 and TBX21 that encodes T-bet (Th1 differentiation), as well as RORC that encodes Rorγt (Th17 differentiation). Conversely, the transcription factor BCL6 was found over-expressed in IPPs compared to JPPs (FC = 1.8 and logFC = 1.83). BCL6 is involved in the differentiation of the follicular helper T (Tfh) cells that have the important function of providing B-cell help for the induction of antigen-specific antibody production⁴³. The Tfh cells contribute to the formation of germinal centers, which are dynamic microenvironments that provide a unique niche for B-cell affinity maturation to occur (review in Ramiscal and Vinuesa)⁴⁴. JCS, known to be expressed on surface of Tfh cells, was however found less expressed in IPPs (FC = -2.6 and logFC = -1.4) than in JPPs as well as IL-21 (FC = -2.8 and logFC = -1.5), known to be produced by Tfh cells⁴³. It has been reported in mice that a fraction of Foxp3+ T cells can become Tfh cells in the PPs after down-regulating FOXP3 expression⁴⁵. This finding is in agreement with our results that revealed opposite differential expression between FOXP3 and BCL6 in IPPs and JPPs. Both JPPs and IPPs are main sites for priming B cells and educating the adaptive immune system while facing the gut microbiota. Tfh cells seem to play a crucial role in sensing microbiota-derived extracellular ATP and shaping commensal microbiota composition⁴⁶. Our results suggest that IPPs are more abundant in Tfh and B cells than JPPs, and this could be linked to the densification of the microbiota along the digestive tract. Even if we cannot exclude that the differential expression found between IPPs and JPPs could be solely due to distinct proportions of B and T cells, our results also suggested that T cells could have different fates in IPPs and JPPs.

Regarding the coherence of our results with previous reports, we observed that *CCR7*-CCL21, which mediate the organization of PPs into follicles, were down-regulated in the IPPs. Moreover, we found that the *TNFSF11* gene involved in the development of GALT was less expressed in IPPs in comparison with JPPs (logFC < -1). These results are consistent with differences in the development kinetics between JPPs and IPPs as reported in two month old piglets, with a slower development of IPPs compared to JPPs⁴⁷. In addition, Gourbeyre et al.¹¹ highlighted differential expression in JPPs and IPPs for *TLR2*, *TLR6*, *TLR7*, *TLR9*, *TLR10*, and *NLRP3*. We also observed a differential expression of *TLR7* and *TLR8* (IPPs < JPPs) and *TLR9* (IPPs > JPPs).
MLNs and PPs participate in the activation of lymphocytes with the presentation of antigens passing through the intestine⁴⁸. Interestingly, *TNFSF11*, a gene involved in the development of lymph nodes and not of early steps in the PP development, was found under-expressed in PPs in comparison with MLNs. Boeker *et al*.³³ explored the distribution of the major lymphocyte populations in the peripheral blood and lymphoid organs in six female Göttingen minipigs. Our results were consistent with their immune-histological data. Indeed, we observed an under-expression of the chemokines *CXCL12*, *CXCL13*, *CXCR5*, and *CX3CR1* that mediate the attachment of cells to the lymphoid organs in PPs by comparison with MLNs.

Beyond describing the differences in gene expression within GALT, we studied the regulation of gene expression by performing antisense transcript and ASE analyses in MLNs, JPPs, IPPs and PB. We did not intend here to predict new antisense transcripts in intergenic regions, but only to quantify antisense transcription in annotated genes. Similarly to the work of He *et al*.⁴⁹ on five different cell types, we observed that antisense transcription was not homogenous along the genome. Antisense transcription also seemed to be more frequent in regions with sense transcription. This might be due to a density reduction of the nucleosome structure where genes are transcribed, which may facilitate genome accessibility for antisense transcription. Chen *et al*.⁵⁰ studied antisense transcription in liver and muscle of a pig F2 generation. They observed that 26% to 32% of liver transcripts and 46% to 54% of muscle transcripts had an antisense transcription, respectively. In our study, we found 24%, 22% and 14% genes with antisense transcription in MLNs, PPs and PB, respectively. The absence of significant correlations between sense and antisense transcription levels suggest that sense and antisense transcriptions are not co-regulated and are subjected to distinct regulation mechanisms.

Regarding ASE, we found several effects across the four tissues (Table 4) and the enrichment analyses showed an overabundance of immune-related functions associated with the genes shared between the four tissues. Overall, 37% of the tested genes displayed an ASE effect in at least one tissue, indicating the extensive *cis*-regulation present in the immunome. A recent study conducted by Chamberlain *et al*.⁵¹ in cattle showed even more *cis*-regulation, with 89% of genes harboring heterozygous SNPs having at least one SNP with significant ASE. Similarly, Crowley *et al*.⁵² reported that 89% of all genes tested in mouse brain show ASE.

We observed an important specificity in the *cis*-genetic control of gene expression between MLNs and PPs, although intriguingly there were not specific enriched functions associated with these genes. For individual tissues, the proportion of genes showing significant ASE is close to 25% for GALT and 17% for blood. These proportions are in agreement with Chamberlain *et al*.⁵¹ who previously reported percentages of genes showing significant ASE for white blood cells (33%), spleen (12%) and intestinal lymph node (31%) in cattle. The GTEx consortium⁵³⁻⁵⁵ found that between 1.5 and 3.7% of 6,385 tested SNPs showed ASE; here we found globally 11% of tested SNPs showed an ASE, confirming a likely significant proportion of regulatory SNPs in the pig genome¹³.

Since the blood transcriptome is widely used as a source of biomarkers, we have performed global correlation analyses to evaluate whether it recapitulates at least in part the transcriptomes of GALT. The PPs, MLNs and PB are complex multicellular tissues, and their gene expression profiles reflect the transcription levels of multiple cell types with an expected over-representation of genes expressed in the most abundant cell subset. As expected, transcripts from PB were dominated by the major blood proteins (hemoglobin and ferritin) and reduced the number of low expression genes that could be detected. Nevertheless, we observed strong positive correlations between the transcription levels of genes expressed in blood and GALT for both sense and antisense transcripts (Table 2), despite the different time of sampling between blood and GALT. These data support the use of blood in pigs as a surrogate tissue and a potent source of biomarkers for functions taking place in other tissues. In addition, we found considerable overlap of ASE effects between GALT and blood, highlighting a shared *cis*-genetic control for the expression of immunity-related genes.

## References

1. Abreu, M. T. Toll-like receptor signalling in the intestinal epithelium: how bacterial recognition shapes intestinal function. *Nat. Rev. Immunol.* **10**, 131–144 (2010).

2. Buettner, M. & Bode, U. Lymph node dissection—understanding the immunological function of lymph nodes. *Clin. Exp. Immunol.* **169**, 205–212 (2012).

3. Reboldi, A. & Cyster, J. G. Peyer's patches: organizing B-cell responses at the intestinal frontier. *Immunol. Rev.* **271**, 230–245 (2016).

4. Chu, R. M. & Liu, C. H. Morphological and functional comparisons of Peyer's patches in different parts of the swine small intestine. *Vet. Immunol. Immunopathol.* **6**, 391–403 (1984).

5. Rothkötter, H. J., Pabst, R. & Bailey, M. Lymphocyte migration in the intestinal mucosa: entry, transit and emigration of lymphoid cells and the influence of antigen. *Vet. Immunol. Immunopathol.* **72**, 157–165 (1999).

6. Bailey, M., Christoforidou, Z. & Lewis, M. C. The evolutionary basis for differences between the immune systems of man, mouse, pig and ruminants. *Vet. Immunol. Immunopathol.* **152**, 13–19 (2013).

7. Andersen, J. K. et al. Systematic characterization of porcine ileal Peyer's patch, I. apoptosis-sensitive immature B cells are the predominant cell type. *Immunology* **98**, 612–621 (1999).

8. Yamamoto, M. et al. Alternate mucosal immune system: organized Peyer's patches are not required for IgA responses in the gastrointestinal tract. *J. Immunol.* **164**, 5184–5191 (2000).

9. Sinkora, M. et al. Ileal Peyer's patches are not necessary for systemic B cell development and maintenance and do not contribute significantly to the overall B cell pool in swine. *J. Immunol.* **187**, 5150–5161 (2011).

10. Levast, B. et al. Differences in transcriptomic profile and IgA repertoire between jejunal and ileal Peyer's patches. *Dev. Comp. Immunol.* **34**, 102–106 (2010).

11. Gourbeyre, P. et al. Pattern recognition receptors in the gut: analysis of their expression along the intestinal tract and the crypt/villus axis. *Physiol Rep* **3** (2015).

12. Mach, N. et al. Extensive expression differences along porcine small intestine evidenced by transcriptome sequencing. *PLoS ONE* **9**, e88515 (2014).

13. Maroilley, T. et al. Deciphering the genetic regulation of peripheral blood transcriptome in pigs through expression genome-wide association study and allele-specific expression analysis. *BMC Genomics* **18**, 967 (2017).

14. Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence data. *Bioinformatics* **30**, 2114–2120 (2014).

15. Yates, A. et al. Ensembl 2016. *Nucleic Acids Res.* **44**, D710–716 (2016).
16. Anders, S., Pyl, P. T. & Huber, W. HTSeq-a Python framework to work with high-throughput sequencing data. *Bioinformatics* **31**, 166–169 (2015).

17. Robinson, M. D., McCarthy, D. J. & Smyth, G. K. edgeR: a bioconductor package for differential expression analysis of digital gene expression data. *Bioinformatics* **26**, 139–140 (2010).

18. Robinson, M. D. & Oshlack, A. A scaling normalization method for differential expression analysis of RNA-seq data. *Genome Biology* **11**, R25 (2010).

19. Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. *Bioinformatics* **29**, 15–21 (2013).

20. Van der Auwera, G. A. et al. From FastQ data to high confidence variant calls: the Genome Analysis Toolkit best practices pipeline. *Curr Protoc Bioinformatics* **43**(11), 10.1–33 (2013).

21. Danecek, P. et al. The variant call format and VCFtools. *Bioinformatics* **27**, 2156–2158 (2011).

22. Harvey, C. T. et al. QuASAR: quantitative allele-specific analysis of reads. *Bioinformatics* **31**, 1235–1242 (2015).

23. Altschul, S. F. et al. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. *Nucleic Acids Res* **25**, 3389–3402 (1997).

24. Eden, E., Navon, R., Steinfeld, I., Lipson, D. & Yakhini, Z. GOrilla: a tool for discovery and visualization of enriched GO terms in ranked gene lists. *BMC Bioinformatics* **10**, 48 (2009).

25. Kanehisa, M. & Goto, S. KEGG: kyoto encyclopedia of genes and genomes. *Nucleic Acids Res.* **28**, 27–30 (2000).

26. Vandesompele, J. et al. Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. *Genome Biol.* **3**, RESEARCH0034 (2002).

27. Meurens, F. et al. Commensal bacteria and expression of two major intestinal chemokines, TECK/CCL25 and MEC/CCL28, and their receptors. *PLoS ONE* **2**, e677 (2007).

28. Meurens, F. et al. Early immune response following Salmonella enterica subspecies enterica serovar Typhimurium infection in porcine jejunal gut loops. *Vet. Res.* **40**, 5 (2009).

29. Bruel, T. et al. Epithelial induction of porcine suppressor of cytokine signaling 2 (SOCS2) gene expression in response to Entamoeba histolytica. *Dev. Comp. Immunol.* **34**, 562–571 (2010).

30. Livak, K. J. & Schmittgen, T. D. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. *Methods* **25**, 402–408 (2001).

31. Lunney, J. K. Advances in swine biomedical model genomics. *Int. J. Biol. Sci.* **3**, 179–184 (2007).

32. Rothkötter, H.-J. Anatomical particularities of the porcine immune system-a physician's view. *Dev. Comp. Immunol.* **33**, 267–272 (2009).

33. Boeker, M., Pabst, R. & Rothkötter, H. J. Quantification of B, T and null lymphocyte subpopulations in the blood and lymphoid organs of the pig. *Immunobiology* **201**, 74–87 (1999).

34. Dvorak, C. M. T. et al. Gene discovery and expression profiling in porcine Peyer's patch. *Vet. Immunol. Immunopathol.* **105**, 301–315 (2005).

35. Machado, J. G., Hyland, K. A., Dvorak, C. M. T. & Murtaugh, M. P. Gene expression profiling of jejunal Peyer's patches in juvenile and adult pigs. *Mamm. Genome* **16**, 599–612 (2005).

36. Rothkötter, H. J., Huber, T., Barman, N. N. & Pabst, R. Lymphoid cells in afferent and efferent intestinal lymph: lymphocyte subpopulations and cell migration. *Clin. Exp. Immunol.* **92**, 317–322 (1993).

37. Rothkötter, H. J., Ulbrich, H. & Pabst, R. The postnatal development of gut lamina propria lymphocytes: number, proliferation, and T and B cell subsets in conventional and germ-free pigs. *Pediatr. Res.* **29**, 237–242 (1991).

38. Okada, T. et al. Chemokine requirements for B cell entry to lymph nodes and Peyer's patches. *J. Exp. Med.* **196**, 65–75 (2002).

39. Agace, W. W. & McCoy, K. D. Regionalized Development and Maintenance of the Intestinal Adaptive Immune Landscape. *Immunity* **46**, 532–548 (2017).

40. Jago, C. B., Yates, J., Câmara, N. O. S., Lechler, R. I. & Lombardi, G. Differential expression of CTLA-4 among T cell subsets. *Clin. Exp. Immunol.* **136**, 463–471 (2004).

41. Zabel, B. A. et al. Human G protein-coupled receptor GPR-9-6/CC chemokine receptor 9 is selectively expressed on intestinal homing T lymphocytes, mucosal lymphocytes, and thymocytes and is required for thymus-expressed chemokine-mediated chemotaxis. *J. Exp. Med.* **190**, 1241–1256 (1999).

42. O'Shea, J. J. & Paul, W. E. Mechanisms underlying lineage commitment and plasticity of helper CD4+ T cells. *Science* **327**, 1098–1102 (2010).

43. Ma, C. S. & Phan, T. G. Here, there and everywhere: T follicular helper cells on the move. *Immunology* **152**, 382–387 (2017).

44. Ramiscal, R. R. & Vinuesa, C. G. T-cell subsets in the germinal center. *Immunol. Rev.* **252**, 146–155 (2013).

45. Tsuji, M. et al. Preferential generation of follicular B helper T cells from Foxp3+ T cells in gut Peyer's patches. *Science* **323**, 1488–1492 (2009).

46. Perruzza, L. et al. T Follicular Helper Cells Promote a Beneficial Gut Ecosystem for Host Metabolic Homeostasis by Sensing Microbiota-Derived Extracellular ATP. *Cell Rep* **18**, 2566–2575 (2017).

47. Barman, N. N., Bianchi, A. T., Zwart, R. J., Pabst, R. & Rothkötter, H. J. Jejunal and ileal Peyer's patches in pigs differ in their postnatal development. *Anat. Embryol.* **195**, 41–50 (1997).

48. Müller-Schoop, J. W. & Good, R. A. Functional studies of Peyer's patches: evidence for their participation in intestinal immune responses. *J. Immunol.* **114**, 1757–1760 (1975).

49. He, Y., Vogelstein, B., Velculescu, V. E., Papadopoulos, N. & Kinzler, K. W. The antisense transcriptomes of human cells. *Science* **322**, 1855–1857 (2008).

50. Chen, C. et al. A genome-wide investigation of expression characteristics of natural antisense transcripts in liver and muscle samples of pigs. *PLoS ONE* **7**, e52433 (2012).

51. Chamberlain, A. J. et al. Extensive variation between tissues in allele specific expression in an outbred mammal. *BMC Genomics* **16**, 993 (2015).

52. Crowley, J. J. et al. Analyses of allele-specific gene expression in highly divergent mouse crosses identifies pervasive allelic imbalance. *Nat. Genet.* **47**, 353–360 (2015).

53. GTEx Consortium. Human genomics. The Genotype-Tissue Expression (GTEx) pilot analysis: multitissue gene regulation in humans. *Science* **348**, 648–660 (2015).

54. GTEx Consortium. The Genotype-Tissue Expression (GTEx) project. *Nat. Genet.* **45**, 580–585 (2013).

55. GTEx, C. et al. Genetic effects on gene expression across human tissues. *Nature* **550**, 204–213 (2017).

## Acknowledgements

We are grateful to the Genotoul bioinformatics platform at INRA Toulouse Midi-Pyrénées (http://bioinfo/genotoul.fr) for providing computing and storage resources. We are grateful to the INRA MIGALE bioinformatics platform (http://migale.jouy.inra.fr) for help and support. We acknowledge the INRA @bridge platform for the secure storage of samples and the team at Le Magneraud (INRA, GENESI) for their implication on the generation of animals and samples. We acknowledge Isabelle Oswald (INRA, TOXALIM, Toulouse) who contributed to the design of the SUS-FLORA project. TM was financially supported by a PhD fellowship co-funded by the INRA
Animal Genetics Division and BIOPORC (association composed of AXIOM, CHOICE GENETICS FRANCE, NUCLEUS and IFIP). Experiments were funded by the SUS-FLORA project (ANR-10-GENM-016) of the French National Agency and internal funds of INRA's Animal Genetics Division. The authors declare they have no competing interests.

**Author Contributions**

C.R.G. coordinated the SUS-FLORA project. J.E., C.R.G., M.B. and F.M. designed the experiments. T.M. performed all computational analyses under the supervision of J.E., A.R. contributed to biostatistical analyses. T.M., J.E. and C.R.G. interpreted all results and drafted the manuscript. M.B., B.B. and F.M. contributed to result interpretation. C.R.G., J.E., M.B., F.M., C.C., S.M., S.V.N., J.J.L., P.L., Y.B. and M.J.M. contributed to animal sampling. M.B., C.C. and S.M. performed RNA extractions and qRT-PCR assays. G.L. and J.L.C. performed the validation by pyrosequencing. D.E. performed RNA-Seq at INRA's Get-PlaGe sequencing platform (http://get.genotoul.fr/en/). All authors have read, improved and approved the manuscript.

**Additional Information**

**Supplementary information** accompanies this paper at https://doi.org/10.1038/s41598-018-27019-7.

**Competing Interests:** The authors declare no competing interests.

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_8.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International

License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2018