Differential Gene Expression Analysis Reveals Global LMTK2 Regulatory Network and Its Role in TGF-β1 Signaling

Lemur tyrosine kinase 2 (LMTK2) is a transmembrane Ser/Thr kinase whose role has been increasingly recognized; however, when compared to other kinases, understanding of the LMTK2 networks and biological functions is still limited. Recent data have shown that transforming growth factor (TGF)-β1 plays a role in modulating LMTK2 function by controlling its endocytic trafficking in human bronchial epithelial cells. Here, we aimed to unveil the LMTK2 regulatory network and elucidate how it affects cellular functions and disease pathways in either TGF-β1 dependent or independent manner. To understand how the LMTK2 and TGF-β1 pathways interconnect, we knocked down (KD) LMTK2 using small(si)RNA-mediated silencing in human bronchial epithelial CFBE41o- cells, treated cells with TGF-β1 or vehicle control, and performed differential gene expression analysis by RNA sequencing (RNAseq). In vehicle-treated cells, LMTK2 KD affected expression of 2,506 genes while it affected 4,162 genes after TGF-β1 stimulation. Bioinformatics analysis shows that LMTK2 is involved in diverse cellular functions and disease pathways, such as cell death and survival, cellular development, and cancer susceptibility. In summary, our study increases current knowledge about the LMTK2 network and its intersection with the TGF-β1 signaling pathway. These findings will serve as basis for future exploration of the predicted LMTK2 interactions and signaling pathways.

LMTK2 is ubiquitously expressed in human tissues (11). Specifically, LMTK2 mRNA was detected in skeletal muscle, brain, and pancreas (5), while LMTK2 protein was detected in prostate and primary differentiated human bronchial epithelial (HBE) cells (6,12). In polarized human bronchial epithelial cells, LMTK2 was detected at the apical and basolateral domains of the plasma membrane. At the apical membrane, LMTK2 coimmunoprecipitated with CFTR, phosphorylated its residue CFTR-Ser 737 , and facilitated CFTR endocytosis (6). The role of LMTK2 in endocytic trafficking was also supported by its interactions with myosin VI found to be necessary for transport of the transferrin receptor from early endosomes to the recycling vesicles (7,8).
The crosstalk between the LMTK2 and transforming growth factor (TGF)-b1 signaling pathways was reported by two studies. First, LMTK2 was shown to phosphorylate kinesin-1 light chain-2 necessary for Smad2 signaling in HeLa cells (4). Smad2 is a transcription factor that mediates the canonical TGF-b1 signaling that regulates variety of cellular processes (13). Second, TGF-b1 recruited LMTK2 to the plasma membrane by a Rab11-dependent mechanism in human bronchial epithelial CFBE41o-cells (14), facilitating the LMTK2-mediated CFTR phosphorylation. These studies are clinically relevant because elevated TGF-b1 levels, observed in many CF patients, compromise the functional rescue of F508del-CFTR by current generation CFTR modulators (15).
The aim of this study was to elucidate a global LMTK2 regulatory network and its connection with the TGF-b1 signaling pathway in human bronchial epithelial cells. We identified candidate genes regulated by LMTK2 in both TGF-b1-dependent or independent manner. Based on the genes most affected by the LMTK2 knockdown (KD), we predict a complex LMTK2 regulatory network affecting a variety of cellular functions and disease pathways in TGF-b1-dependent and independent manner.

Cell Lines, Cell Culture and Reagents
Parental human bronchial epithelial CFBE41o-cells were cultured in Minimum Essential Medium (MEM; Thermo Fisher Scientific, Walthman, MA, USA), enriched with L-glutamine (Thermo Fisher Scientific), fetal bovine serum (FBS; Gemini Bio-Products, West Sacramento, CA, USA) and penicillin-streptomycin (Thermo Fisher Scientific) (6,16,17). CFBE41o-cells were seeded on collagen-coated 6-well plates (Corning, Corning, NY, USA), at a density of~0.3x10 6 cells. FBS was removed from the media 24h before experiments to avoid any residual TGF-b1, to increase cell polarization, and to promote cell cycle synchronization. Human TGF-b1 (Sigma-Aldrich) was resuspended in the vehicle containing 4 mM HCl and 1 mg/mL bovine serum albumin (BSA, Sigma-Aldrich) and used at a concentration 15ng/ml.

