Genome-Wide Analysis of Differentially Expressed miRNAs and Their Associated Regulatory Networks in Lenses Deficient for the Congenital Cataract-Linked Tudor Domain Containing Protein TDRD7

Mutations/deficiency of TDRD7, encoding a tudor domain protein involved in post-transcriptional gene expression control, causes early onset cataract in humans. While Tdrd7 is implicated in the control of key lens mRNAs, the impact of Tdrd7 deficiency on microRNAs (miRNAs) and how this contributes to transcriptome misexpression and to cataracts, is undefined. We address this critical knowledge-gap by investigating Tdrd7-targeted knockout (Tdrd7-/-) mice that exhibit fully penetrant juvenile cataracts. We performed Affymetrix miRNA 3.0 microarray analysis on Tdrd7-/- mouse lenses at postnatal day (P) 4, a stage preceding cataract formation. This analysis identifies 22 miRNAs [14 over-expressed (miR-15a, miR-19a, miR-138, miR-328, miR-339, miR-345, miR-378b, miR-384, miR-467a, miR-1224, miR-1935, miR-1946a, miR-3102, miR-3107), 8 reduced (let-7b, miR-34c, miR-298, miR-382, miR-409, miR-1198, miR-1947, miR-3092)] to be significantly misexpressed (fold-change ≥ ± 1.2, p-value < 0.05) in Tdrd7-/- lenses. To understand how these misexpressed miRNAs impact Tdrd7-/- cataract, we predicted their mRNA targets and examined their misexpression upon Tdrd7-deficiency by performing comparative transcriptomics analysis on P4 and P30 Tdrd7-/- lens. To prioritize these target mRNAs, we used various stringency filters (e.g., fold-change in Tdrd7-/- lens, iSyTE-based lens-enriched expression) and identified 98 reduced and 89 elevated mRNA targets for overexpressed and reduced miRNAs, respectively, which were classified as “top-priority” “high-priority,” and “promising” candidates. For Tdrd7-/- lens overexpressed miRNAs, this approach identified 18 top-priority reduced target mRNAs: Alad, Ankrd46, Ceacam10, Dgat2, Ednrb, H2-Eb1, Klhl22, Lin7a, Loxl1, Lpin1, Npc1, Olfm1, Ppm1e, Ppp1r1a, Rgs8, Shisa4, Snx22 and Wnk2. Majority of these targets were also altered in other gene-specific perturbation mouse models (e.g., Brg1, E2f1/E2f2/E2f3, Foxe3, Hsf4, Klf4, Mafg/Mafk, Notch) of lens defects/cataract, suggesting their importance to lens biology. Gene ontology (GO) provided further insight into their relevance to lens pathology. For example, the Tdrd7-deficient lens capsule defect may be explained by reduced mRNA targets (e.g., Col4a3, Loxl1, Timp2, Timp3) associated with “basement membrane”. GO analysis also identified new genes (e.g., Casz1, Rasgrp1) recently linked to lens biology/pathology. Together, these analyses define a new Tdrd7-downstream miRNA-mRNA network, in turn, uncovering several new mRNA targets and their associated pathways relevant to lens biology and offering molecular insights into the pathology of congenital cataract.


