Transcriptomic analysis revealed increased expression of genes involved in keratinization in the tears of COVID-19 patients

Recent studies have focused their attention on conjunctivitis as one of the symptoms of coronavirus disease 2019 (COVID-19). Therefore, tear samples were taken from COVID-19 patients and the presence of SARS-CoV-2 was evidenced using Real Time reverse transcription polymerase chain reaction. The main aim of this study was to analyze mRNA expression in the tears of patients with COVID-19 compared with healthy subjects using Next Generation Sequencing (NGS). The functional evaluation of the transcriptome highlighted 25 genes that differ statistically between healthy individuals and patients affected by COVID-19. In particular, the NGS analysis identified the presence of several genes involved in B cell signaling and keratinization. In particular, the genes involved in B cell signaling were downregulated in the tears of COVID-19 patients, while those involved in keratinization were upregulated. The results indicated that SARS-CoV-2 may induce a process of ocular keratinization and a defective B cell response.


Results
The demographic and clinical characteristics of the patients suffering from COVID-19 enrolled are summarized in Tables 1 and 2. The mean age of control group was 73.5 ± 10.7 years (75% males; 25% females). No statistical difference was found between COVID-19 and healthy groups in terms of age, gender and race (p > 0.05).

RNA-seq analysis between healthy individuals against Covid-19 patients.
The RNA-seq analysis revealed 25 genes that differ statistically between healthy individuals and COVID-19 patients (Table 3). Among them, 13 genes were upregulated while 12 genes were downregulated. Interestingly, 20 genes, the most of them, had more than a twofold change. Anyway, the remaining 5 genes had a onefold change. The changes in the behavior of the differential expressed genes between all the samples was depicted in the heatmap in Fig. 1. As expected, the dendogram obtained by the plot put the closest association between the control individuals themselves and the COVID-19 patients themselves. Furthermore, in order to inspect the role of the up-and downregulated differentially expressed genes, the PANTHER Classification System 16 for the Gene Ontology Biological Processes was used. The enriched Gene Ontologies of the Biological Processes obtained by PANTHER (Table 4) highlighted the "Keratinization" (False Discovery Rate (FDR) = 1.83 × 10 -3 ) for the up-and the "Regulation of B cell activation" (FDR = 1.49 × 10 -3 ), "Negative regulation of immune system process" (FDR = 2.21 × 10 -3 ) and "Regulation of inflammatory response" (FDR = 1.64 × 10 -3 ) for the down-regulated genes.
To each gene in Gene Symbol column was associated the corresponding name in Gene Name column with the bitr function of the cluster Profiler package of Bioconductor. The expression of each gene in healthy individuals and COVID-19 patients was highlighted in the columns Healthy Expression and COVID-19 Expression, respectively. The column Fold Change shows the changes in the expression level between healthy and COVID-19 individuals while the q-Value column shows significance of the difference in the dataset (q < 0.05).
For each of the most specific biological processes found with PANTHER, the differentially expressed genes that are included were highlighted: a pathway is enriched for the upregulated genes and 5 genes are included; 3 pathways are enriched for the downregulated genes and 6 genes are included. The False Discovery Rate (FDR) of each pathway is lower than 0.05.  18 ) the mRNAs expression level of the genes found in our analysis in B and T cells. As depicted in red bars of the pyramid plot in Fig. 2

