Full-length annotated HLA-B allele sequences were obtained from ImMunoGeneTics (IMGT, version 3.23.0) and aligned using Clustal Omega (12). The multiple sequence alignment was subsequently used to characterize the sequence variability within intron 4 across HLA-B alleles. Sequence logo plots for regions encoding the two mature miRNAs, HSA-miR-6891-5p and HSA-miR-6891-3p, were generated using MATLAB (R2014b) in order to visualize sequence variability within the mature miRNA products. Sequence conservation of the pre-miRNA (HSA-miR-6891) hairpin was determined using BLAST (blastn 2.4.0) against the reference genomic sequence database (refseq_genomic) with the following parameter settings: word size of 28, expected value of 10, hitlist size of 100, match/mismatch scores of 1/−2, gapcosts of 0, low complexity filter on, filter string set to L;m, and genetic code set to 1.

COX cells (13) were obtained from the International Histocompatibility Working Group, Seattle, WA, USA [(IHW09022) http://www.ihwg.org/hla/index.html]. PGF cells (13) were obtained from the Coriell Biorepository (Cat #GM03107). Cells were culured in RPMI-1640 medium with 15% FBS (Sigma Cat #F2442-500ML). HEK 293T cells (a gift from Xianxin Hua in the Department of Cancer Biology at the University of Pennsylvania, Perelman School of Medicine) were cultured in DMEM (Cat #10-013-CV) media with 10% FBS. Primary B-cells were purified from peripheral blood using EasySep™ Direct Human B Cell Isolation Kit from Stemcells Technologies Cat #19674 and directly used for RNA purification. Selective IgA-deficiency patient cell lines (ID 000018, ID000036, ID000038, and ID000057) and cell lines from unaffected, related family members (ID000037 and ID000058) were originally collected and characterized as part of an initiative by the US Immunodeficiency Network and purchased from the Coriell Biorepository.

An “inhibition” lentivirus was generated in cultured HEK 293T cells by transfecting with a pEZX-am03 vector (Genecoepia) containing an HSA-miR-6891-5p antisense insert under the control of a CMV promoter. A lentivirus containing similar insert but with a “scrambled” sequence (i.e., random sequence changes all the bases in the seed region) was similarly generated in HEK 293T cells. Media was discarded after 24 h posttransfection, and packaging media was added to the plate. Scrambled and HSA-miR-6891-5p knockdown viruses were collected every 24 h for 2 days.

For transduction, 1.5 × 10⁵ COX cells were plated in 6-well plates, and 2 mL of fresh scrambled or miR-6891-5p knockdown lentivirus was added along with 4 mg/mL polybrene. The plate was centrifuged at 2,500 rpm for 90 min. After 10 h, 2 mL of additional virus with polybrene was added, and the plate was centrifuged at 2,500 rpm for 90 min. After 16 h, 2 mL of media was discarded, and 2 mL of fresh virus and polybrene were added, and the plate was centrifuged at 2,500 rpm for 90 min. Transduction was allowed to continue for an additional 24 h before cells were collected for RNA extraction. RNA was purified using the miRNeasy kit (Qiagen).

Total RNA extracted from each of the biological replicates of transfected COX cells for both conditions (i.e., inhibition and scrambled control) was used to generate sense-strand cDNA using the Ambion^® WT Expression Kit for Affymetrix^® GeneChip^® Whole Transcript Expression Arrays (P/N 4425209). From each of these reactions, 5.5 μg of sense-strand cDNA was fragmented and labeled using the Affymetrix GeneChip WT Terminal Labeling and Hybridization Kit (PN 702880). Fragmented and labeled sense-strand cDNA (3.25 μg) was hybridized to an Affymetrix Human Gene 2.0ST Array. Arrays were washed on an Affymetrix GeneChip Fluidics Station 450 using fluidics protocol FS450_0002 and scanned on Affymetrix GeneChip Scanner 3000.

Raw data (CEL) files were imported and processed within MATLAB (R2014b). Raw data were first background adjusted using the robust multiarray average procedure, followed by quantile normalization with median polishing and probe level summarization using a custom CDF annotation file (14–18). Those probes on the array with null values for at least one sample were removed so as not to confound subsequent analysis. Principal component analysis (PCA) of the normalized dataset was performed within MATLAB (R2014b) in order to visualize sample clustering and identify sample outliers (19). Differentially expressed transcripts were identified between the miR-6891-5p inhibition samples and control samples (scrambled vector) using significant analysis of microarrays (20). Differentially expressed transcripts were identified using a false discovery rate (FDR) cutoff of 0.05 (21, 22) and a fold change cutoff of 2. Hierarchical clustering of differentially expressed transcripts for each replicate was performed using MATLAB (R2014b). Functional enrichment was performed using DAVID (23). Significant gene ontology (GO) biological processes (BP_FAT) and molecular functions (MF FAT) were determined using a p-value cutoff of 0.05. Raw microarray data are publicly available and may be accessed via NCBI GEO (link provided upon acceptance).

Computationally predicted mRNA targets of HSA-miR-6891-5p were identified throughout the entirety of every annotated gene using miRWalk2.0 with default parameters and every available database, including miRWalk, miRDB, PITA, MicroT4, miRMap, RNA22, miRanda, miRNAMap, RNAhybrid, miRBridge, PICTAR2, and Targetscan (24). The set of genes with a computationally predicted miRNA binding site for miR-6891-5p were then intersected with the set of targets identified by microarray analysis.

Genomic DNA was extracted from the IgA-deficient B-LCLs using the Qiagen Gentra Puregene Blood Kit (Cat No./ID: 158389). Sequencing libraries were generated for each sample using the Omixon Holotype HLA Genotyping Kit as previously described (25). The library was then denatured with NaOH and diluted to a final concentration of 8 pM for optimal cluster density and 600 µL was loaded into the MiSeq reagent cartridge (version 2 500 cycle kit). Samples were demultiplexed on the instrument, and the resulting FASTQ files were used for further analysis. All samples were genotyped at the HLA-B locus using Omixon Target (version 1.8). High-resolution HLA-B genotyping results are found in Table S4 in Supplementary Material.

Total RNA was extracted from cells using a Qiagen miRNeasy kit (Cat #217084) per manufacturer’s protocol (26). Total RNA was reverse transcribed using the Qiagen miRNA Reverse Transcription kit (Cat #218160). qPCR was performed on cDNA generated by reverse transcription using a miSCRIPT SYBR Green PCR kit (Cat #21803). Primers for HSA-miR-6891-5p were obtained from Qiagen (Cat #MS00048202). Primers for IGHA1, IGHA2, β-actin, and HLA-B (sequences provided below) were obtained from IDT. Selective IgA patient cells were HLA genotyped, and qPCR primers were designed to amplify all genotyped HLA-B mRNA transcripts. Selective IgA deficiency cells were harvested, and total RNA was purified using Qiagen miRneasy kit. This RNA was used for qPCR using HLA-B and miR-6891-5p primers. Data were normalized to actin. Significance was assessed using an unpaired one-tailed t-test.

Primer sequences: (i)

IGHA1

COX and PGF cells were cultured in RPMI-1640 media. After 72 and/or 120 h, media were collected, and IgA secretion was analyzed using Ready-SET-Go ELISA kit (27, 28) (Cat #88-50600) from Affymetrix (CA) per manufacturer’s protocol. Significance was assessed using the one-tailed t-test.

The complete (48 nucleotide) 3′UTR sequence of the IGHA1 gene (which is identical to the 3′UTR sequence of the IGHA2 gene), containing the HSA-miR-6891-5p binding site, was synthesized with PmeI and XbaI sites on either end (IDT) and gel purified using a QIAquick Gel Extraction Kit (Qiagen Cat #28704). The product was ligated into the pmiRGLO plasmid (Promega, WI, USA) digested with PmeI and XbaI (New England Biolabs, MA, USA) downstream of the PGK promoter and luciferase gene.

For the luciferase assay, 1 × 10⁶ HEK 293T cells were cultured in multiwell plates and, after 24 h, were transfected with either the wild-type IGHA1 3′UTR or mutant IGHA1 3′UTR construct using Fugene 6 (Promega, Cat #E2691). Some of these cells were also transfected with either HSA-miR-6891-5p antisense or overexpression constructs. After 24 h, the cells were assayed for luciferase activity using the Dual-Luciferase^® Reporter Assay System (Promega, Cat #E1910) (29). For each measurement, firefly luciferase data were normalized to renilla luciferase. Significance was assessed using Student’s t-test.

Following transcript splicing, intron 4 of HLA-B is predicted to form a pre-miRNA hairpin that is further processed by the Dicer enzyme into two mature miRNA products, miR-6891-5p and miR-6891-3p (7) (Figure 1). Given the highly polymorphic nature of the HLA-B locus, we explored miR-6891 sequence variants (isomiRs) by interrogating the sequences of intron 4 among the 384 full-length annotated HLA-B alleles in the international IMGT database (IMGT/HLA, release 3.25) (30). Among those, only eight unique sequence motifs were observed (Figure 2A). Remarkably and despite the very polymorphic nature of the HLA-B gene, there is no sequence variation within miR-6891-5p (Figure 2B) and only two polymorphic sites within the mature miR-6891-3p arm, occurring at positions 6 and 14 of the mature miRNA (Figure 2C). Each of these intronic sequences form stable pre-miRNA hairpin structures with secondary structure minimum free energy values ranging from −43 to −54 kcal/mol. We selected miR-6891-5p for additional study because its conserved sequence suggests an important biological role. The pre-miRNA hairpin sequence of hsa-miR-6891 is evolutionarily conserved, with 90% sequence identity among six primate species including Homo sapiens, Gorilla gorilla, Nomascus leucogenys, Chlorocebus sabaeus, Macaca nemestrina, and Macaca mulatta. In contrast, the closest homolog of hsa-miR-6891 within the mouse genome, which lies within intron 5 of the H2-T10 gene, has only 48% (45/93 base positions identical) sequence conservation with hsa-miR-6891, and there is no annotated miRNA encoded within this locus (miRbase release 21).

To study the function of miR-6891-5p, an appropriate in vitro cell model was first identified by examining the expression level of miR-6891-5p within two B-LCLs, PGF and COX, as well as immortalized HEK293T cells and primary B-lymphocytes (Figure S1 in Supplementary Material). Our qPCR results indicate that miR-6891-5p is expressed in every cell type analyzed, with B-LCLs exhibiting the highest and most uniform expression of miR-6891-5p across biological replicates. For this reason, B-LCLs (COX cells) were selected as our model system to further study the role of miR-6891-5p.

To identify putative target transcripts of miR-6891-5p, we transduced COX cells with a lentiviral construct expressing the antisense transcript of miR-6891-5p to inhibit the activity of miR-6891-5p. Upon transduction, transcripts targeted by miR-6891-5p are expected to be more abundant within the transduced cells since miR-dependent degradation has been inhibited. The experimental design included COX cells transduced with either the lentiviral construct expressing the antisense sequence (inhibition) or scrambled antisense sequence (control) of miR-6891-5p. Adequate and comparable expression of lentiviral constructs from both experimental conditions was observed (Figures S2 and S3 in Supplementary Material). Affymetrix Human Gene 2.0ST Arrays were used to assess transcript expression levels between the miR-6891-5p inhibition and control sample groups. PCA (Figure 3A) of the normalized microarray data demonstrates excellent clustering of the two distinct cell populations, indicating distinct and reproducible mRNA expression profiles among biological replicates. Transcripts with significant differential expression between the miR-6891-5p inhibition and control sample groups were identified. One hundred four upregulated and 99 downregulated transcripts were observed within the miR-6891-5p inhibition sample group when compared with the control group, using a fold change cutoff of ≥2 and a FDR cutoff of 0.05 (Figure 3B; Tables S1 and S2 in Supplementary Material, respectively). Since miRNA are known to bind and downregulate the expression of targeted mRNA transcripts, only those transcripts that were identified as upregulated in the miR-6891-5p inhibition sample group were considered to be putative direct targets of miR-6891-5p (Table S1 in Supplementary Material), whereas the set of downregulated transcripts may be related to indirect effects of miR-6891-5p inhibition (Table S2 in Supplementary Material). The potential binding sites of miR-6891-5p within the 104 upregulated transcripts were identified using an in silico miRNA target prediction algorithm. Among the 104 empirically identified putative targets of miR-6891-5p, 61 (~58%) were found to harbor a computationally predicted miRNA binding site for miR-6891-5p (Table S1 in Supplementary Material).

Functional analysis of differentially expressed transcripts was performed by determining the enriched GO biological processes of upregulated and downregulated transcripts (Table S3 in Supplementary Material). Significantly upregulated transcripts were found to be involved in numerous immunological processes including leukocyte and mast cell activation (GIMAP5, EGR1, NDRG1, and LCP2) and various cellular processes including T-cell antigen receptor-mediated signaling (LCP2) and T-cell quiescence (GIMAP5). Also, the significant upregulated transcripts are 11 DNA-binding proteins and transcription factors (FOS, EGR1, LEF1, TP63, HIST1H2AG, ZFHX4, ZNF730, and ZNF83) including the transcriptional repressor genes SNAI2, PCGF2, and ZNF253. Significantly downregulated transcripts, following miR-6891-5p inhibition, are involved in numerous immunological processes including cytokine production (FCER1G, HMOX1, IFNG, IL10, NFATC2, SIRT1, and TSPAN6), regulation of B-cell-mediated immunity (FCER1G, IFNG, and IL10), inflammation (CCR1, CXCL10, FCER1G, HMOX1, IL10, PNMA1, and PPARG), and the immunoglobulin-mediated immune response (FCER1G, IFNG, and IL10).

The IgA heavy chain encoding transcript was among the most significantly upregulated transcripts following inhibition of miR-6891-5p identified from the microarray analysis (8.5-fold change, FDR = 0.02). To study the role of miR-6891-5p on the abundance of both the IgA mRNA transcript and secreted IgA protein, IgA secreting COX cells were transduced with a lentivirus expressing either the antisense of miR-6891-5p (miR-6891-5p inhibition) or a scrambled antisense sequence of miR-6891-5p (control). Inhibition of miR-6891-5p within COX cells significantly increased the abundance of both the IGHA1 and IGHA2 mRNA transcripts (p = 0.028 and p = 0.007, respectively) (Figure 4A) and secreted IgA protein (p = 0.033) (Figure 4B) compared to cells transduced with the scrambled control. These findings demonstrate that miR-6891-5p inhibits the expression of both IGHA1 and IGHA2. Our in silico molecular modeling of both IGHA1 and IGHA2 transcripts reveals an energetically favorable binding site of miR-6891-5p on the 3′UTR of IGHA1 that is 100% conserved within the 3′UTR of the IGHA2 transcript, suggesting that miR-6891-5p may bind and regulate the expression of both transcripts. The identified non-canonical heteroduplex contains limited base pairing between the miRNA seed region (positions 2–7 of the 5′ end) and the conserved 3′UTR sequence of the IGHA1 and IGHA2 transcripts (Figure 4C).

To validate functional targeting of the modeled miR-6891-5p-binding site within the conserved 3′UTR sequences of both IGHA1 and IGHA2, the UTR sequence was fused to a plasmid-based luciferase reporter and transfected into HEK293T cells. HEK293T cells express miR-6891-5p but not IgA and thus provide a cell model system to study IGHA 3′UTR targeting without competitive binding from endogenously expressed IGHA mRNA. These cells were also transfected with either the miR-6891-5p antisense expression plasmid to inhibit endogenously expressed miR-6891-5p (miR inhibition) or a plasmid expressing miR-6891-5p to increase the level of the endogenously expressed miR-6891-5p (miR overexpression). Inhibition of miR-6891-5p significantly increased luciferase activity (p = 0.013), whereas overexpression of miR-6891-5p significantly attenuated luciferase activity (p = 0.026). Further validation of the binding site was performed by mutating the 3′UTR sequence underlying the binding site of miR-6891-5p (Figure 4D) and fusing it to a plasmid-based luciferase reporter, which was then transfected into HEK293T cells. These cells were also transfected with either the miR-6891-5p antisense expression plasmid to inhibit endogenously expressed miR-6891-5p or a plasmid expressing miR-6891-5p to increase the level of the endogenously expressed miR-6891-5p. In contrast to the wild-type 3′UTR luciferase experiments, no modulation of miR-6891-5p (inhibition or overexpression) was able to affect luciferase activity (Figure 4E), indicating that miR-6891-5p was unable to bind the mutant 3′UTR sequence. Together, these results suggest direct miR-6891-5p targeting on the 3′UTR of both the IGHA1 and IGHA2 transcripts.

Given our findings, we investigated the putative role of miR-6891-5p on the expression and secretion of IgA within B-LCLs obtained from two familial cohorts, consisting of individuals affected by selective IgA deficiency and unaffected relatives (Figure 5A). In order to design effective qPCR primers that amplify the HLA-B mRNA transcripts of each individual, high-resolution HLA genotyping was performed on all affected and unaffected individuals for eight HLA loci (Table S4 in Supplementary Material). Phased MHC haplotypes were inferred from related individuals using the family pedigree when available (ID57, ID58, ID38, ID37, and ID36) or from common MHC haplotypes otherwise (ID18). Expression of HLA-B, miR-6891-5p, IGHA1, and IGHA2 was quantified by sequence-specific qPCR primers (Figure 5B). We observe that IGHA1 is the primarily expressed heavy chain transcript of IgA across all individuals and demonstrate an inverse correlation between miR-6891-5p expression and IGHA1 expression (Pearson correlation −0.87), as well as a strong correlation between HLA-B and miR-6891-5p expression (Pearson correlation 0.96), across all patient samples. Both families showed increased expression of both HLA-B (ID36/ID37 = 18.6×; ID38/ID58 = 4.2×; ID57/ID58 = 4.9×) and miR-6891-5p (ID36/ID37 = 5.3×; ID38/ID58 = 3.5×; ID57/ID58 = 16.8×). In all cases, miR-6891-5p expression was found to be less than that of the host gene, HLA-B. Additionally, inhibition of miR-6891-5p, in an IgA-deficient cell line (ID18), led to a significant increase (~3×) in both IGHA1 and IGHA2 transcript abundance (p = 0.006 and p = 0.043, respectively) as well as a significant increase in the concentration of secreted IgA protein (p = 0.004) (Figure 5C).