RNA-Mediated LMTK2 KD
Transfection of CFBE41o-cells with siRNA targeting human LMTK2 gene (siLMTK2, siGENOME Human LMTK2 siRNA; Dharmacon, Cambridge, UK) or non-targeting siRNA (siCTRL, siGENOME NonTargeting Pool #2; Dharmacon) was performed using Lipofectamine RNAiMAX Transfection Reagent (Thermo Fisher Scientific), according to the manufacturer's instructions. CFBE41o-cells were plated on collagen-coated cell culture plates and incubated with the optimized transfection mixture containing 50nM of siRNA, at 37°C for 24 h. Next day, culture medium was changed to remove FBS and transfection mixture. Experiments were conducted 2 days after siRNA transfection.

Total RNA Sequencing
For gene expression analysis, we performed three biological replicates for each of the four conditions: (1) CFBE41o-cells transfected with siCTRL and treated with vehicle control; (2) CFBE41o-cells transfected with siCTRL and treated with TGF-b1; (3) CFBE41o-cells transfected with siLMTK2 and treated with vehicle control; and (4) CFBE41o-cells transfected with siLMTK2 and treated with TGF-b1. Total RNA was isolated using the Quick-RNA ™ MiniPrep kit (Zymo Research, Irvine, CA, USA), according to the manufacturer's instructions. The starting material, 1 mg of total RNA was used with the Illumina TruSeq Small RNA sample preparation kit (Illumina, San Diego, CA, USA). The RNA quantity and quality were assessed by Qubit 2.0 Fluorometer (Thermo Fisher Scientific) and Agilent Bioanalyzer Tapestation 2200 (Agilent Technologies, Santa Clara, CA, USA). In short, 3′ and 5' adapters were ligated to total input RNA. Reverse transcription followed by polymerase chain reaction (PCR) was used to create cDNA constructs. This process selectively enriches those fragments that have adapter molecules on both ends. PCR was performed with primers that anneal to the ends of the adapters. cDNA libraries were then purified using QIAquick PCR purification kit (Qiagen). The cDNA libraries were validated using KAPA Biosystems primer premix kit with Illumina-compatible DNA primers and Qubit 2.0 fluorometer. Quality was examined using Agilent Bioanalyzer Tapestation 2200. The cDNA libraries were pooled at a final concentration of 1.8pM. Cluster generation and 75 bp pairedend read sequencing was performed on Illumina NextSeq 500's (Illumina).

RNA Sequencing Analysis
RNA sequencing reads were imported into CLC Genomics Workbench 11 (Qiagen, Valencia, CA, USA) and an initial quality control report of each read was performed. Adapter trimming from the sequence reads was performed using the following adapter sequences: Trimmed sequence reads mapping was performed using human GRCh38/hg38 as a reference genome and using the sequence and annotation from ENSEMBL version 86. Expression values were scaled to TPM and only mapped reads were considered. PCA plot was obtained and differential expression was performed using four different groups: (1) siCTRL-transfected cells with vehicle treatment, (2) siCTRL-transfected cells with TGF-b1 treatment, (3) siLMTK2-transfected cells with vehicle treatment, and (4) siLMTK2-transfected cells with TGF-b1 treatment. Pathway analysis was conducted using Ingenuity Pathway Analysis software (QIAGEN Bioinformatics), using two excluding criteria: max group mean (MGM) >1 and false discovery rate (FDR) p value <0.05. Two different bioinformatics analysis were performed: first, we compared the differential gene expression after LMTK2 knockdown in vehicle-treated cells (siCTRL-transfected cells with vehicle treatment versus siLMTK2-transfected cells with vehicle) treatment and then we evaluated gene expression after LMTK2 knockdown in TGF-b1-treated cells (siCTRL-transfected cells with TGF-b1 treatment versus siLMTK2-transfected cells with TGF-b1 treatment).

Statistical Analysis
Within each condition, the gene lists were obtained using an independent Student's t-test and fold-change analysis to identify upregulated and downregulated genes in siLMTK2-transfected CFBE41o-cells, relative to siCTRL. A significance threshold of P-value<0.05 was used to define modulated genes and signaling pathways of interest. A Venn diagram, obtained in http:// bioinformatics.psb.ugent.be/webtools/Venn/, was used to illustrate genes being up-or downregulated simultaneously in both vehicle-and TGF-b1-treated CFBE41o-cells.

LMTK2 Differentially Regulates Gene Expression After Stimulation With TGF-b1
Our initial aim was to examine the LMTK2 KD using anti-LMTK2 siRNA. CFBE41o-cells were transfected with siLMTK2 or siCTRL and then LMTK2 protein abundance was assessed by WB. LMTK2 protein abundance was decreased by 85% in cells transfected with siLMTK2, compared to siCTRL ( Figure 1). Previous studies have demonstrated that TGF-b1 stimulation does not affect total protein abundance of LMTK2 (14).
Next, total RNA was isolated, from CFBE41o-cells transfected with siLMTK2 or siCTRL and treated with TGF-b1 or vehicle control for 24h, using the Quick-RNA ™ MiniPrep kit (Zymo Research), and analyzed on the Illumina platform. The  Figure 2A). In TGF-b1-treated cells, LMTK2 KD affected the expression of 4,162 genes, of which 2,157 were upregulated and 2,005 were downregulated ( Figure 2B). The above data demonstrate that under the vehicle control condition, LMTK2 affected the expression of 9.7% of the genes for which transcripts were identified (N = 2,506 out of 25,685) and after TGF-b1 stimulation, it affected 16.1% of the identified genes (N = 4,162 out of 25,832). Compared to vehicle control, TGF-b1 stimulation did not significantly change the proportion of genes that were up-or down-regulated by LMTK2 KD. Next, we focused on the genes with the highest fold change after LMTK2 KD in CBFE41o-cells treated with vehicle or TGF-b1 ( Table 1). Our RNAseq data confirm that LMTK2 expression was significantly decreased in cells transfected with siLMTK2, compared to siCTRL. Cyclin-dependent kinase 6 (CDK6) and oxidized low-density lipoprotein receptor 1 (OLR1) were significantly upregulated (> 6-fold change), while torsin family 1 member B (TOR1B), coiled-coil-helix-coiled-coil-helix domain containing 1 (CHCHD1), and cellular repressor of E1A stimulated genes 1 (CREG1) expression were significantly downregulated under both treatment conditions (> -5-fold change). Fibroblast Growth Factor 7 Pseudogene 5 (FGF7P5) was the most upregulated gene (> 160-fold change) in the TGF-b1stimulated cells. In turn, AL358113.1, which encodes an uncharacterized protein, was the most downregulated gene (>-1,000-fold change) in vehicle-treated CFBE41o-cells. Many of the LMTK2 targets are involved in the regulation of important cellular mechanisms and distinct disease pathways, suggesting an important regulatory role of LMTK2.