Discussion
The eye has been considered as a potential site for SARS-COV-2 viral infection and dissemination 10,19 . Moreover, involvement of the eye seems to be related more likely to severe forms of COVID-19 disease and often precedes the systemic symptoms or even is the only sign of the disease 15 . In our sample all diseased subjects had a severe SARS-CoV-2 infection; most of them had pneumonia (78.9%) and 70% of cases died with a mean of 26 ± 12 days of hospitalization.
Our study aimed to provide a panel of gene expression in tears of patients affected by SARS-CoV-2 infection comparing with healthy subjects to understand the gene expression pattern in the eye, which can be a starting point site of infection, as well as a concomitant involved organ of the disease. Given the protective function of tear film preserving the homeostasis and health of the conjunctiva and the avascular cornea, abnormalities in concentrations of proteins and inflammatory mediators in lacrimal secretions have been observed in infections, surgery and trauma. Indeed, antimicrobial factors have been well described in tears including lysozyme, lactoferrin, transferrin, ceruloplasmin, IgA, IgG, IgE, complement, glycoprotein, and anti-proteinase, which are found in the aqueous layer of the tear film. As already reported, immunoglobulins play a key role in defense of bacterial, viral, and parasitic infection. IgA, which is the primary immunoglobulin found in tears produced by conjunctiva and lacrimal gland, usually increases during infectious or inflammatory conditions of conjunctiva 20 . Moreover conjunctivitis has shown a rise in inflammatory mediators (IL-1β, TNF-α, and MMP-9) and the activation of proinflammatory mitogen-activated protein kinase (MAP-K) pathways 21 . In detail, tear film has three major layers, such as the inner mucin layer, the middle aqueous layer and the outer lipid layer 22,23 . The inner layer, formed by mucins secreted mainly by the goblet cells in the conjunctival epithelium, has stabilizing role of aqueous layer. It is composed by immunoglobulins, urea, salts, glucose, and proteins as well. The aqueous layer, essential for maintaining hydration and health of the ocular surface, contains proteins, metabolites, inorganic salts, glucose, oxygen, and electrolytes (magnesium, bicarbonate, calcium, urea). The lipid layer, fundamental for controlling tear evaporation, contains cholesterol, wax esters, fatty acids, and phospholipids 23,24 . Other tear film components include lysozyme with its bacteriolytic role, lactoferrin that is able to sequester iron from the bacteria thus stopping their growth. Mucins and glycoproteins secreted by goblet cells have a known role in ocular www.nature.com/scientificreports/ defense from external environment, preventing attachment of pathogens to ocular surface. In our sample, RNAseq analysis revealed 25 genes differing statistically between healthy individuals and patients with a diagnosis of COVID-19. In detail, in the COVID-19 subjects 13 genes were upregulated, while 12 genes were downregulated. We found that different genes involved in the function of B cells were downregulated, as also confirmed by Panther analysis that evidenced that the biological process "Regulation of B cell activation", but also "Regulation of inflammatory response" and "Negative regulation of immune system process", were enriched for the Figure 1. Heatmap of the differentially expressed genes between all the individuals. The heatmap represent the genes differentially expressed after the RNA-seq analysis was concluded. Moreover, the logarithmic correction made on the counts allows to appreciate the differences between the samples. As expected, the dendogram on the columns shows the closest similarity between the healthy subjects each other (CTR) and between the COVID-19 patients themselves (COVID).  26 . IGHG1 decrease was reported also in Fuchs endothelial corneal dystrophy 27 . Moreover, also the genes FCRL4 and FCRL5, that are receptors for IgA and IgG, respectively, were downregulated and have also a role in viral infections 28,29 . DNASE1L3, that was downregulated in COVID-19 group, encoded for DNase γ, member of DNase I family of endonucleases. It was found to be expressed in germinal center B cells and stimulated B cells. It is involved in the somatic hypermutation of immunoglobulin variable region genes that occurs in the germinal center B cells during immune responses, and then contributed to the immunoglobulin V gene diversification 30,31 . DNASE1L3 is also involved in the release of cytokines after activation of the inflammasome 32 .
BANK1 also was found to be downregulated in the tears of COVID-19 patients. It is a positive regulator of B cell signaling through the induction of calcium mobilization after the activation of B-cell antigen receptor 33 . Also CLEC17A, expressed in B cells with function of adhesion to epithelial cells 34 , was downregulated in our analysis.
After sequencing we reported a downregulation of the ATM gene which synthesizes serine/threonine-protein kinase that act in the cell following DNA damage 35 . ATM has an important role in both T and B cells function. In particular, the loss of ATM in B cells caused the reduction in germinal center frequency and size in response to immunization and apoptosis of B cells 36 . Moreover, ATM is involved also in T cell development 37 . Moreover, it was found that RNA viruses can cause DNA damage and genetic instability in host cells modulating components of the DNA damage response, such as ATM 38 . ATM gene mutation were found in ocular adnexal marginal zone B-cell lymphomas and uveal melanoma [39][40][41] . In addition ATM gene inhibition reduces herpes virus corneal infection and particularly epithelial infection and stromal disease.
We found the upregulation of TLE1, MAL and ZBTB16 in the tears of COVID-19 patients. The TLEs family is composed by co-repressors expressed in T cells where they are required for CD8 + T cell lineage choice 42 . MAL is involved in apical transport of proteins in polarized epithelial cells 43 . MAL has been shown to be implicated in lytic plaque formation and viral spread in oligodendrocytes infected with Herpes simplex virus type 1 44 . It is involved also in T cell functions 45,46 . ZBTB16 is a gene also known as PLZF (promyelocyticleukemia zinc finger). It is important in the function and development of immune system and may enhance T cell responses 47 . www.nature.com/scientificreports/ Our seq analysis revealed a significant upregulation of IL19. The upregulation of IL-19 was already found in blood of COVID-19 patients 48 . Li et al. described the gene expression profile from patients with Th cellmediated autoimmune noninfectious uveitis and some cytokines, including IL-19, were found to be highly expressed proving the inflammatory status of the eye 49 . The upregulation of such inflammatory interleukin in patients with the coronaviridae disease shows an inflammation of conjunctiva. Interestingly, IL-19 is produced by macrophages, B-cells but also by keratinocytes, and interestingly, it can also act on keratinocytes, suggesting to be involved also in keratinocyte hyperproliferation 50,51 . This aspect is particular interestingly, considering that we found the upregulation of genes involved in keratinization as confirmed also by the enrichment of the same biological process.

