Human Ubiquitin-Specific Peptidase 18 Is Regulated by microRNAs via the 3'Untranslated Region, A Sequence Duplicated in Long Intergenic Non-coding RNA Genes Residing in chr22q11.21

Ubiquitin-specific peptidase 18 (USP18) acts as gatekeeper of type I interferon (IFN) responses by binding to the IFN receptor subunit IFNAR2 and preventing activation of the downstream JAK/STAT pathway. In any given cell type, the level of USP18 is a key determinant of the output of IFN-stimulated transcripts. How the baseline level of USP18 is finely tuned in different cell types remains ill defined. Here, we identified microRNAs (miRNAs) that efficiently target USP18 through binding to the 3’untranslated region (3’UTR). Among these, three miRNAs are particularly enriched in circulating monocytes which exhibit low baseline USP18. Intriguingly, the USP18 3’UTR sequence is duplicated in human and chimpanzee genomes. In humans, four USP18 3’UTR copies were previously found to be embedded in long intergenic non-coding (linc) RNA genes residing in chr22q11.21 and known as FAM247A-D. Here, we further characterized their sequence and measured their expression profile in human tissues. Importantly, we describe an additional lincRNA bearing USP18 3’UTR (here linc-UR-B1) that is expressed only in testis. RNA-seq data analyses from testicular cell subsets revealed a positive correlation between linc-UR-B1 and USP18 expression in spermatocytes and spermatids. Overall, our findings uncover a set of miRNAs and lincRNAs, which may be part of a network evolved to fine-tune baseline USP18, particularly in cell types where IFN responsiveness needs to be tightly controlled.