LMTK2 Plays a Role in the Canonical TGF-b1 Signaling
Our next aim was to elucidate which signaling pathways were affected by LMTK2 KD and how TGF-b1 modulated the effects. The 'molecular mechanisms of cancer' pathway was highly affected by LMTK2 KD in both treatment groups ( Table 2). Indeed, this observation corroborates several independent studies showing that LMTK2 is associated with the susceptibility and severity of colon, gastric, prostate and lung cancer (18)(19)(20)(21). TGF-b1 affected higher number of genes in the 'molecular mechanisms of cancer' pathway, compared to vehicle control (142 and 91 genes for the TGF-b1 and vehicle treatment, respectively; Table 2, Figure 3). We observed an overlap of 88 genes affected by LMTK2 KD in cells treated with TGF-b1 or vehicle control (Figure 4). Of the 54 genes distinctly affected by TGF-b1 treatment after LMTK2 KD, the most downregulated were Phosphoinositide-3-Kinase Regulatory Subunit 2 (PIK3R2), a lipid kinase that phosphorylates phosphatidylinositol, creating second messengers important in growth signaling, Bone Morphogenetic Protein 8B (BMP8B), a secreted ligand of the TGF-b1 superfamily, and Homeodomain Interacting Protein Kinase 2 (HIPK2), a Ser/Thr kinase involved in transcriptional regulation and apoptosis. These three genes are involved in the 'molecular mechanisms of cancer' pathway and are direct mediators of TGF-b1 signaling. These data demonstrate a previously unknown interplay between LMTK2 and TGF-b1 in the molecular mechanisms of cancer in airway epithelial cells. Other pathways affected by LMTK2 KD co-distributed with the treatment. In vehicle-treated CFBE41o-cells, LMTK2 KD affected the signaling pathways associated with B-cell receptor, hepatocyte growth factor (HGF), phosphoinositide 3-kinase (PI3K), and nerve growth factor (NGF). In TGF-b1-treated cells, LMTK2 KD affected the signaling pathways associated with germ cell-Sertoli cell junction, senescence, eukaryotic initiation factor 2 (EIF2), and insulin-like growth factor 1 (IGF-1). Taken together the above data indicate a specific role of LMTK2 in the global regulatory network and the modulatory effect of TGF-b1 on the LMTK2 function.