Figure 2. Gene expression in T and B cells retrieved in Human Protein Atlas database. The bars represent the mRNA expression level of each gene in T cells (on the left) or in B cells (on the right). The bar is red in the cells in which the
We also found an upregulation of keratin 17 and keratin 78 expression levels in tears of our enrolled patients. A similar KRT17 upregulation was detected by Kulkarni et al. 52 in limbal epithelial stem cells of diabetic patients using deep sequencing analysis if compared with healthy individuals. The authors speculated that KRT17 dysregulation could be involved in the typical corneal modifications of diabetic subjects. From our analysis KRT17 and KRT78 were upregulated likely due to their potential role in ocular surface modifications during SARS-CoV-2 infection.
Our findings showed a significant increase in FLG levels compared to the healthy subjects. FLG has also an important role in maintaining the integrity of stratum corneum of epidermidis and its loss-of-function mutations have been associated with atopic dermatitis, due to a reduced skin hydration and skin barrier function 52 . A previous study indicated that FLG was not expressed in normal conjunctiva, while it increased in moderate and severe forms of parakeratinization of the conjunctiva 53 .
We found the upregulation of both SPRR1B and SPRR3. They are part of the small proline-rich proteins (SPRRs), which constitute cornified cell envelope precursors 50,54 . SPRR1B is a stress-induced transcript on the ocular surface that was shown to be upregulated in conditions of pathologic keratinizationin, in both evaporative and immune-mediated, aqueous-deficient dry eye disease. SPPR1B can be also considered a biomarker for pathology such as Sjögren's syndrome and for the study of the molecular mechanisms of squamous metaplasia. Squamous metaplasia causes pathologic keratinization of the ocular surface in response to disease processes that are autoimmune mediated or caused by infection. Moreover, pro-inflammatory cytokines induced the expression of SPRR1B 55 .
SPDEF, upregulated in our analysis, is a transcription factor that promotes the differentiation of goblet cells in the conjunctiva epithelium 56 . Goblet cells produce and secrete mucins, needed to lubricate the ocular surface. Goblet cell secretions are essential in order to maintain tear stability and ocular surface homeostasis. Goblet cells in the conjunctiva play an essential immunomodulatory role 57 . Increased expression of SPDEF was associated to goblet cell hyperplasia 58 and upregulation of goblet cell-associated genes 59 . Atopic keratoconjunctivitis was associated with goblet cell hyperplasia 60 .
TFF3 was found to be upregulated in our analysis. It is expressed in various ocular tissues, also in goblet cells and has a role in corneal wound healing 61 . Its increase in tear film was found after inflammatory factors or ocular surface damage, such as those in dry eye disease, in experimental models 62 . It increased also in herpetic keratitis 63 .
LTF, produced in the acinar cells of the lacrimal gland, is normally present in tears of humans and is not dependent on age and sex. It has been proved that LTF is reduced in some ocular and systemic disease such as dry eye-related keratopathies, herpes simplex keratitis, chronic irritative conjunctivitis keratoconjunctivitis sicca 64 .
PLA2G2A, which is downregulated in tears of our diseased sample of patients, has been previously described to be involved in tear film stability and the consequent integrity of tear ocular surface 65 .
ALOX15B is a lipoxygenase enzyme. It was shown that LPS bacteria and the proinflammatory cytokines may induce overexpression of ALOX15B 66 . It represents the predominant 15-LOX protein in human cornea, and its product induced apoptosis in a dose dependent manner 67 .
Type II transmembrane serine proteases (TTSPs) have the ability to cleave surface proteins of viruses, including SARS-CoV-2 and influenza viruses, leading to the viral invasion. In particular, SARS-CoV-2 can bind TMPRSS2 of the host cell. This protease facilitates the viral attachment to the surface of targets cells by cleavage and priming of a spike protein of coronaviruses (S protein), at the S1/S2 and the S2' site. This process seems to be the first essential step for the virus entry in the human body, thus demonstrating a key role of TMPRSS2 activity on SARS-CoV-2 spread 68 . Interestingly, in our RNA-seq analysis we found TMPRSS11E upregulated in tears of COVID-19 patients, leading to the hypothesis of a possible additional role of this protease in the Sars-Cov-2 entry in the host cells. Indeed, it was demonstrated that TMPRSS11E also known as DESC1, activated the S proteins of emerging Coronavirus 69 .
PARP15 is a neurodegeneration mediator and its inhibition seems to improve corneal epithelial innervation and wound healing in diabetic rats 70 . It is known that PARPs may play a role in inflammation and virus infections, indeed coronavirus may have the ability to reverse ADP ribosylation induced by PARP in order to counteract host-virus defense response 71,72 .
We can speculate that the down or up regulation of some of the investigated genes already known as involved in ocular inflammation processes are related to the presence of the virus on the ocular surface. In particular the up regulation of IL19 gene, associated to non infectious uveitis and conjunctivitis could be related to SARS-CoV-2 infection in our patients. TFF3 also found to be upregulated in our sample has already been related to corneal inflammation in dry eye disease herpetic keratitis. The inhibition of ATM gene, downregulated in our study, has been found related to a reduction of herpes virus corneal infection.
The main limitation of our study is the lack of a sub analysis differentiating COVID-19 patients with and without conjunctivitis to correlate ocular clinical conditions and the RNA seq analysis.
In conclusion, our study reported in tears of patients suffering from SARS-CoV-2 infection a downregulation of genes involved in B cell signaling, while the genes involved in keratinization were upregulated (Fig. 3). These results are particularly interesting because to our knowledge this is the first study that suggest a possible www.nature.com/scientificreports/ molecular mechanism responsible of the keratoconjunctivitis reported in patients affected by COVID-19 73,74 . It would be interesting to investigate if IL19 increase may link the inflammation process with the keratinization.