INTRODUCTION
Perturbations in lens development results in congenital cataract in humans and animal models (Graw, 2009;Shiels and Hejtmancik, 2019). Studies over the past several decades have led to a detailed understanding of the key signaling and transcriptional regulatory mechanisms that orchestrate the genetic program of lens development (Donner et al., 2006;Lachke and Maas, 2010;Cvekl and Zhang, 2017). However, compared to signaling and transcription, the impact of post-transcriptional gene expression control to organogenesis, in general, and lens development, in particular, remains relatively understudied (Lachke and Maas, 2011;Blackinton and Keene, 2014;Dash et al., 2016;Cvekl and Zhang, 2017). Post-transcriptional control of gene expression is defined as the regulation of any of the different events from the processing of pre-mRNA to the degradation of mRNA (Singh et al., 2015). Indeed, non-coding RNAs such as microRNAs (miRNAs) as well as RNA-binding proteins (RBPs) are involved in various post-transcriptional regulatory processes, including control over translation or decay of mRNA (Pasquinelli, 2012;Manning and Cooper, 2017;Hentze et al., 2018;O'Brien et al., 2018). Thus far, very few RBPs and posttranscriptional regulatory factors, including miRNAs, have been functionally implicated in lens development and cataract (Lorén et al., 2009;Choudhuri et al., 2013;Shaham et al., 2013;Wolf et al., 2013;Xie et al., 2014;Dash et al., 2015Dash et al., , 2020Siddam et al., 2018;Aryal et al., 2020;Barnum et al., 2020;Nakazawa et al., 2020;Shao et al., 2020). This limited information highlights a substantial knowledge-gap in lens and cataract research because post-transcriptional control represents critical mechanisms that allow precise calibration, in terms of dosage and spatiotemporal pattern, of the cellular proteome. Thus, these regulatory mechanisms may be significant for controlling mRNA and protein abundance in lens fiber cells -a cell fate that faces added challenges to regulate these basic processes as they undergo nuclear degradation in terminal differentiation (Dash et al., 2016).
In embryogenesis, expression of Tdrd7 (Tudor domain containing 7) is highly enriched in lens fiber cells and is conserved between aves and mammals, suggesting its critical function in lens development . Indeed, TDRD7 mutations or deficiency results in congenital cataract in humans, mouse and chicken Tanaka et al., 2011;Chen et al., 2017;Tan et al., 2019). Furthermore, single nucleotide polymorphisms in TDRD7 are also linked to age-related cataract (Zheng et al., 2014), making it among the select few genes associated with both early-and late-onset cataract (Shiels and Hejtmancik, 2019). Tdrd7 contains three tudor domains and three OST-HTH (Oskar-Tdrd7-Helix turn helix)/LOTUS motifs (Hosokawa et al., 2007;Anantharaman et al., 2010;Callebaut and Mornon, 2010;Tanaka et al., 2011). Tudor domains are considered to facilitate interaction with methylated arginine or lysine residues within other proteins (Chen et al., 2011;Pek et al., 2012;Gan et al., 2019). The OST-HTH/LOTUS domains in Drosophila protein oskar, predicted to bind RNA (Anantharaman et al., 2010;Callebaut and Mornon, 2010), has been shown to interact with a dead-box helicase (Jeske et al., 2015(Jeske et al., , 2017. While Tdrd-family proteins have been implicated in the control of small RNAs, previous studies have primarily focused on their association with piwi-interacting RNAs (piRNAs) in the context of spermatogenesis (Pek et al., 2012;Gan et al., 2019). Interestingly, while its mutation or deficiency is linked to azoospermia in human and mouse Tanaka et al., 2011;Tan et al., 2019), Tdrd7 has been shown to function in repression of Line-1 retrotransposons but is not found to be essential for production of piRNAs in mouse spermatogenesis (Tanaka et al., 2011). However, the effect of Tdrd7 deletion on other classes of small non-coding RNAs such as miRNAs in the lens, and its impact on lens development and cataract formation has not been addressed.
Previously, we used Tdrd7-deficient mice-which exhibit fully penetrant cataracts and re-capitulate features of the human lens defects-to gain insight into Tdrd7's role in lens development Barnum et al., 2020). We showed that removal of Tdrd7 results in misexpression of several lens expressed mRNAs. We also demonstrated that Tdrd7 protein closely associates with specific mRNAs, for example Hspb1 mRNA, which may enable it to directly control its abundance in the lens. Here, we sought to examine the impact of Tdrd7 deletion on global miRNA abundance in the lens. We performed microarray-based profiling of miRNAs in Tdrd7-/mouse lens at postnatal day (P) 4, which precedes detectable lens defects and overt cataract formation. We identified misexpressed miRNAs and predicted their mRNA targets in the lens. We then performed comparative analysis with Tdrd7-/-lens RNAseq and microarray datasets to identify inversely associated mRNA targets of the misexpressed miRNAs. This allowed us to derive a Tdrd7-downstream miRNA-mRNA network that led to the identification of new candidate genes with potential function in the lens. Importantly, iSyTE analysis showed that similar numbers of overexpressed and reduced mRNA targets in Tdrd7-/-lens were enriched in normal lens, suggesting that Tdrd7 functions in optimal control of a subset of lens-enriched mRNAs likely via regulating miRNAs. Finally, in addition to informing on Tdrd7 deficiency, the global lens miRNA profile generated in this study provides independent support for expression of some of the abundant miRNAs in normal lens development that were described in previous studies, involving for example, in situ hybridization assays (Conte et al., 2010;Karali et al., 2010;Khan et al., 2016). These are: miR-184, miR-26a, let-7b, let-7c, miR-204 and miR-125b, among others. Together these data identify new high-priority Tdrd7-downstream miRNA and mRNA targets, thereby advancing our understanding of how this conserved Tudor family protein functions to fine-tune lens transcriptome in development and how its misregulation impacts cataract pathology.

Mouse Studies
Mice were maintained at the University of Delaware Animal Facility. Experimental protocols followed guidelines based on the Association for Research in Vision and Ophthalmology (ARVO) Statement for the use of animals in ophthalmic and vision research and were approved by the University of Delaware Institutional Animal Care and Use Committee (IACUC). The present studies were performed on Tdrd7 targeted germline knockout (KO) mouse line (Tdrd7 TM1.1Chum , hereafter referred as Tdrd7-/-) that were genotyped as previously described (Tanaka et al., 2011).

Sample Preparation and miRNA Microarray Analyses
For miRNA microarray, microdissected mouse lenses at postnatal day (P) 4 were collected from Tdrd7-/-and control (Tdrd7±, which does not develop cataract) in three biological replicates. Total RNA isolation was performed using the mirVana TM RNA isolation kit (Life Technologies, Grand Island, NY). Global expression profiling for miRNAs was performed using Affymetrix miRNA 3.0 arrays. Analysis of the raw expression datasets was performed under 'R' Statistical environment [http://www.r-project.org/index.html] using "Affy" packages. The datasets were background corrected, normalized and summarized using Robust Multi-array Average (RMA) method (Irizarry et al., 2003a,b). The obtained normalized miRNA expression values were subjected to downstream analysis using 'limma' package. Comparisons of control and Tdrd7-/samples was carried out to identify highly and differentially expressed miRNAs with significant p-value (< 0.05) and absolute fold change (FC) ≥ ± 1.2. The detailed pipeline for Affymetrix microarray dataset analysis after data normalization is published elsewhere Kakrana et al., 2018). The datasets generated in this study were deposited in GSE157061.

qPCR Analysis of miRNAs in the Lens
Selected Tdrd7-/-lens misexpressed candidate miRNAs were analyzed by custom-designed PCR primers using miRCURY LNA miRNA PCR system (Qiagen, Germantown, MD). Total RNA was isolated from P4 Tdrd7-/-and control lenses in three biological replicates using miRNeasy mini kit (Qiagen Catalog: 217004). First-strand cDNA synthesis was performed using miRCURY LNA Universal RT kit (Qiagen Catalog: 339340) and miRNA expression was quantified using miRCURY LNA SYBR Green PCR kit (Qiagen Catalog: 339345) according to the manufacturer's instructions. qPCR was run on BioRad CFX RT-PCR thermal cycler. MiRNA expression was normalized to miR-17, which exhibits robust expression in the lens and is not altered in Tdrd7-/-lens, as well as the housekeeping genes Gapdh and Actb. Relative expression was estimated using the 2 −− CT followed by statistical two-level nested analysis of variance test to calculate p-values.

miRNA Target Prediction and Gene Ontology Enrichment Analysis of miRNA Targets
We performed miRNA target prediction and functional annotation analysis for differentially expressed miRNAs in Tdrd7-/-lenses. To retrieve predicted mRNA targets of miRNAs, the miRDB resource (Chen and Wang, 2020) that uses experimental data (e.g., miRNA overexpression, publicly available CLIP-seq data) in their miRNA-mRNA target prediction algorithm, MirTarget, was used. The miRDB-identified mRNA targets with a probability score ≥ 50 were retrieved for further downstream analysis and were tested for their functional relevance to lens development using gene ontology (GO) analysis in DAVID bioinformatics tool (https://david.ncifcrf.gov). The top clusters with highest enrichment scores with functionally relevant GO categories at p-value ≤ 0.05 were presented. The TargetScan database (www.targetscan.org) was used to search for target mRNA conservation sites that match with the seed region of miRNAs across vertebrates (Lewis et al., 2005).

Identification of Lens Enriched miRNA-mRNA Pairs
The miRNA target mRNAs were curated using lens genediscovery bioinformatics tool-iSyTE (Lachke et al., 2012;Kakrana et al., 2018) to identify candidate miRNA-mRNA pairs (i.e., miRNA and their predicted mRNA targets) associated with lens function and development. In this study, miRNA-mRNA pairs were tested for lens expression (expression ≥ 100) from P4 mRNA microarray dataset in iSyTE . To further assess candidate miRNA-mRNA pairs, we revisited the whole genome RNA-seq (GSE134384) and mRNA microarray datasets on Tdrd7 null and control lenses at P4 (GSE25775)  and performed an analysis to evaluated inverse correlation between differentially expressed 'up' and 'down' miRNA with their mRNA targets. The inverse correlation for elevated and reduced (used in the context of miRNAs) miRNA-mRNA pairs was tested using Pearson method implemented in R-package. Data was presented as miRNA-mRNA plots and miRNA-mRNA regulatory network was visualized as regulatory networks using cytoscape (https:// cytoscape.org).

RESULTS AND DISCUSSION
Global miRNA Profiling of Tdrd7-/-Lens Prior to Formation of Cataract Tdrd7 expression in mouse embryonic lens fiber cells is highly enriched on both the transcript and protein levels (Figures 1A-C) and its germline knockout results in severe cataract defects by postnatal day (P) 22 ( Figure 1D). While the mRNAs misexpressed in Tdrd7-/-mouse lens have been examined on the genomic level in previous studies Barnum et al., 2020), there is no information on the impact of Tdrd7 deficiency on lens miRNA expression. Therefore, we sought to examine Tdrd7-/-lens by global miRNA profiling (Figure 2A). Data on three biological replicates each for Tdrd7-/-and control lenses at postnatal day (P) 4 were generated on Affymetrix miRNA 3.0 microarrays. The stage P4 was selected for this analysis because it preceded cataract formation in Tdrd7-/-lenses ( Figure 1D) and therefore was considered to yield information on the early alterations in miRNA expression, prior to the onset of lens defects, in turn reducing the possibility of detecting secondary miRNA changes. The microarray datasets generated were imported in R-statistical environment for systematic analysis of miRNAs expression in the Tdrd7-/-and control lens samples. The detailed bioinformatics pipeline is outlined in Figure 2A. The quality of miRNA microarray datasets was evaluated using Principal component analysis (PCA) ( Figure 2B) and boxplot analysis ( Figure 2C). Normalized miRNA expression intensity in PCA showed that Tdrd7-/-sample datasets (T1, T2, T3) were distinct from control FIGURE 1 | Tdrd7 is expressed in mouse embryonic lens development and its deficiency results in early onset cataract. (A) In situ hybridization of mouse embryonic day stage E12.5 with Tdrd7-specific probe shows Tdrd7 transcript signal in the eye lens (asterisk), (A') also clearly observed at high-magnification (asterisk). (B) Section in situ hybridization confirms Tdrd7 transcripts to be expressed in the mouse E12.5 lens (le) and demonstrates it to be primarily localized to fiber cells. (C) Immunofluorescence staining using antibody against Tdrd7 demonstrates its protein abundance (red) in the cytoplasm of lens fiber cells at E12.5. (D) Tdrd7-/-germline targeted knockout mouse lens show severe cataract (asterisk) compared to normal control at P22 stage. le, lens; r, retina; e, lens epithelium; fc, lens fiber cells.

Prediction of Downstream mRNA Targets of Misexpressed miRNAs in Tdrd7-/-Lens
MicroRNAs control the cellular proteome by directly binding to their target mRNAs and channeling them to degradation or by inhibiting their translation into protein (Pasquinelli, 2012;Manning and Cooper, 2017;Hentze et al., 2018). In both cases, direct binding of miRNA to their target mRNA is a critical step. Therefore, to gain insight into the downstream impact of the differentially expressed miRNAs in the Tdrd7-/lens transcriptome, we sought to identify miRNA-mRNA target binding pairs using the miRDB resource (Chen and Wang, 2020). This identified potential 11450 and 7142 mRNA targets for 12 elevated and 7 reduced miRNAs, respectively. The database did not contain information on mRNA targets for the differentially expressed miRNAs miR-3107, miR-1935 and miR-3092, and therefore these miRNAs were not included in further downstream analyses. Next, from the numerous potential target mRNAs, we sought to prioritize key candidates. Therefore, we first analyzed these target mRNAs in the context of genomewide RNA-seq data on P4 Tdrd7-/-lenses (Barnum et al., 2020) and normal lens enriched-expression at P4 using the iSyTE database . At P4, RNA-seq showed Tdrd7-/lens exhibited mis-expression of 1982 reduced and 1832 elevated mRNAs ( Table 2). Further, among the 1982 mRNAs reduced in the Tdrd7-/-lens, 574 unique mRNAs (∼29%) were found to be direct targets of 12 miRNAs elevated in the Tdrd7-/-lens (hereafter referred to as miRNA-mRNA pairs). Similarly, among the 1832 elevated mRNAs in the Tdrd7-/-lens, 535 (∼29%) mRNAs were found to be direct targets of 7 miRNAs reduced in the Tdrd7-/-lens (Table 2). Thus, 1109 mRNAs collectively misexpressed and inversely correlated with 19 miRNAs in Tdrd7-/-P4 lenses. Next, the inversely correlated miRNA-regulated mRNAs were tested for lens enriched expression at stage P4 lens using the iSyTE database. We found that the target mRNAs in the 87 miRNA-mRNA pairs (representing 54 unique mRNAs; the reduced number for mRNA targets is due to multiple miRNAs sharing target mRNAs) exhibited lens-enriched expression for miRNAs elevated in the Tdrd7-/-lens ( Figure 3B). The target mRNAs in the 81 miRNA-mRNA pairs (representing 59 unique mRNAs) exhibited lens-enriched expression for miRNAs reduced in the Tdrd7-/-lens ( Figure 3B). These data showed that similar proportions of elevated or reduced miRNA target mRNAs exhibited lens-enriched expression in normal development. We extended this analysis by applying a set of distinct filters and by including lens expression data at a later stage (i.e., P30). This strategy applied differential expression data in Tdrd7-/lens at both P4 and P30 and/or normal lens expression at P4 and P30 to identify miRNA-mRNA pairs as either "toppriority", "high-priority" or "promising" candidates ( Figure 3C).  unique target mRNAs reduced in Tdrd7-/-lens, also termed as belonging to "elevated miRNA group") and 123 reduced miRNA-mRNA pairs (89 unique target mRNAs elevated in Tdrd7-/lens, also termed as belonging to "reduced miRNA group") was validated by performing Pearson correlation analysis. The correlation coefficient (r) for elevated and reduced miRNAs with their target mRNAs was -0.975 at a significant p-value of 0.00001 ( Figure 3D). The "top-priority", "high-priority" and "promising" mRNAs that are targets of elevated and reduced miRNAswhich are also identified among differentially expressed mRNAs in Tdrd7-/-lenses -are presented as heat map in the context of normal lens gene expression (Figures 4, 5).
We next examined the conservation of miRNA seed type in target mRNAs across vertebrates using the TargetScan database. In the elevated miRNA group, the following miRNA target mRNAs were conserved across vertebrates: miR-15a target mRNAs-Egln2 and Aatk; miR-19a target mRNAs-Epn2, Plekhg5 and Timp2; miR-138 target mRNA Nfix. While Egln2 and Plekhg5 have a 7mer conserved seed type, Aatk, Epn2, Timp2 and Nfix have an 8mer conserved seed type. In the reduced miRNA group, the following target mRNAs were conserved across vertebrates: let-7b target mRNAs Sh3glb1, Elovl4, Bcap29; and miR-34c target mRNA Tom1 were conserved at 7mer seed type, while miR-34c target mRNA Tpd52 was conserved at 8mer seed type.

Functional Insights Into mRNA Targets of miRNAs in Tdrd7-/-Lens
The prioritized candidates from the above analysis represent numerous miRNA-target mRNA genes whose functions are relevant to the observed lens defects in Tdrd7-/-mice. To systematically associate cellular function we imported these unique mRNAs (98 from elevated and 89 from reduced miRNA group, please see above) for analysis using the DAVID bioinformatics tool, and identified functional gene clusters (Figure 6 and Supplementary Tables 2, 3). Notably, several significantly enriched gene ontology (GO) categories (enrichment score > 1.0, p < 0.05) for the elevated and reduced miRNA groups are relevant to Tdrd7 function in the lens (Figures 6A,B). For example, for elevated miRNA group (i.e., reduced mRNAs in Tdrd7-/-lens), GO:0042995 "cell projection" contained ten genes. These candidates were Rgs8, Olfm1 (top-priority), Rilpl2, Ncdn, Anks1 (high-priority) and Zfp385a, Plekhg5, Ctnnd2, Mylk, Aatk (promising category). Further, GO:0005604 "basement membrane" contained four genes including Loxl1 (top-priority), Col4a3, Timp2, Timp3 (promising category). This suggests Tdrd7 functions in downstream control of genes involved in cell-cell interaction or connectivity, in turn suggesting their potential impact on Tdrd7-/-cataract pathology, which is associated with vacuolelike gaps in fiber cells and posterior capsule rupture Tanaka et al., 2011;Barnum et al., 2020). Other interesting top clusters with GO categories for elevated miRNA group were GO:0035556 "intracellular signal transduction" (with promising candidates Gna12, Rasgrp1, Rassf5, Plekhm1, Prkcd), UP_keyword "metal binding" (with top-priority candidates Alad, FIGURE 4 | Continued FIGURE 4 | "Top-priority" and "High-priority" mRNA targets of Tdrd7-downstream miRNAs in the lens. Heat-maps representing expression data of "Top-priority" and "High-priority" mRNA targets of miRNAs identified in (A) "Elevated miRNA group" and (B) "Reduced miRNA group" in Tdrd7-/-lens. Differential expression in P4 (RNA-seq and microarray data) and P30 (microarray data) in Tdrd7-/-lens (left) and enriched expression and expression in normal lens as per iSyTE data at postnatal day (P) 4 and P30 (right). "iSyTE Enrichment" represents lens-enriched expression in fold change, while "iSyTE Expression" represents lens expression in fluorescence intensity units at P4 and P30. The mRNA targets are grouped based on the filtering criteria outlined in the text and in Figure 3C. Heatmap keys indicate the following: Yellow and purple color gradient represent low to high differential expression in Tdrd7-/-compared to control lens (in Fold-change) in RNA-seq data and microarray data. Green and red color gradients represent low to high lens-enrichment (in Fold-change) in iSyTE microarray data. Light-red and red gradients represent low to high lens-expression (in fluorescent intensity units) in iSyTE microarray data. FIGURE 5 | "Promising" mRNA targets of Tdrd7-downstream miRNAs in the lens. Heat-maps representing expression of "Promising" mRNA targets of miRNAs identified in (A) "Elevated miRNA group" and (B) "Reduced miRNA group" in Tdrd7-/-lens. The mRNAs identified as "Promising" targets in the "Elevated miRNA" group are identified as reduced in Tdrd7-/-P4 RNA-seq, and are lens-enriched at P4 and/or P30. In the "Reduced miRNA" group the "Promising" targets are identified as elevated in Tdrd7-/-P4 RNA-seq and are lens-enriched at P4 and/or P30. Heatmap keys indicate the following: Yellow and purple color gradient represent low to high differential expression in Tdrd7-/-lens compared to control lens (in Fold-change) in RNA-seq data. Green and red color gradients represent low to high lens-enrichment (in Fold-change) in iSyTE microarray data at P4 and P30. Light-red and red gradients represent low to high lens-expression (in fluorescent intensity units) in iSyTE microarray data at P4 and P30.
FIGURE 6 | Functional gene ontology (GO) assessment of "Top-priority," "High-priority," and "Promising" mRNA targets of differentially expressed miRNAs in Tdrd7-/-lens. Functional GO analysis of (A) 154 miRNA-mRNA pairs representing 98 unique mRNA targets for elevated miRNAs and (B) 123 miRNA-mRNA pairs representing 89 unique mRNA targets for reduced miRNA using DAVID bioinformatics tool. Candidate mRNAs were grouped in GO clusters, which are presented from high to low cluster enrichment score (left panel). One or more relevant GO category in a cluster is presented as bar graph (middle panel). Number of candidate mRNAs in each cluster are grouped as top, high and promising candidates (right panel).
Ppm1e, Loxl1, high-priority candidates Ypel1, Adarb1, Egln2, Hdac4), GO:0005811 "lipid particle" (with top-priority candidate Dgat2, and promising candidates Cyb5r3 and Pnpla2). On the other hand, the top clusters (enrichment score > 1.2, p < 0.05) in reduced miRNA group (i.e., elevated mRNAs in Tdrd7-/-lens) were enriched in GO categories such as GO:0044822 "poly(A) RNA binding", GO:0006915 "apoptotic process", GO:0031625 "ubiquitin protein ligase binding", GO:0006810 "transport" and GO:0006417 "regulation of translation" (Figure 6B). The upregulated genes (i.e., reduced miRNA group) such as Pum2, Ddx3x are commonly found in the GO categories of poly(A) RNA binding and regulation of translation and Ddx3x is also present in the GO for apoptotic process. These GO categories defining the differentially expressed target mRNAs could help explain how misregulation of genes result in cataract. For example, these elevated targets could represent an attempt by lens cells to correct for the defects in post-transcriptional control resulting from Tdrd7-deficiency. Furthermore, elevated transcripts for genes associated with apoptotic process could help explain the lens defects such as large "gaps" (sometime referred to as "vacuoles") in the fiber cell compartment of Tdrd7-/-mice.

Tdrd7-Downstream miRNA-Based Coordinated Control in the Lens
Notably, miRNAs can function individually or in a coordinated manner to mediate regulation of their target mRNAs. A single miRNA can regulate multiple target mRNAs in a specific pathway or multiple miRNAs can converge on a single target mRNA to mediate its control. To gain insights into such coordinated control events potentially mediated by Tdrd7downstream miRNAs in the lens, we examined coregulatory relationships among the "top-priority", "high-priority" and "promising" mRNAs (as defined above). An individual miRNA can have multiple targets that are commonly identified in the same GO category. For example, the upregulated miRNA, miR-138, targets the mRNAs encoding Plekhm1 and Rassf5 that are commonly found in GO category "signal transduction", while miR-345 targets Col4a3 and Timp3 mRNAs that are commonly found in the GO category "basal membrane". On the other hand, we found single miRNA targeting mRNAs involved in varying cellular function. For example, miR-1224 targets eight mRNAs with diverse functions, namely, the enzymes G6pdx (dehydrogenase), Lpin1 (phosphohydrolase) and Wnk2 (kinase), an ER-localized protein with unknown function (Olfm1), a member of Shisha family (Shisa4), a lysosome morphology regulator (Rilpl2), a transmembrane protein (Lemd2) and a solute carrier family protein (Slc35e3). Similarly, miR-15a alone targets eight mRNAs such as Aatk (apoptosis associated tyrosine kinase), Ankrd46 and Anks1 (both associated with cytoskeletal regulations), Egln2 (oxygen homeostasis), Gna12 (signaling), Pacsin3 (protein kinase C involved in linking actin cytoskeleton with vesicle formation), Stub1 (ubiquitin ligase/co-chaperone) and Usp42 (de-ubiquitination). Other such miRNAs that alone target multiple mRNAs are miR-19a, miR-138, miR-328, miR-384, miR-3102 and miR-467a. Conversely, in the elevated miRNA group (i.e., with reduced target mRNAs in Tdrd7-/lens), 39 mRNAs were a common target for multiple miRNAs (Figures 4, 5). For example, both miR-384 and miR-15a have two targets in common, namely those encoding a protein phosphatase (Ppm1e) and a G-protein regulator (Rgs8) in the GO categories "cell projection" and "metal binding", respectively. Further, six mRNAs with varying cellular functions such as Clic5 (involved in actin-based cytoskeletal structures), Ednrb (a receptor molecule), Lin7a (involved in maintaining cell membrane receptors and channels), Ppm1e (serine/threonineprotein phosphatases), Snx11 (a member of the sorting nexin family), 4931406P16Rik (unknown function) are common targets of different combinations of four miRNAs that are elevated in the Tdrd7-/-lens. Interestingly, for upregulated miRNAs, cohorts of mRNA identified in key GO categories were found to be targets of multiple miRNAs. For example, mRNAs involved in Ras signaling (Rasgrp1, Rasl11b, Rassf5, Rhob, Rras) were targets of multiple combinations of miRNAs (e.g., miR-138, miR-15a, miR-19a, miR-3102 and miR-467a). Similarly, mRNAs encoding proteins in solute carrier family (Slc24a3, Slc2a1, Slc35a2, Slc35e3), metallopeptidases (Timp2, Timp3), proteins related to ubiquitin (Ube2v1, Usp2, Usp45) and zinc-finger proteins (Zfp385a, Zfyve21, Zxdc) were targets of multiple combinations of miRNAs.

Derivation of Common Regulatory Networks Between Tdrd7-Downstream miRNA-mRNAs and Other Key Lens Regulators
We next sought to examine which downstream mRNAs are common between Tdrd7-regulated miRNA targets and those of other key regulators implicated in lens development and cataract. We focused on the top-priority Tdrd7-downstream mRNA targets for both elevated and reduced miRNA groups (i.e., mRNA targets of miRNAs elevated or reduced in Tdrd7-/-lens; see above) and examined which of these were also misexpressed in different gene-specific perturbation (either gene-specific knockout, dominant negative or overexpression, representing The edges color/shape represent positive regulation (dark green, arrow) and negative regulation (red, inhibitory sign) between miRNAs or transcription/signaling factors and mRNA targets. The transcription/signaling factor perturbation data in various mouse models is obtained from independent studies using gain-of-function (pink) or loss-of-function (pale green) approaches (details in the Methods section). Node shapes represents miRNA (elevated/reduced), mRNA and transcription/signaling factors.
Frontiers in Cell and Developmental Biology | www.frontiersin.org 13 February 2021 | Volume 9 | Article 615761 FIGURE 8 | Continued FIGURE 8 | Regulatory networks for Tdrd7-downstream negatively regulated miRNAs and their mRNA targets and the relationship of these mRNA target with key lens regulators. (A) Network depicting Tdrd7 downstream negatively regulated miRNAs and the relationship with their predicted target mRNAs based on interpretation of Tdrd7-/-lens transcriptomics data. Deletion of Tdrd7 (pale green node) results in elevation of downstream miRNAs and therefore the relationship between Tdrd7 and these miRNA in normal lens is shown by red color inhibitory edge. In the Tdrd7-/-lens, this elevation of miRNAs results in repression of their target mRNAs, which are found to be reduced. Therefore, the relationship between the miRNAs and these target mRNAs in normal lens is indicated by red color inhibitory edge. Similarly derived networks based on lens microarray data on loss-of-function conditions of (B) Hsf4, (C) Notch, (D) Klf4, (E) E2f1/E2f2/E2f3, (F) Mafg/Mafk, (G) Brg1, and gain-of-function (in fiber cells) of (H) Foxe3 shows the regulatory relationship between these lens regulators and the Tdrd7-downstream negatively regulated miRNAs in the lens. The key provides information on the Edges, Nodes and gene expression conditions and directionality.
loss-of-function or gain-of-function conditions) mouse models of lens defects/cataract that were subjected to meta-analysis in iSyTE . For key lens regulators, we focused on various transcription factors (e.g., Brg1, E2f1/2/3, Foxe3, Hsf4, Klf4, Mafg, Mafk) and signaling pathways (e.g., Notch) (Landgren et al., 2008;He et al., 2010;Saravanamuthu et al., 2012;Gupta et al., 2013;Agrawal et al., 2015;Anand et al., 2015). Based on this analysis, we derived a regulatory network which demonstrates how signaling, transcription and post-transcriptional regulatory pathways mediate combinatorial control over expression of key genes in the lens (Figure 7). To gain insight into the connectivity between Tdrd7-downstream miRNA-mRNA pairs and the above key lens defects/cataract-linked genes, we derived individual regulatory modules that effectively discern the nodeedge relationship for individual gene perturbation conditions (Figures 8, 9). The regulatory network module for Tdrd7 shows its relationship with 18 downstream reduced mRNAs (Alad, FGF2-induced rat lens explants (Wolf et al., 2013). To identify miRNA common to these pathways, we next compared Tdrd7downstream miRNAs (present study) with the FGF2-induced lens explant miRNAs dataset (Wolf et al., 2013). In both studies, miRNAs miR-138, miR-328, miR-345 were commonly found to be significantly elevated, indicating that Fgf and Tdrd7 regulatory pathways potentially converge to mediate control over select downstream miRNAs in lens fiber cells.

Microarray Profiling Identifies Highly Expressed miRNAs in Early Postnatal Lens
Finally, the miRNA profiling by microarray allows an opportunity to assemble a global catalog of the different miRNAs that are robustly expressed in the early postnatal mouse lens. Therefore, we next examined the highly expressed miRNAs (n = 31) defined as those having expression intensity ≥ 500 (p ≤ 0.05) in the control P4 lens (Table 3). This analysis identified 26 new highly expressed miRNAs in the lens while also validating the high expression of previously identified miRNAs (Supplementary Table 4). For example, the miRNAs miR-5105, miR-5109, miR-1298 and miR-378 were newly identified to be highly expressed in the lens. Further, this work offered independent support for the high expression of numerous miRNAs that were previously described in the lens. For example, the present study validates miRNA expression in the mouse lens in agreement with previous studies using microarrays, RNA-sequencing, and in situ hybridization (Supplementary Table 4) (Conte et al., 2010;Karali et al., 2010;Khan et al., 2016). Interestingly this analyses also identified several miRNAs (miR-184, miR-26a, miR-204, let-7b and let-7c) ( Table 3) that were previously found to be misexpressed in human cataractous lenses, again offering independent support that these miRNAs are of significance to lens biology and cataract (Wu et al., 2012(Wu et al., , 2017. Further, many FGF2-regulated miRNAs that were described in a previous study (Wolf et al., 2013) were also found to be highly expressed in the lens in the present study, offering independent support that miRNA function is important in lens development. Thus, this study identifies many miRNAs with high expression in the lens, which can be candidates for future investigations in lens development and cataract pathology.

CONCLUSION
These findings suggest that Tdrd7-downstream miRNAs function to maintain optimal levels and specificity of the mRNA FIGURE 9 | Continued FIGURE 9 | Regulatory network for Tdrd7-downstream positively regulated miRNAs and their mRNA targets and the relationship of these mRNA target with key lens regulators. (A) Network depicting Tdrd7 downstream positively regulated miRNAs and the relationship with their predicted target mRNAs based on interpretation of Tdrd7-/-lens transcriptomics data. Deletion of Tdrd7 (pale green node) results in reduction of downstream miRNAs, thus representing positive control by Tdrd7 in normal lens (therefore indicated by green arrow edge). This results in release of repression of the predicted target mRNAs of these miRNAs in the Tdrd7-/-lens (where these mRNA are found to be elevated) and thus the relationship between miRNA and these mRNAs in normal lens is indicated by red color inhibitory edge. Similarly derived networks based on lens microarray data on loss-of-function conditions of (B) Hsf4, (C) Notch, (D) Klf4, (E) E2f1/E2f2/E2f3, (F) Mafg/Mafk, (G) Brg1, and gain-of-function (in fiber cells) of (H) Foxe3 shows the regulatory relationship between these lens regulators and the Tdrd7-downstream positively regulated miRNAs in the lens. The key provides information on the Edges, Nodes and gene expression conditions and directionality.
transcriptome in the lens, the misexpression of which may contribute to cataract pathology. This is in line with the notion that many small changes orchestrated by several miRNAs may contribute toward fine-tuning gene expression in a cell/tissue (Lim et al., 2005). Regulatory connections between many genes relevant to lens biology and pathology were identified in this study. For example, candidates associated with cell projection, such as Mylk, Olfm1, Plekhg5 and Rgs8, among others, were in the reduced mRNAs (elevated miRNAs) category in the Tdrd7-/-lens. Further, in "intracellular signal transduction" category, Rasgrp1 was identified among other candidates. This is interesting because in an independent study we recently found Rasgrp1 to be involved in the rescue of fiber cell defects in Fgfr2:Pten compound null lenses (Padula et al., 2019). In the GO category "basement membrane", Col4a3, Loxl1, Timp2 and Timp3 were identified as reduced mRNA candidates, which is relevant to the lens capsular defects observed in Tdrd7-defcient lenses, especially because COL4A3 (Collagen Type IV Alpha-3) and LOXL1 (lysyl oxidase-like 1) mutations are associated with cataract and glaucoma in humans (Thorleifsson et al., 2007;Uzak et al., 2013). Further, this analysis also identified -among the promising candidates -the transcription factor CASZ1, which was previously predicted by iSyTE as potentially important in lens  and recently found in a genomewide association study (GWAS) to be associated with cataract in humans (Choquet et al., 2020). Finally, several genes in the GO category "apoptotic process" were found to be elevated (reduced miRNAs) in the Tdrd7-/-lens, offering an explanation for the large vacuole-like gap defects observed in the fiber cells of Tdrd7defecient cataractous lenses. To test individual or combinatorial contributions in the lens of the cohort of miRNAs described here, simultaneous gain/loss of function experiments involving these miRNAs will have be to performed in the future. Also, the present study does not assess the impact of these miRNAs on the level of proteins in the lens. This can be addressed in future studies determining the proteome profile of Tdrd7-/-lens at a stage prior to the onset of lens defects. It will also be interesting to examine the nature of control (whether direct or indirect) that is mediated by Tdrd7 over miRNAs in the lens. In sum, the data presented here indicate that Tdrd7 coordinates distinct downstream regulatory events-either through miRNA-mRNA interactions or through protein-mRNA interactions-to mediate post-transcriptional gene expression control in lens development, misregulation of which causes lens defects and congenital cataract.

DATA AVAILABILITY STATEMENT
The miRNA microarray datasets generated for this study can be found in the Gene Expression Omnibus database under accession number GSE157061.

ETHICS STATEMENT
The animal study was reviewed and approved by University of Delaware Institutional Animal Care and Use Committee (IACUC).

AUTHOR CONTRIBUTIONS
DA, SA, CB, SS, SC, and SL contributed to the generation of the data. DA, SA, and SL analyzed the data. DA and SL wrote the manuscript. All the authors contributed to the revision of the manuscript.