LMTK2 and TGF-b1 Co-mediate Molecular Mechanisms and Cellular Functions
Next, we examined the molecular mechanisms and cellular functions associated with the differential gene expression after LMTK2 KD in CFBE41o-cells treated with TGF-b1 or vehicle control. We observed that LMTK2 KD affected the cell cycle, cell death and survival, cellular development, and cellular growth and proliferation (Table 3, Figure 5). TGF-b1 modulated regulation of these cellular functions by LMTK2. Indeed, the number of genes affected after LMTK2 KD increased in all mechanisms after TGF-b1 treatment. Furthermore, we observed that the cellular movement mechanisms were dependent on LMTK2 specifically in vehicle-treated cells, while the protein synthesis mechanisms were dependent on LMTK2 in TGF-b1-treated CFBE41o-cells.

LMTK2 and TGF-b1 Co-mediate Diseases Pathways
Lastly, we aimed to examine which pathological processes and organ systems are more susceptible to LMTK2 KD and how it is modulated by TGF-b1 treatment. LMTK2 KD was associated with susceptibility to cancer, organismal injury, and endocrine, gastrointestinal and reproductive system diseases in the TGF-b1 and vehicle control group (Table 4, Figure 5). Although LMTK2 KD was associated with similar disorders in both treatment groups, TGF-b1 affected more genes associated with each disorder, suggesting that it changes the susceptibility to the disease pathways related to LMTK2.