Materials and methods
Patients enrollment. In this observational study, approved by our Institutional Review Board (prot. n. The following anamnestic data about all patients enrolled were collected: • Demographic data (age, sex, race);  • Symptoms (shortness of breath, cough, fever, myalgia, pneumonia); • SaO2; • Evidence of conjunctivitis (presence of bilateral conjunctivitis was defined as red eyes: macroscopic signs of conjunctival congestion and was confirmed by using a portable slit lamp (SL-17; Kowa Ltd., Tokyo, Japan) and if present the onset of the ocular disease was registered: pre-admission, on admission, during hospitalization); • Necessity of non-invasive (facemask oxygen)/invasive mechanical ventilation; • Clinical outcomes (recovery, discharge from hospital, death).
Tear sampling. Tear secretion was measured by using the Schirmer's Test (Schirmer strips; Whatman, Maidstone, UK). Schirmer's papers have been kept in the lower lid margin for 5 min. Tear secretion was measured as the length of the wet strip (in millimeters). Schirmer's tear fluid collection was carried out simultaneously in both eyes. The paper strips from eyes were placed in Eppendorf tubes containing 500 μl of water with RNase inhibitors (Diethyl pyrocarbonate (DEPC)-treated water, Ambion, CA, United States) and strongly vortexed for 5 min and frozen in dry ice. Upon arrival to the laboratory both samples have been stored at − 80 °C. RNA extraction and library preparation. The library was prepared in according with Illumina protocols and how already described 76 . Briefly, the total RNA was extracted with the Maxwell RSC simplyRNA Cells Kit (Promega, Madison, WI, USA) following the manufacturer's instructions. The library preparation was carried out according to the TruSeq RNA Exome protocol (Illumina, San Diego, CA, USA) 77 .
The cDNA was synthesized with the SuperScript II Reverse transcriptase (Invitrogen, Carlsbad, CA, USA) and the 3' end was later adenylated. The library was amplified with PCR and validated. Two steps of hybridization followed the first PCR and later one more PCR was performed. The MiSeq Instrument (Illumina) was used to sequence the libraries.
Bioinformatics analysis. The obtained raw data was checked for quality using the fastQC tool and the needed trimming for adapters and low quality bases scores was performed by Trimmomatic (Usadel Lab, Aachen, Germany) 78 (version 0.38). The reads were than aligned to the version GRCh38 of the human reference genome with the Spliced Transcripts Alignment to a Reference (STAR) RNA-seq aligner 79 . After sorting the reads with STAR, the counting of the reads was made using the htseq-count python package 80 . The differentially expressed genes were finally obtained using the Bioconductor package DESeq2 in R language 81 . The Benjamini-Hochberg post-hoc procedure with threshold 0.05 was used to filter out the false positives. The bitr function of the clusterProfiler package of Bioconductor associated each gene symbol with its corresponding name 82 . Panther was then used to enrich the genes with the biological processes in which they are involved. Finally, the Human Protein Atlas database was used to find the mean level of expression of the mRNAs coding for these genes in B cells and in T cells. After logaritmic conversion, the values were depicted in R using the function ggplot2 83 .

Data availability
The datasets generated during and/or analyzed during the current study will be available from the corresponding author on reasonable request.