Aire Gene Influences the Length of the 3′ UTR of mRNAs in Medullary Thymic Epithelial Cells

Aire is a transcriptional controller in medullary thymic epithelial cells (mTECs) modulating a set of peripheral tissue antigens (PTAs) and non-PTA mRNAs as well as miRNAs. Even miRNAs exerting posttranscriptional control of mRNAs in mTECs, the composition of miRNA-mRNA networks may differ. Under reduction in Aire expression, networks exhibited greater miRNA diversity controlling mRNAs. Variations in the number of 3'UTR binding sites of Aire-dependent mRNAs may represent a crucial factor that influence the miRNA interaction. To test this hypothesis, we analyzed through bioinformatics the length of 3'UTRs of a large set of Aire-dependent mRNAs. The data were obtained from existing RNA-seq of mTECs of wild type or Aire-knockout (KO) mice. We used computational algorithms as FASTQC, STAR and HTSEQ for sequence alignment and counting reads, DESEQ2 for the differential expression, 3USS for the alternative 3'UTRs and TAPAS for the alternative polyadenylation sites. We identified 152 differentially expressed mRNAs between these samples comprising those that encode PTAs as well as transcription regulators. In Aire KO mTECs, most of these mRNAs featured an increase in the length of their 3'UTRs originating additional miRNA binding sites and new miRNA controllers. Results from the in silico analysis were statistically significant and the predicted miRNA-mRNA interactions were thermodynamically stable. Even with no in vivo or in vitro experiments, they were adequate to show that lack of Aire in mTECs might favor the downregulation of PTA mRNAs and transcription regulators via miRNA control. This could unbalance the overall transcriptional activity in mTECs and thus the self-representation.