INTRODUCTION
Ubiquitin-specific peptidase 18 (USP18) is an interferon (IFN)stimulated gene (ISG) exerting a specific and non-redundant role in the negative regulation of type I IFN (here, IFN) responses (Honke et al., 2016;Basters et al., 2018). USP18 is also an isopeptidase (de-ISGylase) that non-redundantly cleaves the ubiquitin-like ISG15 from protein conjugates (Malakhova et al., 2006;Francois-Newton et al., 2011;Arimoto et al., 2017). USP18 contains a conserved cysteine protease catalytic domain that is flanked by short N-and C-terminal regions providing key binding surfaces. Through the N-ter region (amino acids 36-51) and the C-ter region (amino acids 313-371, encoded by exons 9-11), human USP18 is recruited to the IFNAR2 receptor subunit at the plasma membrane. This interaction is reinforced by adjacent regions binding to STAT2 (Arimoto et al., 2017). The complex interferes with JAK1/TYK2 activation and attenuates STAT-mediated IFN-stimulated genes ISGs induction. The role of the C-ter region is supported by the dysregulated IFN response phenotype of the Usp18 Ity9/Ity9 mouse strain bearing a single homozygous mutation (Leu361Phe) in this region (Dauphinee et al., 2014). Deficiency of USP18 causes a severe type I interferonopathy resulting in perinatal death with serious brain malformations due to spontaneous microglia activation, which most likely results from unrestrained response to constitutive IFNβ (Meuwissen et al., 2016). Indeed, high baseline USP18 maintains microglia quiescence and prevents sub-threshold activation (Goldmann et al., 2015;Schwabenland et al., 2019). In mouse models of VSV and LCMV infection, it was shown that the ability of CD169+ spleen macrophages and dendritic cells to present viral antigens and elicit an innate and adaptive immunity relies on the expression of Usp18 (Honke et al., 2011(Honke et al., , 2013. Altogether, these findings point to USP18 as a key determinant of cell responsiveness to IFN, including in the context of constitutive low IFNβ levels. We and others have shown that in humans, but not in mice, free ISG15 sustains the level of USP18, by preventing its proteasomal degradation by the S-phase kinase-associated protein 2 (SKP2; Zhang et al., 2015;Speer et al., 2016;Vuillier et al., 2019). Yet, little is known on how baseline USP18 levels are set in different cell lineages.
MicroRNAs (miRNAs) act as fine-tuners of gene expression and protein output. Each gene target can be regulated by numerous miRNAs, which typically bind to specific sites in 3' untranslated regions (3'UTRs) of messenger RNAs (mRNAs), influencing mRNA stability and translation. Hence, the amount of target varies from cell to cell depending on the identity and the expression level of each miRNA (Bartel, 2004). MiRNAs operate in different physiological and pathological processes (Ebert and Sharp, 2012) and also participate to the complex network that regulates immune cell development, function, and response to stimuli (Lindsay, 2008). As for many immune pathways, miRNAs have been shown to regulate the IFN system both at the level of production and response. Hence, some of the key players of the type I IFN signaling pathway are under the control of miRNAs (Forster et al., 2015), but no miRNA targeting USP18 has been described so far.
Here, we investigated whether miRNAs contribute to tune down baseline USP18 in a cell-context specific manner. We identified four miRNAs that negatively regulate USP18 through binding the 3'UTR. Several copies of the USP18 3'UTR are present in the human genome. These duplications reside on chr22q11.21 and are embedded in long intergenic non-coding (linc) RNA (lincRNA) genes that we found expressed in different human tissues. In particular, we describe a lincRNA bearing USP18 3'UTR with a remarkable cell type-specific expression in testicular germ cells. The positive correlation of the expression of this lincRNA with USP18 raises the possibility that the USP18 3'UTR not only acts as a cis-regulatory sequence targeted by miRNAs, but also acts in trans when embedded in lincRNA molecules expressed in particular cell types.

MATERIALS AND METHODS
Computational Prediction of miRNAs Targeting the USP18 3'UTR Bibiserv 1 and miRWalk 2.0 2 interfaces were used. The human USP18 3'UTR sequence [580 nucleotide (nt)] was downloaded from the UCSC genome browser (hg38; https://genome.ucsc. edu/) and used as input. In Bibiserv, RNA hybrid algorithm was used with the following parameters: energy threshold = −20, no G:U in the seed, helix constraint from nt 2 to 8 (7-mer seed), and approximate p-value for 3utr_human. In miRWalk, TargetScan, Pita, miRDB, and RNA22 algorithms were used with the following parameters: minimum seed length 7, position 2 (nt 2) as start position of the miRNA seed, and p value 0.05.

Testicular Cells Isolation
Normal human adult testes were obtained either after orchidectomy from prostate cancer patients who had not received any hormone therapy or at autopsy. The procedure was approved by Ethics Commitee Ouest V, Rennes, France (authorization DC-2016-2783  and stored at −80°C until RNA extraction. Pachytene spermatocytes, round spermatids, Leydig cells, and peritubular cells were isolated and cultured according to previously described procedures (Rolland et al., 2019) before freezing at −80°C for RNA extraction. Primary Sertoli cells were purchased from Lonza (Walkersville, MD, United States) and cultured as previously described (Rolland et al., 2019).

RNA Extraction and RT-qPCR Quantification
Total RNA was extracted from Trizol lysates of freshly isolated monocytes or PBL, testis explants, and testicular populations, using miRNeasy mini kit (Qiagen, Germantown, MD, United States) and following the manufacturer's instructions. Total liver RNA was purchased from Thermo Fisher scientific. Quantification and purity were assessed by a Nanodrop spectrophotometer (Nanodrop2000, Thermo Fisher Scientific  -191-5p, miR-24-3p, miR-423-5p, and miR-532-3p were purchased from Horizon Discovery Ltd. A negative control -miRIDIAN miRNA Mimic Negative Control #1 -(Horizon Discovery Ltd) was used to assess the specificity of the effect driven by the specific miRNA sequences. Mimic reverse transfection (50 nM) was performed in HeLa S3 cells cultured in DMEM supplemented with 10% heat inactivated FCS and antibiotics. Transfection was carried out by using Lipofectamine RNAiMAX (Thermo Fisher Scientific) according to the manufacturer's instructions. Briefly, the mimic-RNAiMAX complexes were prepared in Opti-MEM serum-free medium by mixing 50 nM of mimic with a dilution of 1:1,000 of lipofectamine RNAiMAX (calculated on the final volume of the culture). After 15 min at room temperature, the complexes were distributed in a p24-(for RNA isolation and luciferase assays) or a p60-(for protein extraction) well plate and 1 × 10 5 or 9 × 10 5 cells, respectively, and were plated on top of the complex. After 48 h of transfection, cells were harvested and lysed for RNA isolation, luciferase assay, or protein extraction.

Identification and Analysis of Duplicated USP18 Sequences
To identify USP18 copies, we performed BLAT 3 analysis using human USP18 gene sequence. The sequence was downloaded from UCSC genome browser (hg38; https://genome.ucsc.edu/). We selected only sequences with more than 95% identity to USP18 sequence and longer than 100 bp. To find whether the USP18 copies overlapped genes, we looked at the coordinates of these repeats on the UCSC genome browser (hg38). The genes overlapping USP18 copies were identified, and their genomic sequences and predicted mRNA sequences were downloaded. To measure the expression of FAM247A,C,D and linc-UR-B1, we designed specific primers, targeting the unique exon-exon11 junctions to avoid detection of USP18 mRNA. Next, we performed qPCR analysis as described above. All qPCR products were verified by sequencing.

Analysis of Public Available Datasets
Data on miRNA expression, shown in Figure 2 and Supplementary Figure S2 were retrieved from the FANTOM5 dataset 4 and from (Allantaz et al., 2012; Cohort Roche, GEO accession: GSE28492). The FANTOM5 dataset was downloaded and analyzed using Qlucore Omics Explorer software. RNA-seq datasets covering a panel of 32 tissues 5 were analyzed, and the results are shown in Supplementary Figure S5B. PolyA+ RNA-seq track from testis provided by the ENCODE project (GEO accession: GSM2453457-GSM2453458) and shown in Supplementary Figure S6A were visualized using autoscale mode in IGV browser. The data of the single cell RNA-seq on testis used in Figure 6 and Supplementary Figure S7 were retrieved from (Guo et al., 2018). Reads were aligned on the hg38 assembly of the human genome and only unique reads were filtered.
After 24 h, plasmid transfection was carried out using FuGENE HD (Promega) according to the manufacturer's instructions. Three biological replicates were prepared for each condition. Cells were lysed 24 h after transfection and analyzed with Dual-Luciferase ® Reporter Assay System (Promega) according to the manufacturer's instructions. Firefly luciferase was measured as normalizer.

Western Blot Analysis and Antibodies
Cells were lysed in modified RIPA buffer (50 mM Tris-HCl pH 8, 200 mM NaCl, 1% NP40, 0.5% DOC, 0.05% SDS, and 2 mM EDTA) with 100 mM PMSF and phosSTOP and a cocktail of antiproteases (Roche). A total of 40 μg proteins was separated by SDS-PAGE and analyzed by western blot. Membranes were cut horizontally according to molecular size markers, and stripes were incubated with different Abs. Immunoblots were analyzed with the ECL Western blotting Reagent (Pierce) or the more sensitive Western Lightning Chemiluminescence Reagent Plus (PerkinElmer), and bands were quantified with Fuji LAS-4000. For reprobing, blots were stripped in 0.2 M glycine (pH 2.5) for 20 min at rt. The following Abs were used: rabbit anti-USP18 (D4E7, Cell Signaling Technology, Beverly, MA, United States) and mouse anti-actin B (Sigma-Aldrich, St. Louis, MO, United States).

3' of cDNA Ends and PCR
3' of cDNA Ends (3'RACE) analysis was performed on testis and liver RNA using the 3'RACE System kit (Thermo Fisher Scientific) according to manufacturer's instructions. PCR was performed with Platinum SuperFi II DNA Polymerase (Thermo Fisher Scientific).

Primers
For the sequences of all primers used refer to Supplementary  Table S2.

Sequences
The sequences of human USP18 exon11, FAM247A/C/D, and a partial sequence of linc-UR-B1 are available as Supplementary sequences.

The 3'UTR of Human USP18 Is Targeted by at Least Four miRNAs
Using the bioinformatic prediction programs miRWalk and RNAhybrid, we identified 27 miRNAs predicted to target with high score the 580 nt-long 3'UTR of the human USP18 mRNA ( Figure 1A). USP18 is basally detected in most human tissues, albeit at variable levels (Supplementary Figure S1). We, therefore, analyzed the expression of the 27 miRNA candidates in the 90 cell types of the FANTOM5 dataset (de Rie et al., 2017). Fourteen miRNAs were found to be expressed at >50 counts per million mapped reads (CPM) in at least one cell type ( Figure 1A). We next tested the functional impact of these 14 miRNAs on USP18. For this, we expressed in HeLa S3 cells, each miRNA as mimic and measured the level of USP18 mRNA by qPCR. As shown in Figure 1B, seven miRNAs led to a reduction of endogenous USP18, four of them exhibiting a stronger effect. Consistently, this also resulted in lower USP18 protein levels, at least for those four miRNAs with stronger effects at the RNA level ( Figure 1C).
In order to validate the binding of each miRNA to the predicted site on USP18 3'UTR ( Figures 1D,E), this latter sequence was subcloned downstream of the Renilla luciferase gene in the psiCHECK2 vector and the impact of the four miRNAs on luciferase was measured in transfected HeLa S3 cells. As shown in Figure 1F, miR-191-5p, miR-24-3p, miR-423-5p, and miR-532-3p, expressed as mimics, significantly decreased luciferase activity with respect to the control mimic. Moreover, the effect of each mimic was abrogated when a 5 nt mutation was introduced in the seed-matched sequence of the predicted binding site on the 3'UTR ( Figure 1F). Altogether these data demonstrate that the USP18 3'UTR sequence can be directly targeted by miR-191-5p, miR-24-3p, miR-423-5p, and miR-532-3p, with an effect on USP18 abundance.
To determine in which cell type(s) this post-transcriptional control may operate, we performed a principal component analysis (PCA) based on the expression of the seven active miRNAs identified above in the 90 cell types (Figure 2A). The analysis revealed a different distribution of immune vs. non-immune cells, with the first three PC vectors explaining 75% of the total variance. While the distribution of non-immune cells was homogenous, the distribution of immune cells was more scattered and mostly driven by miR-191-5p, miR-24-3p, miR-532-3p, and miR-361-5p (Figure 2A). These four miRNAs were expressed variably across the 90 cell types, but were enriched in immune cells (Supplementary Figure S2A). We therefore performed a second PCA based on expression of these four miRNAs exclusively in circulating PBMCs (CD14 + monocytes, CD8 + T cells, CD4 + T cells, NK cells, and CD19 + B cells). These miRNAs appeared to contribute to the distribution of the cell populations analyzed ( Figure 2B). In particular,, and miR-532-3p were found to be enriched in monocytes, which clustered apart from all lymphoid cells (Supplementary Figure S2B). To reinforce these data, we surveyed an independent dataset (Cohort Roche, GSE28492; Allantaz et al., 2012), which confirmed a higher level of the three miRNAs in monocytes with respect to NK, T, and B cells ( Figure 2C). Of note, miR-191-5p and miR-24-3p were more abundant than miR-532-3p in monocytes. -191-5p, miR-24-3p, and miR-532-3p in monocytes prompted us to search for a correlation with USP18. For this, CD14+ monocytes and monocyte-depleted  -191-5p, miR-24-3p, miR-423-5p, and miR-532-3p in the USP18 3'UTR. Seed-matched sequences are in red. Pairing of miRNA-USP18 sequences was obtained with the RNAhybrid algorithm. (F) HeLa S3 cells were transfected with the indicated miRNAs and 24 h later with the psiCHECK2-USP18 3'UTR reporter plasmid, WT, or mutated in the miRNA binding site (BS mut). Each mutant plasmid contains a 5 nt mutation (nt 2-6 of the seed-matched sequence shown in E). Luciferase activity was measured 24 h after plasmid transfection (n = 4) and expressed as Renilla/ Firefly ratio (± SEM).

The enrichment of miR
Frontiers in Genetics | www.frontiersin.org 6 February 2021 | Volume 11 | Article 627007 PBMCs (here PBL for peripheral blood lymphocytes) were freshly isolated from blood of four donors and levels of the three miRNAs and of USP18 were measured by RT-qPCR. The three miRNAs were more expressed in monocytes than PBL (Figure 3A), while USP18 exhibited an opposite profile ( Figure 3B). In addition, a significant negative correlation between each miRNA and USP18 was observed ( Figure 3C). Interestingly, the opposite was observed for four other ISGs, i.e., OAS1, IRF7, IFIT1, and STAT2, which were found to be more abundant in monocytes than PBL, suggesting a specific profile for USP18 as compared to other ISGs ( Figure 3D).
Altogether, these results suggest that miR-191-5p, miR-24-3p, and miR-532-3p contribute to restraining baseline USP18 in human monocytes. Although not enriched in monocytes as compared to PBL, miR-423-5p may also contribute as it was relatively highly expressed among the studied miRNAs (Supplementary Figure S3).

Copies of USP18 3'UTR Embedded in lincRNA Genes
The results described above suggest that the level of human USP18 can be post-transcriptionally fine-tuned by four miRNAs A B C  are enriched in monocytes. (A) Principal component analysis (PCA) performed on 90 human cell types (FANTOM5 project dataset, three donors per cell population). The distribution is based on expression of the seven microRNAs (miRNAs) significantly targeting ubiquitin-specific peptidase 18 (USP18; see Figure 1B). Red and gray circles: immune and non-immune cells, respectively. The PCA plot shown captures 75% of the total variance within the selected data set (PCA1 35%, PCA2 23%, and PCA3 17%). value of (ANOVA FDR adjusted p-value) q < 0.05. The five USP18-targeting miRNA candidates significantly contributing to the data distribution are shown on the three PC axes of the plot. (B) PCA performed on circulating immune cell populations (FANTOM5 project data set). Analysis was restricted to the four USP18-targeting miRNAs enriched in immune cells. Color legend on the left. The PCA plot captures 97% of the total variance within the selected data set (PCA1 65%, PCA2 28%, and PCA3 6%). value of (ANOVA FDR adjusted p-value) q < 0.05. The four miRNAs significantly contributing to the data distribution are shown on the three principal component axes of the PCA plot. (C) Normalized expression (log2, ±SEM) of miR- 532-3p, miR-191-5p, and miR-24-3p in circulating immune cell subsets. Data retrieved from (Allantaz et al., 2012), Cohort Roche, GSE28492. Ordinary one-way ANOVA (Dunnett's multiple comparison -compared to monocytes) was performed, adjusted p-values shown as stars *p < 0.05, **p < 0.01, ***p < 0.001.
Frontiers in Genetics | www.frontiersin.org 7 February 2021 | Volume 11 | Article 627007 binding its 3'UTR. These miRNAs and their genes are conserved down to the mouse (Supplementary Figure S4A), but this is not the case for their binding sites on USP18 3'UTR. Some of these miRNA-USP18 interactions (as for miR-191-5p and miR-423-5p) may be recent, while others (miR-24-3p and miR-532-3p) may be conserved throughout evolution, as suggested by the conservation of the seed-matched sequences in the USP18 3'UTR of other mammals (Supplementary Figures S4B,C). In searching for sequence conservation, we found several copies of the USP18 3'UTR in the human and chimpanzee genomes ( Table 1), raising the possibility that transcripts other than the bona fide USP18 protein-coding mRNA contain the USP18 3'UTR and bind the same miRNAs.
The USP18 gene and all copies of the 3'UTR sequence map on chr22q11.21, a highly complex region of the human genome that contains four segmental duplications or low-copy repeats (LCR22-A to D; Figure 4A; Morrow et al., 2018). The high sequence identity and copy number variation in the LCR22s have indeed hampered accurate sequencing and gene annotation in this region and, even in the most recent version of the human reference genome, LCR22s are rich in gaps and assembly errors (Demaerel et al., 2019). By performing a BLAT analysis of USP18 gene sequences on the last hg38 assembly, we identified six copies of the 3'UTR in positive or negative strand orientation ( Figure 4A; Table 1). Four copies -named A1-4 here -contain most of intron 10 and the entire exon 11 (i.e., the USP18 3'UTR) of the USP18 gene. The other two copies, named B1 and B2, contain a small segment of intron 10 and the entire exon 11. All six copies reside in LCR-A and LCR-D (Figures 4A,B). Importantly, with the exception of B2, the other copies overlap each an annotated lincRNA gene (see list in Table 2). The  -191-5p, miR-24-3p, and miR-532-3p (A) and of USP18 (B) in ex vivo monocytes and monocyte-depleted peripheral blood mononuclear cells (PBMCs; here PBL; four donors). Results (± SEM) shown as expression (2 -ΔCt ) relative to U6 or 18S, used as housekeeping genes. Paired t-test values of p shown as stars **p < 0.01, ***p < 0.001. (C) Negative correlation between   lincRNA genes spanning copies A1, A3, and A4 encode three annotated transcripts of nearly identical sequence, which were previously described by Delihas (2020b) and named FAM247D, FAM247C, and FAM247A, respectively (HUGO nomenclature; Figure 4C). The annotated transcript spanning A2 (FAM247B) contains only few nucleotides of the 3'UTR and was not further studied. The lincRNA gene spanning B1 (here named linc-UR-B1) encodes two transcript isoforms that are annotated as TCONS_00029753 and TCONS_00029754 (UCSC genome browser) and differ of only 4 nt at the Frontiers in Genetics | www.frontiersin.org 9 February 2021 | Volume 11 | Article 627007 exon-exon junction ( Figure 4D; Supplementary Figure S5A).
As presently annotated, all these transcripts contain only two exons (Figures 4C,D). Next, we measured the expression of these transcripts by qPCR using a commercial panel of cDNAs from polyA+ RNAs of 20 human tissues. In the NCBI database, FAM247A, FAM247C, and FAM247D appear to be widely expressed in human tissues, as measured by RNA sequencing (BioProjects: PRJEB4337, PRJEB2445, PRJNA280600, and PRJNA270632). However, given the repeated nature of their sequences, a precise quantification of their expression requires either the use of specific qPCR primers or additional filters to extract uniquely mapping reads in the RNA-seq data. Thus, we designed specific primers across the unique exon-exon junction of FAM247A, C, and D and across the unique exon-exon junction of linc-UR-B1 (Figures 4C,D; primers in Supplementary Table S2). FAM247 transcripts (FAM247A/C/D) were lowly expressed in most tissues and showed the highest expression in fetal and adult liver, kidney, and thymus ( Figure 4E). On the other hand, linc-UR-B1 was uniquely and abundantly expressed in testis ( Figure 4F). To confirm the testis-restricted expression of linc-UR-B1, we interrogated RNA-seq datasets covering a panel of 32 tissues (E-MTAB-2836). To distinguish linc-UR-B1 from all other USP18containing transcripts, uniquely mapped reads were selected and aligned again to the human genome (hg38). This analysis revealed that linc-UR-B1 is expressed only in testis as a single isoform corresponding to the TCONS_00029754 annotation (Supplementary Figure S5B). This was further confirmed by sequencing the only product amplified by qPCR on testis cDNA (Supplementary Figure S5C). High inter-individual variability in the level of linc-UR-B1 was noticed in the RNA-seq dataset (eight donors; Supplementary Figure S5B). Our qPCR analyses on testis fragments from six donors displaying normal spermatogenesis and one donor with impaired spermatogenesis (as observed by transillumination of seminiferous tubules) confirmed this variability and showed the lowest linc-UR-B1 expression in the latter donor (Supplementary Figure S5D).
As opposed to protein-coding genes, the annotations of lincRNAs are far from being complete (Lagarde et al., 2016). To validate the existing annotations, we performed 3'RACE on RNAs from liver and testis. We confirmed the 3' end of FAM247A/C/D as annotated (Figure 5A, right panel) and showed that the 3' end of linc-UR-B1 extends further than annotated and actually contains the entire USP18 3'UTR (580 nt; Figure 5B, right panel). Mapping the 5' end of these transcripts proved to be challenging, since the sequences that are annotated as exon 1 (i.e., upstream of the 3'UTR) are themselves duplicated, with several copies residing on chromosome 22 and other chromosomes ( Supplementary Table S1). Hence, to enrich for transcripts bearing the USP18 3'UTR, we performed a RT PCR reaction on RNAs from liver and testis using a specific primer (spRT) complementary to the 3'UTR and using forward primers specific to the annotated 5' ends ( Figures 5A,B, left panels). These analyses confirmed the annotated sequences for FAM247A/C/D and linc-UR-B1. Moreover, we could not amplify FAM247A/C/D sequence further using primers designed in the upstream genomic region, suggesting that these RNAs do not extend further.
The genomic region upstream of linc-UR-B1 contains a pseudogene annotated as NR_135922 and belonging to the large POM121 family (Supplementary Figure S6A). Similar to linc-UR-B1, NR_135922 is expressed only in testis (Supplementary Figure S6B). To study whether linc-UR-B1 transcript extends 5' into NR_135922 sequences, we designed primers either in the last annotated exon of NR_135922 or in the intergenic region, and combined them with a reverse primer in the USP18 3'UTR. For all primer pairs, the expected products could be amplified from testis RNA (Figure 5C left panel and Figure 5D), indicating that NR_135922 and linc-UR-B1 are part of the same gene. Indeed, 3'RACE on testis RNA yielded the expected NR_135922 product (537 nt) and a longer isoform (linc-UR-B1) that terminates with USP18 3'UTR ( Figure 5C, right panel and Figure 5D).

Linc-UR-B1 Is Uniquely Expressed in Male Germ Cells and Positively Correlates With USP18
Given the high and specific expression of linc-UR-B1 in the testis, we sought to identify in which testicular cell type(s) this RNA is expressed. We first performed qPCR analysis on cDNAs from freshly isolated testicular germ cell populations (spermatocytes and spermatids) and from primary somatic peritubular, Leydig and Sertoli cells and detected linc-UR-B1 only in spermatocytes and spermatids ( Figure 6A). To obtain further information on linc-UR-B1 expression, we surveyed a  Complementary DNA (cDNA) was obtained by oligo-dT reverse transcription (RT) of polyA+ liver RNA. PCR was performed on this cDNA using the forward (A_FW1) and reverse (AUAP_RV) primers indicated on the schematic above. Several bands were obtained, including the expected 666 bp. Left gel, specific RT-PCR (spRT-PCR). cDNA was obtained from polyA+ liver RNA using a primer designed on the USP18 3'UTR (3UTR_GSP1). PCR was then performed using the forward primer (A_FW1 or A_FW2) and a reverse primer designed on the 3'UTR (3UTR_RV) and internal to 3UTR_GSP1. The length in nt is indicated above each amplified product.
(B) Analysis of linc-UR-B1 transcripts in testis. The strategies described in (A) were used on testis poly A+ RNA. Primers are indicated on the schematic above. The length in nt is indicated above each amplified band. (C) Study of linc-UR-B1 transcript extending on the 5' end. Three forward primers were designed on the annotated intergenic region (B_FW3-4-5). One primer (B_FW6) was designed on the sequence further upstream, corresponding to exon 6 of the annotated NR_135922 transcript. Right panel, 3'RACE on testis RNA. The 537 nt product obtained using the B_FW6 primer corresponds to the annotated NR_135922. Longer products indicated by arrowheads were obtained when using the three primers designed in the intergenic region. Left panel, spRT-PCR was performed on testis RNA using the indicated forward primers and the 3UTR_RV primer. PCR products of the expected lengths were obtained for all primers tested.
Frontiers in Genetics | www.frontiersin.org single-cell RNA-seq dataset of testicular cells from three individuals (Guo et al., 2018). The expression profile of linc-UR-B1 was analyzed by tracking the sequence of the unique exon-exon junction between the 3'UTR and the upstream exon. As shown in Figure 6B, linc-UR-B1 was primarily detected in specific germ cell populations (clusters 4-6). Interestingly, the induction of linc-UR-B1 expression occured in the transition between meiotic early primary and late primary spermatocytes. This burst was followed by a gradual decrease in post-meiotic round and elongated differentiating spermatids. No or low expression was detected in somatic cell populations (clusters 9-13). To identify more precisely the meiotic phase in which linc-UR-B1 induction occurs, we analyzed the dataset of early and late primary spermatocytes reclustered for the five different stages of meiotic prophase I. Linc-UR-B1 was first detectable in the late pachytene stage and was maximal during the diplotene stage ( Figure 6C). Importantly, we observed that the expression of USP18 paralleled that of linc-UR-B1 in germ cell populations and during prophase I (Figures 6D,E). Since USP18 is IFN-inducible, its expression may be related to the presence of local IFN. To test this, we analyzed expression of two ISGs, IFIT1 and OAS1, but no or minimal levels of these transcripts were detected in germ cells (Supplementary Figure S7). In fact, IFIT1 and OAS1 were higher in cells other than germ cells, whereas USP18 was highest in meiotic and post-meiotic germ cells (Figure 6D). The positive correlation between linc-UR-B1 and USP18 expression during spermatogenesis ( Figure 6F) and, in particular, in meiotic cells ( Figure 6G) suggests that these transcripts -which share the same 3'UTRmay cross-talk. Interestingly, we found that testis, and in particular germ cells, express moderate to high levels of three USP18-targeting miRNAs, namely miR-191-5p, miR-24-3p, and miR-423-5p (Supplementary Figure S8). These data hint at the possibility of a miRNA-dependent cross-talk between USP18 and linc-UR-B1 in germ cells, with the latter possibly acting as a sponge for miRNAs targeting USP18.

DISCUSSION
Ubiquitin-specific peptidase 18 determines the threshold of activation of the type I IFN signaling pathway and ultimately the amount of ISGs in a given cell. To exert its function as a negative regulator, USP18 binds to the IFN receptor/JAK complex and attenuates the magnitude of the response (Francois-Newton et al., 2011;Wilmes et al., 2015;Arimoto et al., 2017). The abundance of other IFN-stimulated components (i.e., STAT2/1, IRF9, and SOCS1) will also impact cell contextspecific ISG activation, but in distinct and less specific manners (Kok et al., 2020). Given the critical non-redundant role of USP18, we conceived this study to assess whether USP18 can be post-transcriptionally regulated by miRNAs directed to the 3'UTR. Through the use of target prediction tools, functional validation, public data mining, and correlative analyses, we identified four miRNAs (miR-191-5p, miR-24-3p, miR-423-5p, and miR-532-3p) that directly pair to sites within the 580 nt-long 3'UTR and tune down endogenous USP18 at both mRNA and protein levels. These miRNAs are expressed at different levels in numerous immune and non-immune cell types. We do not know whether all four miRNAs collaborate to robustly target USP18. Yet, are enriched in circulating human monocytes, which we found exhibit low baseline USP18 with respect to PBL.
Monocytes are a subset of leukocytes involved in antimicrobial defenses, anti-tumor immune responses, and in tissue homeostasis (Shi and Pamer, 2011). Interestingly, Uccellini et al. using an ISRE-dependent reporter mouse, showed that Ly6C hi inflammatory monocytes -corresponding to the classical human monocytes and representing ≈ 90% of the circulating ones -display a high basal IFN response, further enhanced upon infection (Uccellini and Garcia-Sastre, 2018). Under physiological conditions, circulating monocytes are primed by constitutive IFNβ, in order to be armed with adequate levels of key proteins (transcription factors, i.e., STATs/IRFs, some  Table S1). (D) Schematic of the 3' end of the two NR_135922 isoforms identified here. Top, the NR_135922 transcript terminating with exon 6, as in the annotation. Bottom, longer isoform of NR_135922 terminating with the USP18 3'UTR, here called linc-UR-B1. The sequences of FAM247A/C/D and linc-UR-B1 are provided as Supplementary Information.  Figure S5A).

All coordinates, genes and transcripts refers to annotations found in the UCSC genome browser.
Frontiers in Genetics | www.frontiersin.org ISGs, and immune response genes products; Gough et al., 2010;Lindqvist et al., 2016). During acute viral infection, in some auto-immune and inflammatory diseases, monocytes encounter higher levels of IFN that in turn promotes their activation and their differentiation into dendritic-like cells with a potent antigenpresenting capacity (Gerlini et al., 2008). Low USP18 may be critical to maintain high responsiveness of these cells to IFN. In line with this, a recent study on ISG15-deficient patients showed that, among PBMCs, monocytes displayed the highest IFN signature (Martin-Fernandez et al., 2020). Thus, the targeting of USP18 by miR-191-5p, miR-24-3p, and miR-532-3p, and possibly miR-423-5p, may assist the priming of monocytes toward a robust inflammatory and immune response.
The conservation analysis of the USP18 3'UTR sequence led us to identify several duplications of the sequences spanning intron 10-exon 11. We found six copies of this sequence mapping in human chr22q11.21. We named these copies as A or B depending on the breakpoint in intron 10 ( Figure 4B). In chr22q11.2 there are eight LCRs, of which four (LCR22A-D) span most of chr22q11.21, a 3 Mb region that is deleted in 90% of patients affected by the 22q11.2 deletion syndrome. First described by DiGeorge in the 1960s, it represents the most frequent chromosomal microdeletion syndrome, affecting 1 per 3,000-6,000 live births. Patients show heterogeneous clinical features, multi-organ dysfunction, cognitive deficits and neuropsychiatric illness, immunodeficiency, and cardiac and palatal abnormalities (Mcdonald-Mcginn et al., 2015;Morrow et al., 2018). Such recurrent genomic rearrangements have fueled many studies to resolve the complex architecture of the LCR22s, identify breakpoints, and study the high variability among patients (Demaerel et al., 2019). The LCR22A-D duplications are believed to have occurred recently, in the hominoid-lineage after their divergence from the Old World monkeys (Babcock et al., 2007;Delihas, 2020a). The highly dynamic nature of LCR22s, as of other LCRs, is believed to be a driving force in genome evolution and adaptation, as their rearrangements can give rise to new coding and non-coding genes (Babcock et al., 2003;Delihas, 2018). Interestingly, the USP18 protein-coding gene resides at the centromeric boundary of LCR22A. Moreover, the gene USP41 resides in the LCR22B ( Figure 4A). In databases, USP41 is annotated as a proteincoding gene. However, this gene contains a sequence nearly identical to USP18, but lacks 5' and 3'UTRs. We could not find experimental evidence that USP41 is transcribed in human tissues. Thus, despite its coding potential, USP41 is likely an unprocessed pseudogene. Previous in-depth studies of the LCR22s did report on several copies of USP18 intron 10-exon 11 and suggested that Alu elements present at their 5' breakpoint (Supplementary Figure S9A) may have driven duplication events (Babcock et al., 2003;Delihas, 2018). Delihas referred to "USP18-linked" sequences as part of a repeat unit present in human and chimpanzee genomes, which includes upstream sequences related to the gamma-glutamyl-transferase (GGT) gene family and downstream sequences related to the FAM230 gene family (Delihas, 2020a). Accordingly, this is the exact genomic context that surrounds the A copies of USP18 exon 11 (Supplementary Figure S9B). Members of the FAM230 family are also found downstream of the two B copies and of the USP18 proteincoding gene, while members of the POM121 family are found upstream (Supplementary Figure S9B). Altogether these observations reinforce the hypothesis that the USP18 intron 10-exon 11 sequence was trapped, together with FAM230 sequences, in a gene-forming element, and this event may have contributed to the generation of new genes (Delihas, 2020a). Interestingly, four of the six copies of USP18 exon 11 are embedded in expressed lincRNA genes, three of which were previously described as FAM247A, C, and D and retain an intact USP18 3'UTR (Delihas, 2020b). Here, we identified an additional lincRNA bearing USP18 exon 11, i.e., the entire 3'UTR, to which we refer as linc-UR-B1 (Figures 4C,D). By tracking the unique junction between USP18 exon 11 and the upstream exon, we showed that FAM247A/C/D are more abundant in tissues like liver, kidney, and thymus, while linc-UR-B1 is uniquely and highly expressed in testis, notably in spermatocytes and spermatids. To our knowledge, this is the first case known of a 3'UTR that is found duplicated independently from the rest of the ancestral gene and embedded in expressed non-coding RNAs. Few non-coding transcripts that bear the 3'UTR sequence of a protein-coding gene have been described. These transcripts are expressed from pseudogenes, notably PTENP1, KRASP1, and BRAFP1. These latter originated by a single duplication of the ancestral protein-coding genes (coding region and UTRs), they are transcribed but not translated due to mutations creating premature stop codons or frameshifts (Poliseno et al., 2010;Karreth et al., 2015). Conversely, the lincRNAs that we have identified bear only the 3'UTR of USP18.
In recent years many non-coding RNAs have been discovered, but only few have been assigned functions (Yao et al., 2019). Moreover, a large number of lncRNA genes are transcriptionally active in testis during meiosis and spermatogenesis (Cabili et al., 2011;Soumillon et al., 2013;Rolland et al., 2019). The lincRNAs we describe here may represent junk DNA, proliferating in a selfish manner. Yet, evolution may have repurposed some of them to confer an advantage, be transmitted, and potentially fixed in the population (Dawkins, 1976;Eddy, 2012). It is a fact that the inclusion of exon 11 -possibly favored by a strong constitutive acceptor site (ttctctctagGCAGGAAACT) at the intron 10-exon 11 junction -provides to these transcripts a canonical AATAAA signal for poly-adenylation. The best studied non-coding transcript bearing the 3'UTR of a proteincoding gene is PTENP1. This pseudogene undergoes copy number loss in human cancers, this correlating with a decrease in levels of the tumor suppressor PTEN. It was proposed that PTENP1 exerts its function by sponging miRNAs, thus sustaining PTEN levels (Poliseno et al., 2010). Likewise, lincRNAs bearing the USP18 3'UTR may act as decoys titrating away USP18targeting miRNAs. Interestingly, our analysis of public RNA-seq data from testicular cell subsets revealed a positive correlation between linc-UR-B1 and USP18 in spermatocytes and spermatids. Moreover, we found that three USP18-targeting miRNAs are expressed in whole testis and in germ cells. Further work is needed to investigate the possibility that in human germ cells linc-UR-B1 and USP18 compete for the binding of these miRNAs.
Of note, the lincRNAs described here contain a small ORF of 43 nt that could encode the last 14 amino acids of USP18 (aa 359-372), eight of which are conserved from zebrafish to primates (Delihas, 2020b). In USP18, these C-ter residues are involved in binding IFNAR2 (Arimoto et al., 2017;Basters et al., 2018). Intriguingly, in linc-UR-B1 an ATG codonpresent at the junction of exon 11 and the upstream exon -may serve as translation initiation codon for the 15 amino acids peptide. The existence of such peptide and its potential impact will have to be investigated in germ cells. Moreover, since linc-UR-B1 is a hybrid of NR_135922 (highly similar to POM121L8P) and USP18 sequences, regions other than exon 11 are likely to contribute to its function(s).
Interestingly, studies in the mouse have shown that high/ continuous IFNα/β affects spermatogenesis by causing apoptosis of germ cells and induces sterility (Hekman et al., 1988;Iwakura et al., 1988;Satie et al., 2011). It is therefore tempting to speculate that the maintenance of baseline USP18 represents a recently evolved mechanism to attenuate IFN signaling and protect the precursors of spermatozoa. In this scenario, human spermatocytes and spermatids may be a good shelter for viruses. Indeed, the testis is a reservoir for a number of viruses, such as Zika and Ebola viruses which can be sexually transmitted by infected men after recovery (Matusali et al., 2018;Le Tortorec et al., 2020). We recently revealed the prolonged seminal excretion of testicular germ cells persistently infected by Zika virus . Future studies will be necessary to determine whether linc-UR-B1 contributes to maintain USP18 in human male germ cells and indirectly favors viral persistence. For this, we will privilege approaches capable of solving RNA-RNA or RNA-protein interactions of the endogenous linc-UR-B1 in the unique context of germ cells.
In conclusion, our work reveals the existence of miRNAs regulating USP18 through its 3'UTR, and of several lincRNAs containing the USP18 3'UTR. Combined, these molecules may form a non-coding network tuning USP18 levels in cell types where IFN responsiveness needs to be tightly controlled.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Ethics Commitee Ouest V, Rennes, France (authorization DC-2016-2783) and the French National Agency for Biomedical Research (authorization PFS09-015). Written informed consent for participation was not required for this study in accordance with the national legislation and the institutional requirements.

AUTHOR CONTRIBUTIONS
ER and MC designed and performed experiments and analyzed the data. ÖV, LV, and GG performed experiments under supervision. NT and ER analyzed the RNA-seq and sequence conservation data. AT, AR, and ND-R provided and performed experiments on testis fragments and purified testicular populations. ND-R and FM provided expert advice on experiments. ER, MC, and SP wrote the manuscript. All authors revised the manuscript. SP supervised the work. All authors contributed to the article and approved the submitted version.

FUNDING
Research in the Unit of Cytokine Signaling is funded by the Institut Pasteur, the Fondation pour la Recherche Médicale (Equipe FRM DEQ20170336741) and the Institut National de la Santé et de la Recherche Médicale (Inserm). ER was supported by Sorbonne Université and by FRM. MC was supported by the FRM grant above. GG and ÖV were supported by the Erasmus plus programme of the EU commission. Work at Rennes University is funded by Inserm and Défis Scientifiques Emergents.