DISCUSSION
The RNAseq analysis in human bronchial epithelial cells elucidated effects of LMTK2 on gene expression and signaling pathways and demonstrated a complex interplay between the LMTK2 and TGF-b1 signaling pathways.
LMTK2 is an integral membrane protein associated with endosomal membranes (7,14). Endocytic trafficking of membrane proteins can affect the intracellular signal transduction, including effects on the activity of various nuclear transcription factors. Indeed, LMTK2 has been previously shown to interact with different protein partners   involved in several cell signaling events and protein transcription (3)(4)(5). LMTK2 activity was shown to promote the kinesin-1mediated intracellular transport of Smad2, a transcription factor in the canonical TGF-b signaling pathway. LMTK2 interacts with Protein Phosphatase 1 (PP1c), leading to its phosphorylation and subsequent inactivation (3)(4)(5). In turn, PP1c inactivates the TGF-b1 signaling pathway, through inhibitory phosphorylation of type I TGF-b1 receptor. Our data provide novel evidence for the role of intracellular trafficking events regulated by LMTK2 in modulating the TGF-b1 signaling during cell death and survival, cellular development, and cancer susceptibility. LMTK2 has been associated with several disorders, including CF (6), neurodegeneration (22)(23)(24), male infertility (25,26) and different types of cancer (10,18,20,21,27). The role of LMTK2 in prostate cancer has been well established. LMTK2 was proposed as a potential biomarker and therapeutic target for prostate cancer (10). Indeed, LMTK2 gene polymorphisms associated with decreased LMTK2 mRNA and protein levels were found in prostate cancer tissues (10,19,28). Furthermore, studies demonstrated that LMTK2 interacts with the androgen receptor (AR) and inhibits its activity (10,12). The role of LMTK2 in other pathological processes, such as the lung, gastric, and colorectal cancer as well as infertility has not been well established. In this study, we demonstrated the role of LMTK2 in several signaling networks involved in the pathophysiology of these diseases.
Our results confirm the role of LMTK2 in cancer d e v e l o p m e n t a n d s u s c e p t i b i l i t y . I n f a c t , L M T K 2 polymorphisms that increase LMTK2 protein levels have been previously reported in cancer in different organs, including the lung (21,29), gastrointestinal tract (18,20), and prostate (19,28,30,31). We observed that LMTK2 influenced the 'molecular mechanism of cancer' pathways and that treatment with TGF-b1 enhanced these effects. Indeed, the role of TGF-b1, a known mediator of cell growth and death, has also been investigated in many different cancer types. TGF-b1 acts primarily as a tumor suppressor, inhibiting cell proliferation and inducing apoptosis of premalignant epithelial cells (32,33). However, at certain stages of carcinogenesis, TGF-b1 may function as a metastasis promoter, inducing epithelial-mesenchymal transition (EMT), cell invasion, and expression of genes that facilitate metastatic colonization of secondary sites (32,33). Currently, it is unknown how and when TGF-b1 changes its function, from a suppressor to a promoter of cancer (32). For this reason, TGF-b1-directed therapies have several limitations (33). Besides 'molecular mechanisms of cancer', LMTK2 was also shown to regulate B cell receptor signaling, HGF signaling, PI3K signaling in B lymphocytes, IGF-1 signaling and senescence pathways, all of which have been previously described to play a role in cancer development (34)(35)(36)(37)(38). Thus, future studies may identify potential checkpoints of the LMTK2 pathway amenable for cancer therapeutic targeting.
Furthermore, we identified multiple molecular signaling pathways and cellular mechanisms that have been associated with distinct pathological disorders after LMTK2 KD. Signaling networks, such as NGF and eIF2, play a role in cellular development, growth, proliferation, and movement. Indeed, NGF signaling was previously shown to downregulate LMTK2 activity, leading to neurite outgrowth, which indicates that LMTK2 inhibits neuronal differentiation (39). Endogenous LMTK2 was found to be phosphorylated upon stimulation with NGF in PC12 cells through a protein kinase C (PKC)dependent mechanism; however, a direct interaction between LMTK2 and PKC was not demonstrated (39). Our study increases the current knowledge on NGF-LMTK2 interaction,  demonstrating that LMTK2 affects the NGF signaling homeostasis through direct regulation of 38 genes (33.3% of the signaling mediators), in vehicle-treated CFBE31o-cells.
Our data suggest that LMTK2 dysregulation may improve the susceptibility to several other pathologies, including organismal injury and abnormalities in the endocrine, gastrointestinal and reproductive systems. Our data shed light on the mechanisms of disease development. However, because our data were generated in human bronchial epithelial cells, further validation of the role of LMTK2 and TGF-b1 in the suggested disorders outside of human airway would be needed in the respective tissue models.
In summary, our data increase the knowledge on the interplay between LMTK2 and TGF-b1 regulatory networks in certain pathological processes affecting different organ systems, including distinct types of cancer. Thus, the study will help to direct future studies to validate the candidates as novel LMTK2 interactors, understand the role of LMTK2 in several disease processes, and identify novel treatments.

DATA AVAILABILITY STATEMENT
The RNAseq data have been deposited in NCBI's Gene Expression Omnibus (GEO) and are accessible through GEO accession number GSE155987.

AUTHOR CONTRIBUTIONS
CF and AS-U designed the study, analyzed the data, and wrote and reviewed the manuscript. DC performed experiments, analyzed the data, and wrote and reviewed the manuscript. NM analyzed the data, and wrote and reviewed the  manuscript. FM analyzed the data and reviewed the manuscript.