INTRODUCTION
Autoimmune regulator (Aire) gene encodes a protein that play its role as a transcriptional modulator rather than a classical transcription factor in medullary thymic epithelial cells (mTECs). AIRE protein is considered a non-classical transcription factor due two characteristics found in its mode of action. Firstly, it associates with several protein partners in pushing RNA Pol II on the chromatin, favoring the elongation phase of transcription, and second due to the large set of downstream genes it controls (1)(2)(3)(4)(5)(6).
In addition to AIRE, Forebrain embryonic zinc finger-like protein 2 (FEZF2) also play its role in mTECs but as a classical transcription factor (7)(8)(9). Most of genes controlled by AIRE encode peripheral tissue antigens (PTAs) that once expressed by mTECs ensure the self-representation in the thymus triggering the elimination of auto-reactive thymocytes though negative selection. Due to the diversity of autoantigens represented, this phenomenon is referred to as promiscuous gene expression (PGE) (3)(4)(5)(10)(11)(12). AIRE is considered the main controller of PGE/self-tolerance in the thymus and its disruption caused by mutations results in an aggressive autoimmune disease referred to as APECED or APS-1 (13)(14)(15).
The meaning of the Aire expression is relative to the immunological tolerance that occurs in the thymus, since the expression of this gene is preferentially in that organ. Other organs such as pancreas, adrenal cortex, testes, lymph nodes, fetal liver and appendix also express Aire, without however knowing its exact function in these organs (6,(16)(17)(18)(19).
The human AIRE gene maps on chromosome 21q22.3 and in the mouse on chromosome 10 position 39.72 cM and the primary structure of Aire gene suggests that the distribution of its domains encodes a protein that interacts with the chromatin and regulates transcriptional processes (2,(20)(21)(22)(23)(24)(25)(26)(27).
Nevertheless, AIRE protein does not act as a classical transcription factor, which binds directly to gene promoters' consensus sequences. Instead of, AIRE forms a multimolecular complex in which participates the so termed AIRE-partners (CBPBP, SIRT-1, DNA-PK, PARP-1, DNA-TOP2a, BRD4, P-TEFb) that work together in the releasing of stalled RNA-Pol II in the chromatin (4,5,28). Not yet fully resolved if acetylation of lysine residues of AIRE by CREB-binding protein (CBPBP) stabilizes intracellular localization of AIRE protein (29,30), but group-III histone-deacetylase SIRT-1 preserves AIRE-dependent PGE by deacetylation of such residues (31,32).
Of note, the DNA-topoisomerase (DNA-TOP2a) protein, which is a strong AIRE partner acts on the DNA topology by removing positive and negative supercoils when it causes transient DNA breaks. DNA breaks result in a local relaxation of chromatin that facilitates the initiation and post-initiation of gene transcription (6,33) and attracts DNA repair elements such as DNA-PK and PARP1 (6).
However, it is the interaction with PTEFb that resolves that AIRE play its role in the post-initiation phase of transcription (34). As previously discussed (6), gene transcription in eukaryotic cells is annulled if PTEFb does not proceed with the elongation and splicing of pre-mRNA into mature mRNA and subsequent release of stalled RNA Pol II (35). In turn, bromodomaincontaining protein 4 (BRD4) interacts with PTEFb being that balanced phosphorylation and acetylation of AIRE are needed to preserve such interaction (36).
Evidence suggests that AIRE and its partners bind to chromatin super-enhancers, which include the transcription start sites of most Aire-dependent genes (37). This multimolecular complex enables the transcription elongation step of a myriad of genes including those encoding autoantigens (4,6,28).
AIRE also associates to proteins involved with pre-mRNA processing (1), alternative splicing (38) and polyadenylation factors (HNRNPL/F, SFRS1, and CLP1) (39). It is plausible that this might be directly or indirectly related to the processing of Aire-dependent mRNAs in mTECs.
By means of classical gene knockout (KO) (40)(41)(42) or gene KO through Crispr-Cas9 (43) it was demonstrated that Aire affects the abundance of mRNA isoforms in mTECs. Previously we observed that miRNAs exert posttranscriptional effect on Airedependent mRNAs (44). However, this action is not uniform and a greater set of miRNAs exert their effect on mRNAs when Aire is down-regulated in mTECs (45).
These findings prompted us to investigate whether Aire could affect isoforms of mRNAs involving the length of 3'UTR end as we have previously suspected (3). Whether this process exists among the Aire-dependent mRNAs, this may be influencing the control of the PGE and consequently of the self-representation in the thymus. Accordingly, the aim of this study was to evaluate whether Aire affects the length of the 3'UTR region of mRNAs it controls.
To test this we analyzed the existing RNA-seq data of mTECs isolated from Aire wild-type or Aire-knockout (KO) mice (41) identifying and quantifying those mRNAs differentially expressed between these two samples and that showed changes in the length of the 3'UTR region and thus in the miRNA binding sites.

RNA-Seq Data Analysis
We obtained the sequences from raw RNA-seq data deposited in the Gene Expression Omnibus GEO repository (https://www. ncbi.nlm.nih.gov/geo/) and available under accession number GSE65617 (41), which consisted of RNA-seq of CD45 − EPCAM + sorted as cTECs (Ly51 + UEA1 − ), Ly51 − UEA + CD80 lo MHC-II lo sorted as mTEC lo , Ly51 − UEA + CD80 hi MHC-II hi sorted as mTEC hi . The analyzed sequences were from mTEC hi cells (Ly51 − UEA + CD80 hi MHC II hi ) isolated of wild type (WT) or Aire knockout (KO) C57BL/6 mice. Each biological replicate (n = 3 for each cell type) consisted of a pool of mTECs (FACS sorted) of four isogenic animals (two males and two females from the same litter) WT or KO. Sequencing was performed on an Illumina HiSeq 2000 RNA-seq platform with paired-end libraries with fragments of 100 bp (2 × 100 bp).
Analyses of RNA sequences were performed according to the following pipeline adopted for this study. The Fastq sequences provided by GEO were initially checked for quality through the FASTQC program (http://www.bioinformatics.babraham.ac. uk/projects/fastqc/). Then, the Fastq sequences were aligned against the Mus musculus reference genome (mm10) (Mouse Genome Reference Consortium https://www.ncbi.nlm.nih.gov/ grc/mouse) through the Spliced Transcripts Alignment to a Reference (STAR 2.5) program (https://github.com/alexdobin/ STAR) (46), which outputs a BAM file containing both the sequence and its genomic coordinate.
Together with a General Transfer Format (GTF) the file containing all the genomic annotations were used to determine the total number of reads per transcript using the HTSeq Count program (https://htseq.readthedocs.io/en/release_0.11.1/) (27).
Information about genome assembly and gene annotation corresponds to the mm10 reference genome that was submitted by the Genome Reference Consortium (mouse genome reference GRCm38). The complete description of this annotation as for example, base pair counts, gene counts among other information, can be found available at (m.ensembl.org/Mus_musculus/Info/Annotation).
For each sample analyzed, we obtained a list of genes with their respective counts, which were analyzed to determine the genes differentially expressed by using the DESEq2 package (https://bioconductor.org/) (47), which is carried out through the platform R (https://www.r-project.org/). A heat-map was constructed in order to hierarchically cluster samples and mRNAs based on Euclidean distance and the complete linkage method.
In this work, genes with a p < 0.05, Benjamini-Hochberg FDR correction and fold change ≥ 1.5, were considered as differentially expressed.

mRNA 3'UTR and Alternative Polyadenylation Site (APA) Analysis
For evaluation of changes in the 3'UTR sequences of mRNAs expressed in WT or Aire KO mTEC cells, we used the 3'UTR Sequence Seeker (3USS) (https://bio.tools/3uss) program. The 3USS program automatically identifies alternative 3'UTR sequences based on alternative polyadenylation sites (APA) obtained from actual RNA-seq data (long or short sequences) in relation to the mouse reference genome and transcriptome annotated in the UCSC (https://genome.ucsc.edu), NCBI (http:// www.ncbi.nlm.nih.gov/), Ensembl (http://www.ensembl.org/ index.html) and Gencode (http://www.gencodegenes.org/) databases. In addition, this program attempts to identify new 3'UTR sequences, which have not yet been annotated in public databases (48).
We then classify the 3'UTR alternative sequences into two categories, Aire-dependent or Aire-independent considering the differentially expressed genes between WT and Aire KO mTECs.
Next, we focused on the 3'UTR alternative sequences of Aire-dependent peripheral tissue antigens (PTAs) mRNAs. The identification of PTA mRNAs was done based on the classification of Sansom (40) that use both the HSGeneAtlas (https://dintor.eurac.edu/doc/doc/HSGeneAtlas.html) and the mTEC transcriptome obtained by these authors as a parameter. Alignment between the reference sequences and 3'UTR alternative sequences of mRNAs encoding Aire-dependent PTAs was done through the BLAST program (http://blast.ncbi.nlm. nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastHome) (49).
In addition, we use the TAPAS tool (https://github.com/ arefeen/TAPAS) for the determination of the number of APA sites from RNA-seq data, which works with two or more APA sites or with APA sites located prior to the last exon of a RNA target sequence (50). This tool compares RNA-seq data of pairs of biological samples that in this study were samples from Aire WT or Aire KO mTECs. The tool identify those differentially expressed APA sites as well as mRNAs that contain 3'UTRs with shortening/ lengthening events.

In silico Validation of miRNA-mRNA Interactions
The new and predicted miRNA-mRNA interactions were validated using the RNA-Hybrid algorithm (http://bibiserv. techfak.uni-bielefeld.de/rnahybrid) (51). It is a bioinformatics tool that analyzes the secondary structure, and determines the most favorable hybridization site between a given miRNA and its target mRNA. The methodology involves a dynamic programming algorithm that calculates the location for the most energy-efficient hybridization. This methodology is based on the thermodynamic property of miRNA-mRNA interaction, in which the RNA double strand is more stable, that is, the miRNA-mRNA binding is stronger when the energy is lower. The interaction force is determined by the quantification of the minimum free energy (mfe) normalized by the length of the target mRNA. In this study, we selected interactions with mfe ≤ −20 kcal / mol, which are considered thermodynamically stable.

Functional Enrichment of Non-PTA Aire-dependent mRNAs
To better understand the functional annotation of Airedependent mRNAs other than those encoding PTAs, but that featured alterations in the length of their 3 ′ UTRs, we use the DAVID (david.ncifcrf.gov/) considering a P-value EASE score ≤ 0.1 to examine the significance of gene-term enrichment with a modified Fisher's exact test. The false discovery rate (FDR) under certain rate (e.g., ≤ 0.05) provided by DAVID is calculated on the basis of the multiple testing correction as Bonferroni and Benjamini methods.

3'UTR Shortening / Lengthening Events Among the Differentially Expressed mRNAs
For identification of possible changes in 3 ′ UTR sequences of mRNAs expressed in mTEC cells, transcripts assembled from the WT (146,539 transcripts) or KO mTECs (84,945 transcripts) were used as input to the 3USS program. This program automatically identified alternative 3 ′ UTR sequences, using as reference the Mus musculus genome, version GRCm38/mm10. The set of transcripts assembled from the WT mTECs comprised 30,014 and the set from KO mTECs comprised 29,705 protein encoding mRNAs. Of the mRNAs expressed by the WT mTECs, 1,143 showed changes in their 3'UTRs (shortening / lengthening) and of the mRNAs expressed by the KO mTECs, 830 showed changes in their 3'UTRs ( Figure 1).
Then, WT or KO mTEC transcripts, which had the modified 3 ′ UTRs sequences (short or long 3 ′ UTRs) were regrouped into three categories: (1) those specific of WT mTECs (642 transcripts); (2) those specific of KO mTECs (359 transcripts); and (3) common between WT and KO mTECs (500 transcripts) (Figure 2). From these categories, we selected only the 359 KO specific transcripts of which 152 were differentially expressed when comparing WT vs. KO mTECs. Among these transcripts, 41 represent PTA mRNAs and 111 non-PTA mRNAs. This set was considered as differentially expressed Aire-dependent transcripts, which were hierarchically clustered and selected for further analysis (Figure 3, Table 1 and Supplementary Table 1). For further information about gene symbol, full name, chromosomal localization and biological processes, see Supplementary Table 2).
Next, we focused on the mRNAs that encode PTAs to evaluate the type of change that occurred in their 3'UTRs. We observed that most of them (39/41) showed increase in their 3'UTR sequences, ranging from +651 nt (Cxcr6 mRNA) to +7,724 nt (Rcvrn mRNA). Regarding the shortening, only Ttc9 and Nos1 PTA mRNAs showed reduction of their of the 3'UTR sequences corresponding to −2.004 and −2.136 nucleotide, respectively. Table 2 summarizes the nucleotide increase counts among transcripts analyzed.

Interaction Between miRNAs and Aire-dependent PTA mRNAs
In order to evaluate whether the increase or reduction in the PTA mRNA 3 ′ UTRs affected their predicted interactions with miRNAs, we used the miRDB database that initially evaluates the potential complementarity of a given miRNA seed region to the 3 ′ UTR of a target mRNA. Firstly, we identified the miRNAs that interact with the 3'UTR regions of PTA mRNAs of the  mouse reference transcriptome that is available in the GenBank NCBI (https://www.ncbi.nlm.nih.gov/genbank/) database, and these results were considered as a reference for further analysis. We then used the 3'UTR mRNA sequences from the mTEC RNAseq data to compare the same PTA mRNAs from this study that showed their increased or reduced 3' UTRs. This type of analysis allowed us to qualitatively and quantitatively evaluate the miRNA-mRNA interactions as well as to predict new interactions ( Table 3). All miRNA-mRNA interactions showed in Table 3 were validated in silico through the RNA-Hybrid program. The results revealed that the majority of interactions show perfect complementarity between the miRNA seed region and the 3'UTR sequence of the target mRNA with mfe ≤ −20 kcal / mol. Here we have exemplified only the miRNA-mRNA interactions of the two mRNA targets that showed reduction of the 3'UTR sequence i.e., the Nos1 and Ttc9 mRNAs and of the two targets that showed increase of their 3'UTRs, i.e., Cxcr6 and Rcvrn mRNAs, respectively (Supplementary Figures 1-4).
Supplementary Figure 1 show that when we compared the interaction between the seed region of miRNA 4661-5p or miRNA 6924-5p with the 3'UTR region of the GenBank reference mRNA Cxr6 with the result obtained with mTECs (this study) for Cxcr6 mRNA, we observed that increased 3 'UTR sequence of this PTA mRNA does not unbalance the miRNA-mRNA interactions. This increase favors energetically stable interactions presenting mfe of −25.6 kcal/mol (miR-4661-5p and Cxcr6) and−26 kcal/mol (miR-6924-5p and Cxcr6).
With respect to the Rcvrn mRNA, Supplementary Figure 2

PTA Self-Representation
The 41 mRNAs that encode PTAs and the featured 3'UTR alterations were assigned to 11 anatomic functional body systems as follows (in descending order): central nervous, glands, immune and lymphatic, reproductive, ocular, endocrine, digestive, renal and urinary, muscular, integumentary and respiratory (Figure 4).

Functional Enrichment of Non-PTA Aire-dependent mRNAs
We performed the functional enrichment of mRNAs that do not encode PTAs, but which are Aire-dependent and which have altered the size of the 3 ′ UTR region ( Table 2). The biological process most represented by these mRNAs is relative to the regulation of transcription (Figure 5) including the transcription factors (Tbx6, Zfpm1, Crebbp, Olig2, Ikzf3, Kdm4d, Ar, Zfp605, Pura, and Zfp382).

Identification of Alternative Polyadenylation (APA) Sites in
Aire-dependent PTA-mRNAs Table 4 shows the positioning of APA sites throughout the nucleotide sequence of the Aire-dependent PTA mRNAs. Those known APA sites for each mRNA, those most frequent and those APA sites not previously described, but which were identified in the mTEC cell samples studied, could be located. We observed that the number of APA sites does not change significantly between WT and Aire KO mTEC cells. However, we observed new APA sites in PTA mRNAs expressed by this cell type.

DISCUSSION
In this study, we analyzed two large groups of mRNAs that were expressed in medullary thymic epithelial cells (mTECs) isolated from autoimmune regulator (Aire) wild type (WT) or Aireknockout (KO) mice (41). The first group of mRNAs encode peripheral tissue antigens (PTAs), and the other encode proteins involved with transcription control including transcription factors (TFs). The analysis was done to assess the hypothesis that Aire regulates the mRNA 3'UTR length, the preferential region where are found the binding sites of miRNAs.
In previous studies, we have observed that Aire regulates a set of miRNAs in mTECs (52), which establish miRNA-mRNA network and regulate a set of downstream Aire-dependent mRNAs (45). Interestingly, reduction in Aire in these cells imbalance miRNA-mRNA interactions and a larger set of miRNAs now act on the target mRNAs (45).
These observations raised questions about the factor (s) that lead to increased miRNA diversity acting on mRNAs (3,44,45). In our view, factors such as stoichiometry between miRNAs and mRNA targets or stability of base pairing may have influenced the miRNA interaction. The stoichiometry between these RNA species is expected to vary during the cell cycle and does not  explain the augmentation of miRNA diversity that interact with their mRNA targets.
The Gibbs free energy or minimum free energy (mfe) ≤ −20 is considered optimal for maintaining miRNA-mRNA pairing and this value is widely used as a cutoff in the search for stable miRNA interactions. Even so, different miRNAs could still play a role at different mfe values (51). In our view, a plausible explanation for the increase in the number of miRNAs that interact with target mRNAs would be the modification in the target, most likely in the 3'UTR region. To act, these miRNAs would have to have more binding sites in the target 3'UTR region.
Variations in the structure of Aire-dependent PTA mRNAs involving the number of 3'UTR miRNA binding sites could be the most obvious cause influencing the greater variety of controller miRNAs in Aire-downregulated mTEC cells. Variations in the number of 3'UTR binding sites of Aire-dependent PTA mRNAs may represent a possible factor that influence the miRNA interaction (3,52) an aspect that has been pursued in our group and independently by the group of M. Giraud (C. Guyon) at the INSERM (U1016) in Paris and at the University of Nantes, France (3).
As we previously asked (3), what would be the consequences of miRNA action on these PTAs? Could the different miRNAs in mTECs control autoantigen protein synthesis and consequently control the self-tolerance induction? Could changes in the length of mRNA 3 ′ UTRs influence the action of miRNAs? Previously, we stated that there might be changes in length of the 3 ′ UTR sequence of Aire-dependent mRNAs expressed in mTECs (3,45,52).
To test this hypothesis, we analyzed the length of 3 ′ UTRs of a large set of Aire-dependent mRNAs whose data were obtained from existing RNA-seq data of mTECs isolated of wild type or Aire-knockout (KO) mice. In addition to PTAs, we also analyze the variations of the 3 ′ UTRs of mRNAs that encode transcription factors, as these elements are at the basis of transcriptional gene expression control, which may influence self-representation by mTECs.
The use of the 3USS software allowed us to identify the existence of alternative 3'UTR sequences when comparing the mTEC wild type cell transcript set with the Aire KO mTECs. Both samples showed shortening or lengthening in this region and this demonstrated that Aire influences the generation of mRNA isoforms with respect to the 3'UTR region.
The implication of this is associated to immune homeostasis, which is largely associated to self-non-self discrimination by T-cells, a process that occurs in the thymic stroma during the development of thymocytes. These cells in population are exposed to a myriad of peripheral tissue antigen (PTA) peptides that are expressed and presented via MHC-II by medullary thymic epithelial cells (mTECs) to thymocytes (3-5, 10, 11). The greater the diversity of autoantigens presented by mTECs, the more comprehensive and at the same time accurate will be the self-tolerance of the T-cell population able to migrate from thymus to the periphery (3,53,54). This process must be finely controlled, if not, loss of tolerance may occur with consequences in the emergence of autoimmune reactions.
For example, our results showed that in the absence of Aire, the mRNA encoding the PTA Rcvnr, increased by 7,729 nucleotides in its 3'UTR and, consequently, the generation of two new miRNA binding sites (for miR-7076-5p and miR−7075-5p). As this PTA represents retinal pigment epithelium, retina and iris, which are target during autoimmune reactivity, this suggests a mechanism of autoimmunity of the eye structures. Similarly, the mRNA that encodes the PTA Syde-1 had an increase of 5,800 nucleotides in its 3'UTR with the appearance of a new miR-1187 binding site. This PTA represents ovary and adrenal glands that are target organs for autoimmune reactivity in APECED syndrome. This illustrates in what way changes in the 3'UTR of mRNAs could clinically affect aggressive autoimmunity.
Two other mRNAs that encode proteins involved in the biology of mTECs, such as Crebbp (an AIRE partner) and Ago-2 (a key element in the RISC complex during the miRNA mechanism of action), also had an increase in their 3 ′ UTRs and consequently the generation of new miRNA binding sites. CREBBP being negatively regulated by miRNAs may destabilize the AIRE complex in the chromatin of mTECs with consequences for PGE control. In addition, Ago-2 being negatively regulated might cause a general imbalance in the posttranscriptional control.
Aire gene/protein is the main transcriptional regulator of PTAs that acts in synergy with Forebrain embryonic zinc fingerlike protein 2 (FEZF2) in mTECs (4)(5)(6)(7)(8)(9). Although Aire is an upstream transcriptional controller, this gene / protein has a link to posttranscriptional control in mTECs (11,45,52). In addition to evaluating variations in the length of 3'UTR sequences of the Aire-dependent mRNAs, we also examined whether these variations could affect miRNA-mRNA interactions.
Initially, it was necessary to establish a reference for comparing miRNA interactions with known target PTA mRNAs by making use of miRDB database. This approach enable us to compare and to evaluate the effect of Aire on both the length of mRNA 3'UTRs and the number of different miRNAs they interact, i.e., we treated the data qualitatively and quantitatively. We use as a criterion for thermodynamic stability of the interactions between the seed regions of the miRNAs and the 3'UTR of the target mRNAs the score of minimal free energy (mfe ≤ −20). This criterion allowed the development of a thermodynamically favorable miRNA-mRNA interaction prediction algorithm (51), which is a widely used method for different model-systems.
We show that Aire interferes with the 3'UTR length of Airedependent mRNAs and that its absence modify the composition of miRNA-mRNA interactions creating new miRNA sites on PTA mRNAs. This might have influence in reduced PTA protein translation. We evidence that absence of Aire not only interferes directly with transcription of PTAs, but also with the structure of its respective mRNAs.
As previously shown (41,43,55) Aire also controls the expression of non-PTA mRNAs. In this study, it was found that those mRNAs associated with transcription control, including those encoding transcription factors, were also under Aire control that affected the size of their respective 3'UTRs. This certainly leads to an imbalance in the overall gene expression of mTECs with implications for self-representation. Although FEZF2 plays its role as a TF in mTECs that works in synergy to AIRE, we were not able to evidentiate its transcriptional expression in the samples analyzed.
It is well-known that the maturation of the 3 ′ ends of most eukaryotic mRNAs is dependent on cleavage and polyadenylation (C/P) and that this RNA species has several C/P sites. Processing of mRNAs by C/P, which involves endonucleolytic cleavage of the nascent RNA and then addition of the poly (A) tail is one of the processes of control of gene expression in eukaryotes (56)(57)(58)(59)(60).
The site where polyadenylation occurs, the poly(A) site, is so important that mutations at this region are associated with human diseases such as thrombophilia to some forms of thalassemias and cancer as well (57,59,61,62). In addition, shortening or lengthening and polymorphisms of 3'untranslated region (3'UTR) corresponds to a source of mRNA isoforms and this process is associated to some human diseases including autoimmune diseases (60,(63)(64)(65).
However, under physiological conditions eukaryotic genes express isoforms that show differences in the 3' end of their respective mRNAs. The mRNAs from a given gene but differing at their 3' end are translated into proteins that show structural and functional differences. In addition, 3' end processing may involve the half-life of the transcribed mRNA due to inclusion or exclusion of regulatory elements, as a rule miRNA binding sites. The absence or presence of these elements implies the biological properties of mRNA/protein such as export, abundance, stability, location, and transport and translation rate (59,62).
AIRE protein is also involved with these processes either interacting with proteins of pre-mRNA processing (1), alternative splicing (38) and polyadenylation factors (39). Then, we investigated whether Aire could affect the number of APA sites of Aire-dependent mRNAs. Our results show that rather reflecting the presence or absence of Aire, the new APA sites in mRNAs is more associated to mTEC cell type. This is still an open question to search for genes/proteins involved in this process.
Overall, these results point toward the existence of a controlling mechanism in mTECs in which Aire gene/protein affects the length of 3'UTR of Aire-dependent mRNAs. Under absence of Aire, the mRNA targets feature increased sequence were is found miRNA binding sites. The miRNAs then exert posttranscriptional control over these mRNAs and might decrease protein translation that in the case of PTAs influence self-representation in the thymus.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the GEO acc number GSE65617.

AUTHOR CONTRIBUTIONS
EO: conceived the study and performed data analysis. AA: performed data analysis. CS-H: performed data analysis. MD: performed data analysis. GP: conceived the study, chose the analysis tools, and wrote the manuscript.

ACKNOWLEDGMENTS
We thank Mrs Aline F. Tiballi (School of Dentistry of Ribeirão Preto, USP, Ribeirão Preto, Brazil) for technical assistance.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu. 2020.01039/full#supplementary-material Supplemental Figure 1 | Molecular interactions between miRNAs and mRNA targets that feature reduction or increase in their 3'UTRs and their respective thermodynamic minimal free energy (mfe). Comparison of the interaction between the seed region of miRNA 4661-5p or miRNA 6924-5p with the reference 3'UTR (GenBank NCBI) of mRNA Cxr6 with the respective interaction obtained with mTECs (this study) for Cxcr6 mRNA.
Supplemental Figure 2 | Molecular interactions between miRNAs and mRNA targets that feature reduction or increase in their 3'UTRs and their respective thermodynamic minimal free energy (mfe). The increase of 3'UTR of mRNA Rcvrn observed in mTECs created new interactions for miRNA 7076-5p and miRNA 7075-5p.
Supplemental Figure 3 | Molecular interactions between miRNAs and mRNA targets that feature reduction or increase in their 3'UTRs and their respective thermodynamic minimal free energy (mfe). The reduction of Nos1 mRNA 3'UTR in mTECs created interaction with two new miRNAs, miRNA-7073-5p and miRNA-6992-5p.
Supplemental Figure 4 | Molecular interactions between miRNAs and mRNA targets that feature reduction or increase in their 3'UTRs and their respective thermodynamic minimal free energy (mfe). The 3'UTR reduction of the Ttc9 mRNA observed in mTECs did not alter the interaction with miRNA 3057-30, but created two new possibilities of interaction involving miRNA-705 and miRNA-7008-5p.