Introduction

Tear fluid, produced by the lacrimal glands and spread across the ocular surface, serves functions beyond basic lubrication of the eyes; it also serves as an indicator of the biochemical and physiological condition of the conjunctiva and cornea1,2. Increasing attention has been given to tear fluid as a source of biomarkers relevant not only to ocular pathologies but also to systemic conditions. Traditionally, the discovery of diagnostic biomarkers has concentrated on blood, serum, and plasma, yet these fluids often require complex preprocessing because of blood cells and their intricate proteome3,4. By contrast, tears are relatively easy to obtain through noninvasive methods, such as Schirmer strips or microcapillary tubes, making tear collection more convenient for both clinicians and patients5. While tears contain proteins, lipids, and metabolites, the detection of nucleic acids, especially RNA species, offers a sensitive snapshot of early gene expression changes6. This early window of detection is crucial because shifts in gene expression often precede protein-level alterations, potentially enabling earlier diagnosis and intervention in ocular diseases7,8,9.

Dry eye disease (DED) is a prevalent, multifaceted condition marked by tear film instability and high osmolarity, triggering an inflammatory response on the ocular surface10,11,12. Several studies suggest that inflammatory mediators, including pro-inflammatory cytokines and chemokines, are crucial in causing and perpetuating ocular surface damage in DED13,14,15. Similarly, in glaucoma management, benzalkonium chloride (BAK), a common preservative in topical medications, can exacerbate ocular surface toxicity, highlighting the need for early biomarker detection16. Therefore, a thorough RNA analysis of tear fluid could improve specificity and sensitivity in diagnostics, potentially enhancing current clinical evaluations for DED and BAK-related toxicity.

The advantage of RNA as a biomarker lies in its dynamic reflection of gene transcription within cells, which contributes to the tear film. Both messenger RNA (mRNA) and noncoding RNAs (including microRNAs) have been detected in tears17,18,19,20,21. These molecules often reflect disease severity, making them promising biomarkers for individualized (patient-specific) diagnostics and for monitoring therapeutic outcomes. mRNA in bodily fluids, such as peripheral blood, saliva, and seminal fluid, exhibits tissue-specific expression profiles. This enables precise fluid identification through quantitative RT-PCR, as evidenced in forensic applications leveraging stable mRNA biomarkers22. This molecular specificity highlights the diagnostic utility of tear fluid mRNA in detecting early transcriptional alterations associated with ocular surface pathologies, enabling non-invasive biomarker discovery for conditions like DED and glaucoma. However, to achieve reliable RNA-based testing, consistent sample collection, stabilization, and quantification methods are paramount, an endeavor complicated by the limited volume and heterogeneous composition of tear fluid23.

Although considerable interest exists in tear fluid RNA, analytical variability presents a major challenge. This variability can arise during sample collection (for example, reflex tearing triggered by external stimuli), storage, and RNA extraction methods24,25. Moreover, identifying stable reference genes, commonly referred to as housekeeping genes, for data normalization is essential for the reliability of gene expression analyses26. Traditionally, genes like GAPDH, ACTB, and 18S have been utilized across various tissues and cell lines. However, their expression may not remain consistent in disease-specific contexts or in the unique environment of tear fluid27,28. Even in neuroinflammation models, systematic validations have indicated that so-called “universal” housekeeping genes can exhibit significant variability under various stressors or treatments29,30,31. The same caution should be applied to tear fluid, highlighting the necessity for a thorough evaluation of potential reference genes before drawing any conclusions about changes in target gene expression.

In inflammatory ocular conditions such as DED and glaucoma patients using BAK-preserved topical treatments, Caspase-1 (Casp1) and Apoptosis-associated speck-like protein containing a CARD (ASC) are increasingly recognized as critical mediators of inflammasome activation14,32. These molecules are crucial for the immune response, facilitates the maturation of the pro-inflammatory cytokine interleukin-1β, while ASC acts as a scaffold protein essential for inflammasome assembly32. Although protein-level assays most directly reflect functional activity, changes in Casp1 and ASC mRNA might serve as an early indicator of inflammasome activation, often preceding overt protein-level modifications in conditions like DED or glaucoma14. Consequently, evaluating the RNA levels of these key molecules in tears may enable the detection of early or subclinical stages of ocular inflammation, potentially assisting in the development of more targeted therapeutic strategies.

The objective of this study was to systematically evaluate the expression stability of seven candidate reference genes in tear fluid collected from individuals with DED, glaucoma patients using benzalkonium chloride (BAK)-preserved topical medications, and healthy controls. By employing five established computational algorithms (geNorm, NormFinder, BestKeeper, the comparative ∆CT method, and RefFinder), we aimed to identify reliable internal controls for normalizing RNA expression in tear samples, thereby establishing a robust foundation for accurate transcriptomic analysis in ocular surface diseases.

Results

This study systematically evaluated the stability of seven candidate reference genes (18S, RER1, ACTB, GAPDH, PGK1, UBC, AP3D1) for RT-qPCR analysis of tear fluid RNA from 24 participants (n = 8 per group: DED, glaucoma patients using BAK-preserved topical medications, healthy controls). Tear samples were collected noninvasively using Weck-Cel® spears under standardized conditions, with total RNA extracted and quantified (Supplementary Table 3). Candidate genes were selected based on their consistently high expression in transcriptomic data from conjunctiva, cornea, and eyelid tissues33 and prior use in ocular surface studies34, ensuring relevance to tear fluid RNA, primarily from exfoliated conjunctival and corneal epithelial cells. Gene expression stability was assessed using five established algorithms (geNorm, NormFinder, BestKeeper, comparative ΔCT method, RefFinder) to identify reliable normalization controls for tear-based transcriptomic studies.

Demographic characteristics and clinical parameters

This study evaluated a total of 24 eyes from 24 participants who met strict inclusion criteria: the DED group included individuals with an OSDI score ≥ 13 and/or corneal staining (NEI) 3; the glaucoma group included patients treated with BAK-preserved hypotensive medications; and the control group included asymptomatic individuals without ocular surface disease (Table 1). The mean age was 36.75 ± 11.34 years in the control group (n = 8), 56.63 ± 13.89 years in the DED group (n = 8), and 72.25 ± 7.01 years in the glaucoma group (n = 8). Female participants represented 50% of the control group and 37.5% of both DED and glaucoma groups. Supplementary Tables 1 and 2 provide clinical characteristics of the DED and glaucoma groups.

Table 1 Demographic characteristics and ocular surface parameters in patients with DED, glaucoma, and healthy controls.

Clinically, OSDI scores were significantly higher in both the DED (31.25 ± 15.84, p = 0.0024) and glaucoma (38.89 ± 17.53, p = 0.0002) groups compared to controls (5.35 ± 4.93). NEI corneal staining scores were markedly increased in the glaucoma group (4.67 ± 2.81, p = 0.0003), with more moderate elevations observed in the DED group (2.57 ± 0.98, p = 0.0227) relative to controls (0.25 ± 0.46). Tear film stability, measured by the first non-invasive tear break-up time (NIKBUT), was significantly reduced in the DED group (5.23 ± 2.65 s, p = 0.0118) compared to controls (8.69 ± 1.21 s), whereas values in the glaucoma group (7.11 ± 2.71 s, p = 0.348) did not differ significantly. Ocular redness scores were significantly higher in both DED (1.65 ± 0.46, p = 0.0013) and glaucoma (1.78 ± 0.44, p = 0.0011) groups compared to controls (0.86 ± 0.24). Meibography scores also demonstrated increased gland dropout in both DED (1.57 ± 1.13, p = 0.0279) and glaucoma (1.44 ± 1.12, p = 0.0413) groups compared to controls (0.25 ± 0.53).

Expression level of the candidate reference genes

We performed RT-qPCR to assess the transcriptional expression levels of seven candidate reference genes in all 24 samples from DED (n = 8) and glaucoma (n = 8) patients, as well as controls (n = 8), as outlined in the methods section. Expression levels are represented by the raw quantification cycle (Cq) values. The Cq values for all seven reference genes across the samples varied from 19.01 to 35.19. Notably, 18S showed the lowest Cq value, while AP3D1 had the highest, indicating that 18S is the most prevalent reference gene in the collected tear samples, as shown in Fig. 1.

Fig. 1
figure 1

Quantification Cycle (Cq) values of the reference genes across all samples. Bars represent mean ± SD of the Cq values in control (n = 8), DED (n = 8), and glaucoma (n = 8) tear samples.

The stability of these seven genes was evaluated using five analytical methods: geNorm, , ΔCt method, NormFinder, BestKeeper, and RefFinder. The results are displayed in Table 2, reflecting the stability parameters for each platform. 18S consistently emerges as the top-ranked or second-ranked reference gene across all metrics. RER1 and ACTB are similarly high-ranking, while UBC and AP3D1 exhibit poorer performance. AP3D1 demonstrates a mixed ranking; BestKeeper identifies it as stable, yet other methods place it among the least stable.

Table 2 Stability values and ranking assessed using the ∆CT method (mean SD values), BestKeeper (SD values), NormFinder (stability values), geNorm (M values), and RefFinder (Geomean) for the seven candidate reference genes from the control (n = 8), DED (n = 8), and glaucoma (n = 8) tear samples.

Analysis of gene expression stability

geNorm

The geNorm algorithm evaluates gene stability by calculating the average pairwise expression ratio, known as the M value28. The analysis conducted using qbase+35, showed that 18S ribosomal RNA (18S) exhibited the lowest M value (0.436), identifying it as the most stable reference gene. Conversely, AP3D1 demonstrated the highest M value of 1.036, categorizing it as the least stable for normalization purposes. According to geNorm criteria, M values below 0.5 indicate high stability, values ranging from 0.5 to 1.0 denote moderate stability, and values exceeding 1.0 signify low stability.

In our dataset (Table 2), all evaluated reference genes, except for AP3D1 (M = 1.036), showed M values below 1, indicating generally acceptable stability. Notably, both 18S and RER1 achieved the lowest M value of 0.436, further supporting their suitability as the most stable reference genes according to geNorm analysis (Fig. 2A). geNorm ranked the stability of expression from most to least stable reference genes as follows: 18S < RER1 < ACTB < GAPDH < PGK1 < UBC < AP3D1. Additionally, geNorm assesses pairwise variation (V) to determine the necessity of incorporating additional reference genes for optimal normalization. A V value below 0.15 is considered acceptable for reliable normalization36. In our analysis, the pairwise variation V3/4 was calculated to be 0.135, indicating that the inclusion of the three most stable reference genes (18S, RER1, and ACTB) is required for accurate normalization (Fig. 2B).

Fig. 2
figure 2

Gene expression stability and the optimal number of reference genes needed for normalization using GeNorm. (A) Ranking of reference genes, the least stable gene is identified by the highest M value, with lower M values indicating greater stability. (B) The analysis of pairwise variation (Vn/Vn + 1) identifies the optimal number of genes for effective normalization, with a V value of less than 0.15 considered acceptable. The calculated value of V3/4 is 0.135, indicating that incorporating the three most stable reference genes is crucial for accuracy normalization. Control (n = 8), DED (n = 8), and glaucoma (n = 8) tear samples.

Comparative ∆CT method

The comparative ∆CT method was employed to evaluate the relative expression of gene pairs within each sample, facilitating the identification of suitable reference genes for expression analyses37. This approach involves calculating the ∆CT values, defined as the differences in quantification cycle (Cq) values between pairs of reference genes, alongside their corresponding Standard Deviation (SD) values. Minimal variability in ∆CT indicates stable expression between gene pairs, ensuring reliable normalization across diverse samples and experimental conditions. Conversely, significant fluctuations in ∆CT suggest variability in the expression of at least one gene within the pair, undermining normalization accuracy.

In our analysis, the ∆CT method ranked the stability of reference genes from most to least stable as follows: 18S < RER1 < ACTB < GAPDH < PGK1 < UBC < AP3D1, which is consistent with the rankings obtained via geNorm (Fig. 3). Table 2 shows the average SD values resulting from the ∆CT method for each reference gene. Notably, 18S ribosomal RNA (18S) exhibited the lowest average SD of 0.750, designating it as the most stable gene according to the ∆CT approach. In contrast, AP3D1 demonstrated the highest average SD of 1.645, categorizing it as the least stable reference gene.

Fig. 3
figure 3

Stability ranking of reference genes assessed by the comparative ∆CT method. A high delta Ct value indicates instability in the expression of the gene. The dots illustrate the average Standard Deviation (SD) values for each candidate reference gene. Reference genes are ordered from most stable (lowest SD) to least stable (highest SD). Control (n = 8), DED (n = 8), and glaucoma (n = 8) tear samples.

BestKeeper

The BestKeeper software computes raw Cq-based parameters, including standard deviation (SD), coefficient of variance (CV), and the BestKeeper Index, to assess gene stability. Table 2 presents the BestKeeper-derived SD values for the assessed reference genes27. Typically, genes with an SD greater than 1 are considered unstable and should be avoided in further analyses. In our study, AP3D1 has the lowest SD (0.575), indicating high stability, while UBC has the highest SD (1.060), showing lower stability (Fig. 4). According to the rankings derived from this method, AP3D1 holds the top position, whereas UBC is identified as the least stable, suggesting it should be excluded from subsequent normalization procedures based on BestKeeper standards. Expression stability by BestKeeper was ranked as follows: AP3D1 < 18S < GAPDH < ACTB < RER1 < PGK1 < UBC.

Fig. 4
figure 4

Stability ranking of candidate reference genes based on BestKeeper analysis. Each reference gene’s standard deviation (SD) values are plotted, with lower SD values indicating greater expression stability. Genes are ordered from the most stable on the left to the least stable on the right, highlighting AP3D1 as the most stable and UBC as the least stable reference gene. Control (n = 8), DED (n = 8), and glaucoma (n = 8) tear samples.

NormFinder

NormFinder software evaluates both intra- and inter-group variations to generate a stability value, thereby facilitating the identification of optimal reference genes across various experimental conditions38. In our initial analysis, 18S (S = 0.218) and RER1 (S = 0.395) emerged as the most stable reference genes, demonstrating minimal expression variability. Conversely, UBC and AP3D1 exhibited lower stability, as shown in Fig. 5. The overall ranking of gene stability as determined by NormFinder is as follows: 18S < RER1 < ACTB < GAPDH < PGK1 < UBC < AP3D1. Detailed stability values (S) for each reference gene are presented in Table 2.

Fig. 5
figure 5

Stability ranking of candidate reference genes based on NormFinder analysis. The stability values (S) for each reference gene are plotted, with lower values representing greater stability. Genes are ranked from most stable (left) to least stable (right), with 18S emerging as the most stable and AP3D1 as the least stable reference gene. Control (n = 8), DED (n = 8), and glaucoma (n = 8) tear samples.

RefFinder

RefFinder is a comprehensive online tool that integrates results from the comparative ∆CT method, BestKeeper, NormFinder, and geNorm to generate an overall stability ranking for reference39. It calculates a geometric mean (Geomean) for each candidate reference gene based on the weighted scores obtained from these four methodologies, with smaller Geomean values indicating greater stability. As illustrated in Table 2, RefFinder ranked 18S as the most stable reference gene with the lowest Geomean value of 1.189, followed by RER1 at 2.115. In contrast, UBC (6.236) and PGK1 (5.233) were identified as the least stable reference genes across all analytical tools. Notably, AP3D1, which was ranked highest by BestKeeper alone, received a moderate RefFinder rank of 5 (Geomean = 4.304), underscoring the enhanced insights achieved through combined evaluations of multiple methods (Fig. 6).

Fig. 6
figure 6

Overall stability ranking of candidate reference genes as determined by RefFinder analysis. The geometric mean (Geomean) values for each reference gene are depicted, with lower values indicating higher stability. Genes are ordered from most stable (18S) to least stable (UBC), reflecting the integrated assessment from multiple analytical methods. Control (n = 8), DED (n = 8), and glaucoma (n = 8) tear samples.

To derive an overall ranking that encompasses each analytical approach, the geometric mean of each reference gene’s position across all individual programs was calculated. The comprehensive expression stability determined by RefFinder was: 18S < RER1 < ACTB < GAPDH < AP3D1 < PGK1 < UBC. Our final ranking and RefFinder both confirmed 18S, RER1, and ACTB as the best reference gene combination in the tear samples of controls, DED, and glaucoma groups.

Reference gene validation analysis

To experimentally validate our overall ranking results, we focused on two key markers involved in pyroptosis activation: ASC and Casp-1. Pyroptosis is distinct from the well-known processes of apoptosis and necrosis40. Unlike apoptosis, which is generally non-inflammatory, and necrosis, characterized by uncontrolled cell death, pyroptosis is uniquely driven by inflammation32,40. This process is mediated by the inflammasome, a protein complex central to the innate immune system’s response to cellular stress and infection. The inflammasome comprises key components such as ASC, NOD-like receptor pyrin-containing proteins (NLRP), Caspase-1, and gasdermin-D41,42.

In previous studies, we observed that patients undergoing long-term glaucoma therapy exhibited significantly elevated Caspase-1 protein levels compared to both DED patients and healthy controls. Specifically, the glaucoma group had Caspase-1 levels averaging 109.20 ± 42.59 pg/mL, whereas DED group averaged 91.62 ± 43.86 pg/mL, and healthy controls averaged 54.88 ± 23.04 pg/mL. These differences were statistically significant, with p-values of 0.001 and 0.003 when comparing glaucoma and DED groups to controls, respectively14. These findings highlight the potential of inflammasome markers in elucidating inflammatory mechanisms in these pathologies and underscore possibly their value in developing new diagnostic biomarkers.

To evaluate the expression levels of ASC and Casp1 genes in the tears of control, DED, and glaucoma patients using BAK-preserved topical treatments, we utilized both the most stable reference genes (18S, RER1, and ACTB) and the least stable reference genes (UBC and AP3D1) identified in our validation analysis for normalization purposes. Normalization using the top-ranked reference genes revealed a statistically significant upregulation of both ASC and Casp1 in the DED group (2.47 ± 0.91, p = 0.0001 and 2.56 ± 0.62, p < 0.0001, respectively) and in the glaucoma group (2.384 ± 0.83, p = 0.0004 and 2.020 ± 0.71, p = 0.0103, respectively) compared to the control group (ASC: 1.016 ± 0.20 and Casp1: 1.019 ± 0.20) (Fig. 7A). Conversely, normalization with the combination of the two least stable reference genes did not reveal a significant expression trend, with the DED group showing p-values of 0.1283 for ASC and 0.0886 for Casp1, and the glaucoma group showing p-values of 0.5235 for ASC and 0.9276 for Casp1 (Fig. 7B). These findings underscore the importance of selecting appropriate reference genes, as the use of stable reference genes enhances the detection of relevant gene expression regulations.

Fig. 7
figure 7

Expression levels of ASC and Casp-1 genes in tears of control, DED, and glaucoma groups. (A) Gene expression was normalized using the geometric mean of the Cq values of the three most stable reference genes: 18S, RER1, and ACTB. (B) Gene expression was normalized using the geometric mean of the Cq values of the two least stable reference genes: UBC and AP3D1. Bars indicate mean ± SD. Statistical significance was assessed via the One-way ANOVA test, with * p < 0.05, ** p < 0.01, **p < 0.001, and **** p < 0.0001 indicating comparisons to the control group. Control (n = 8), DED (n = 8), and glaucoma (n = 8) tear samples.

Discussion

Our study aimed to identify and validate stable reference genes suitable for normalizing RNA expression in tear samples from individuals with DED and glaucoma patients using BAK-preserved topical treatments and healthy controls. Through multiple analytical algorithms, namely the comparative ∆CT method, BestKeeper, NormFinder, geNorm, and RefFinder, 18S, RER1, and ACTB consistently emerged as the top candidates, exhibiting high expression stability. Conversely, UBC and AP3D1 demonstrated comparatively higher variability and proved less suitable for normalization. Such findings underscore the necessity of systematically validating reference genes when performing transcriptomic analyses on tear fluid, a medium known for its limited volume and complex composition.

Our analysis revealed that 18S, RER1, and ACTB exhibited consistent stability across various algorithms. The 18S ribosomal RNA gene is essential for ribosome assembly and protein synthesis, and it has been extensively documented as stable across different tissues and conditions28,43. Similarly, RER1, a receptor involved in the retrieval of endoplasmic reticulum (ER) membrane proteins, showed high stability, likely due to its vital function in preserving ER homeostasis44. ACTB, which encodes beta-actin, is a recognized housekeeping gene; its structural role in the cytoskeleton typically ensures stable expression levels34,45.

In contrast, UBC and AP3D1 showed greater variability. UBC encodes the polyubiquitin precursor, a key player in the ubiquitin–proteasome system that controls protein breakdown and turnover. The expression of UBC is naturally dynamic since the ubiquitination process reacts significantly to cellular stress, inflammation, and metabolic shifts46. Conditions affecting the ocular surface, like DED and patients treated with BAK, frequently experience varying stress and inflammatory signals in the cellular environment13,14,15. These fluctuations can modify ubiquitination, which in turn affects UBC expression levels. As a result, UBC may not provide the stable expression necessary for consistent normalization in tear fluid RNA analysis.

Similarly, AP3D1, a subunit of the adaptor protein complex involved in vesicle-mediated transport, was found to be less stable. The dynamic secretory processes may influence AP3D1’s expression in ocular surface cells. The tear film is continuously replenished and modulated by vesicular transport mechanisms that can be affected by both basal secretion and reflex tearing47, especially in pathological states such as DED or under the influence of toxic preservatives like BAK. The variability observed in AP3D1 may reflect fluctuations in vesicle trafficking activity, which is more pronounced during inflammatory or stress responses48,49. Notably, AP3D1 ranked highly stable in the BestKeeper algorithm, yet consistently emerged as one of the least stable genes in NormFinder and geNorm analyses. This discrepancy likely arises from fundamental methodological differences: while BestKeeper evaluates stability using the standard deviation of raw Cq values, it does not account for group stratification or biological variance, which are critical components of NormFinder and geNorm. Thus, the divergent ranking of AP3D1 may reflect limitations inherent to BestKeeper’s reliance on ungrouped Ct variation, rather than true biological stability. These findings underscore the importance of employing multiple algorithms in reference gene validation to avoid misleading conclusions based on a single computational metric.

Existing literature has long recognized the importance of rigorous reference gene validation, especially in clinically relevant samples such as ocular tissues and fluids3,4. Our finding that GAPDH exhibits varying stability aligns with previous reports showing that even reliable housekeeping genes may not consistently express uniformly in certain pathological scenarios26,27,28,50. The inconsistencies observed in AP3D1, which was considered highly stable by BestKeeper but received low rankings in other algorithms, mirror similar discrepancies reported in other tissue studies, further demonstrating that no single universal reference gene exists43,51.

We hypothesized that at least one canonical housekeeping gene, such as ACTB or GAPDH, would emerge as a robust reference. Our comprehensive approach indeed placed ACTB among the top-performing genes, corroborating part of our hypothesis. However, GAPDH exhibited moderate stability rather than top-tier performance, partially refuting our initial expectation. Moreover, we anticipated that genes like UBC might be less stable due to their involvement in dynamic processes such as ubiquitination; this presumption was confirmed, as UBC consistently ranked low across multiple methods.

Our findings support research indicating that depending on just one reference gene may result in biased or inaccurate data interpretations7,8. Multiple lines of investigation have similarly documented that 18S often exhibits robust stability in ocular tissues, highlighting its broad utility52,53. Additionally, the minimal variability shown by RER1 agrees with separate findings that link endoplasmic reticulum-related genes to stable expression profiles under various stress conditions54,55. These similarities enhance the dependability of our multi-algorithm method for reliably identifying stable reference genes.

A key aspect of our study is the use of five distinct algorithms: geNorm, ∆CT method, NormFinder, BestKeeper, and RefFinder, to evaluate reference genes performance. Each tool targets different facets of expression variability, from pairwise gene comparisons (geNorm) to estimates of intra- and inter-group variance (NormFinder)27,28,37,38,39. Utilizing RefFinder to integrate outcomes has provided us with an in-depth understanding of gene stability. This multifaceted approach reduces potential biases associated with single-algorithm techniques, confirming that the identified reference genes are reliably strong across various analytical systems.

The proven stability of 18S, RER1, and ACTB in tears has significant implications for discovering biomarkers and clinical diagnostics. As interest in tear-based transcriptomics for early disease detection rises, dependable normalization becomes crucial. When target gene expression levels accurately reflect real biological changes instead of technical errors, clinicians and researchers can interpret variations in inflammatory mediators with greater confidence, opening up opportunities for earlier interventions in conditions such as DED and glaucoma patients using BAK-preserved topical medications. To further elucidate the clinical relevance of our findings, the validated reference genes enable precise normalization of inflammasome-related genes like ASC and Caspase-1, which showed significant upregulation in DED and glaucoma patients. This upregulation is consistent with prior protein-level findings14 underscores the potential of tear RNA to detect early inflammatory changes, particularly in BAK-related toxicity.

However, the translational potential of tear RNA as a biomarker depends on addressing its correlation with tissue-specific expression and methodological challenges. Tear RNA, derived primarily from exfoliated conjunctival and corneal epithelial cells, infiltrating immune cells, and lacrimal gland exosomes, offers a non-invasive snapshot of ocular surface health but may not fully reflect localized transcriptomes due to its heterogeneous cellular origins1,56. For instance, conjunctival or corneal biopsies may reveal specific inflammatory or stress-response gene signatures that differ from tear RNA profiles, which aggregate contributions from multiple ocular sources and are modulated by tear film dynamics, as observed in comparative tear-tissue studies5. This complexity poses a translational challenge, as direct correlations with tissue-specific expression remain understudied. Future investigations incorporating paired tissue biopsies and tear samples from the same subjects, alongside spatial transcriptomics or in situ hybridization, are warranted to validate the extent to which tear-based expression mirrors intra-tissue transcriptional dynamics.

Additionally, the reliability of tear RNA measurements is highly susceptible to methodological heterogeneity in sample collection. Schirmer strips, which often induce reflex tearing and enrich epithelial cell content, contrast with microcapillary tubes, which selectively capture basal tears but yield limited RNA, and the Weck-Cel® spears, which balance sample volume with minimal invasiveness1,5,56. This variability in collection techniques can alter RNA profiles by differentially sampling cellular and exosomal components, complicating biomarker consistency. Clinical factors, including disease subtype heterogeneity, severity, and medication use, further confound analyses by introducing patient-specific transcriptional variability. To address these hurdles, future studies should leverage single-cell RNA sequencing or spatial transcriptomics to correlate tear RNA with conjunctival and corneal biopsies, validating biomarker specificity. Establishing standardized collection protocols will be critical to ensure reproducibility across diverse clinical cohorts, enhancing the diagnostic potential of tear-based transcriptomics.

Tear fluid mRNA primarily originates from exfoliated conjunctival and corneal epithelial cells, with contributions from infiltrating immune cells and extracellular vesicles such as exosomes, contrasting with the homogeneous cellular RNA in cell culture17,56. This heterogeneity influences reference gene selection, as housekeeping genes validated in cell culture (e.g., 18S, ACTB, GAPDH) may not always be stable in extracellular fluids due to post-transcriptional regulation or selective exosomal packaging. Recent studies, such as Boychev et al.18, demonstrate that ACTB and GAPDH are stable in tear fluid collected via contact lenses or Schirmer strips in a rabbit model, with GAPDH more stable than ACTB. In our human study, ACTB’s high stability aligns with these findings, though GAPDH’s moderate stability (geNorm M = 0.644) may reflect microenvironmental changes from inflammation or BAK exposure in DED and glaucoma. Unlike cell culture’s controlled conditions, tear fluid’s dynamic microenvironment may modulate gene expression, yet the consistent stability of 18S, RER1, and ACTB supports their suitability for tear fluid mRNA normalization, particularly for mRNA targets like ASC and Caspase-1 (Fig. 7). Exosomal mRNA, likely a minor component in our samples, may require alternative reference genes for exosome-specific analyses. Future studies using paired conjunctival cells and tear samples, or single-vesicle RNA sequencing, should clarify the relative contributions of exosomal versus cellular mRNA.

Although our results suggest the most stable reference genes in tear fluid from DED and glaucoma patients using BAK-preserved topical treatments, several limitations should be acknowledged. Firstly, our sample size of 16 participants (8 glaucoma and 8 controls) restricts the generalizability of these results. Larger and more diverse cohorts are likely to yield more profound insights regarding the variations in reference gene expression that may occur in relation to different disease severities and demographic factors. Second, our study focused exclusively on seven genes identified through the Human Eye Transcriptome Atlas Project33, and additional candidates might exist that yield even more robust normalization. Finally, we centered on tear samples and did not investigate parallel ocular tissues such as the cornea or conjunctiva, thus restricting the scope of inference regarding gene stability in other microenvironments of the eye.

In conclusion, our study establishes 18S, RER1, and ACTB as the most reliable reference genes for normalizing RT-qPCR data from tear fluid samples in individuals with DED and glaucoma patients using BAK-preserved topical medications, as well as healthy controls. Integrating various analytic algorithms enhances the credibility of our findings and establishes a framework for future ocular biomarkers research. Continual validation across different patient populations, disease states, and ocular tissues remains essential to refining normalization practices. By employing these methodologies, we can utilize the diagnostic capabilities of tear fluid more effectively for the early and precise identification of ocular surface inflammation, thereby directing focused therapeutic strategies for conditions such as DED and ocular toxicity associated with therapies preserved by BAK.

Methods

Study population

This study was approved by the Institutional Review Board of the University of Miami Miller School of Medicine (protocol #20190334) and was conducted in compliance with the United States Health Insurance Portability and Accountability Act and the Declaration of Helsinki principles. Patients were evaluated at the Bascom Palmer Eye Institute and recruited from either the Dry Eye Clinic or the Glaucoma Division (for the case groups) or from optometry clinics (for the control group). For the DED group, participants were included if they exhibited symptoms or clinical signs of dry eye, specifically an Ocular Surface Disease Index (OSDI) score of 13 or higher and/or corneal staining (CS) scores of 3 or greater, as graded using the National Eye Institute (NEI) scale. The glaucoma group consisted of individuals undergoing treatment with BAK-containing topical hypotensive medications, regardless of whether they exhibited ocular surface symptoms. The control group included asymptomatic individuals with OSDI scores below 13 and no signs of ocular surface damage (CS < 3). Exclusion criteria for all groups included a history of ocular radiotherapy, pregnancy, age outside the range of 21 to 90 years, allergic diseases, active infections, and any previous ocular surgery. In total, 16 patients (8 individuals in the DED group, 8 in the glaucoma group) and 8 controls were enrolled in the study (Table 1). Detailed clinical characteristics of the DED and glaucoma groups are provided in Supplementary Tables 1 and 2, respectively.

The sample size was primarily determined based on the availability of clinical samples and precedent from similar ocular biomarker studies2,5. A retrospective power analysis was also conducted using effect size data from comparable published datasets using G*Power 3.1.9.7 software57. Effect size estimates were derived from previously published Caspase-1 tear expression data in comparable patient populations14. The analysis indicated that, for a control versus glaucoma comparison (Cohen’s d = 1.47), a minimum of 9 participants per group would be required to detect statistically significant differences with 80% power at α = 0.05. Although our groups included 8 subjects each, this number closely approximates the target and is considered adequate for exploratory gene expression analysis in this pilot study.

Data and tear sample collection

After obtaining informed consent, participants first completed the OSDI questionnaire, followed by the collection of tear samples using Weck-Cel® eye spears (BVI Medical). Although previous research indicates minimal diurnal fluctuations in tear composition58,59. All samples in this study were collected within a standardized time window (9:00 AM–2:00 PM) to further minimize potential variability. For glaucoma patients, tears were collected at least two hours after the last BAK-preserved medication dose to reduce acute effects, with confirmed use on the day of and the day before sampling. The sampling procedure was conducted by trained personnel, who gently placed the spear in the inferior lateral tear meniscus for 10 s, minimizing contact with the ocular surface and lid margin. No topical anesthetic was administered prior to sampling. Immediately, after tear collection, each participant underwent a detailed clinical evaluation. The ocular surface was assessed with the Oculus Keratograph 5M Topographer (Oculus Optikgeräte GmbH, Wetzlar, Germany). Important parameters measured included non-invasive tear break-up time (NIBUT) and bulbar redness (graded using the Jenvis scale). Topical fluorescein was applied, and corneal staining (CS) was evaluated using the NEI grading scale60.

Sample processing

Tear samples were collected from one eye using the Weck-Cel® Eye Spear for 10 s, ensuring a non-traumatic collection from the lateral canthus to reduce reflex tearing. For glaucoma patients using BAK-preserved topical treatments, samples were obtained from the affected eye if treatment was solely unilateral. If bilateral, the use of medication the day of and the day prior to sampling was confirmed. The samples were then placed into separate sterile 1.5 ml collection tubes containing 300 µl of PureLink lysis buffer (PureLink RNA Mini Kit, Invitrogen). They were kept cool during collection and were stored at − 80 °C until processed. The cellulose composition of the Weck-Cel® spears facilitates rapid and effective collection of tear fluid.

RNA extraction

We transferred Weck-Cel® spears containing the tear sample and PureLink lysis buffer to a new collection tube equipped with a homemade spin column constructed by piercing the bottom of a sterile 0.5 mL Eppendorf tube with a syringe needle and placing it within a 1.5 mL collection tube. This configuration allowed the Weck-Cel® sponge to be retained while permitting efficient flow-through of the lysis buffer during centrifugation. The use of lysis buffer for RNA extraction from tear film has been described in previous studies17, and it facilitates efficient recovery of RNA from the tear fluid, which primarily contains extracellular RNA with only minimal contributions from exfoliated ocular surface cells. Although lysis buffer is often associated with extracting cellular RNA, our method is optimized for tear fluid, and we did not perform a direct cell count; thus, the exact cellular contribution remains unquantified. The mixture was centrifuged at 12,000×g for 5 min at room temperature.

Following centrifugation, the supernatant containing the sample and buffer was recovered for RNA extraction using the PureLink® RNA Mini Kit (Invitrogen) in accordance with the manufacturer’s instructions. An RNase-Free DNase kit (Qiagen) was employed during the extraction process to eliminate potential DNA contamination. The RNA was then eluted in 20 µL of RNase-free water and stored at − 80 °C for future analysis. We assessed the quantity of total RNA by measuring absorbance at 260 nm using a Nanodrop ND-1000 UV–Vis Spectrophotometer, with purity confirmed by checking the 260/280 nm and 260/230 nm ratios. RNA concentration and purity data for each sample are provided in Supplementary Table 3. After that, 50 ng of purified RNA was reverse transcribed into complementary DNA (cDNA) using the iScript Advanced cDNA Synthesis Kit (Bio-Rad Laboratories), following the manufacturer’s guidelines.

Reference gene selection and stability assessment

Accurate normalization of gene expression data is crucial in quantitative real-time PCR (RT-qPCR) studies to ensure reliable and reproducible results. The selection of appropriate reference genes, also known as housekeeping genes, is critical for mitigating technical variability and accurately quantifying target gene expression. In the context of ocular research, particularly when analyzing RNA from tear fluid, identifying stable reference genes is essential due to the limited volume and complex composition of tear samples.

Six candidate reference genes (RER1, ACTB, GAPDH, PGK1, UBC, AP3D1) were selected based on their documented high expression in conjunctiva, cornea, and eyelid tissues (the primary cellular sources of tear fluid RNA) in the Human Eye Transcriptome Atlas Project (https://www.eye-transcriptome.com)33, and their established or potential use in ocular surface transcriptomic studies34. Additionally, biological functions were considered to ensure relevance to the cellular and inflammatory contexts of DED and BAK-treated glaucoma.18S (Assay ID: Hs03003631_g1) was chosen for its high abundance in cellular RNA, constituting approximately 20% of total cellular RNA content across conjunctiva-derived, limbal-derived, and cultured ocular surface epithelial cells34,45. ACTB (Assay ID: Hs01060665_g1) and GAPDH (Assay ID: Hs02786624_g1) were selected for their widespread use as housekeeping genes in RT-qPCR normalization. PGK1 (Assay ID: Hs00943178_g1) was included for its role in stable metabolic processes, and RER1 (Assay ID: Hs00199824_m1) for its function in endoplasmic reticulum homeostasis, hypothesized to be stable in ocular surface contexts. UBC (Assay ID: Hs05002522_g1) and AP3D1 (Assay ID: Hs00926919_m1) were chosen to evaluate their stability in the inflammatory contexts of DED and BAK-treated glaucoma, despite their roles in dynamic processes like ubiquitination and vesicular transport, respectively. Table 3 provides details about the seven reference genes, including their GenBank accession numbers, assay identification numbers, and specific functions.

Table 3 List of selected candidate reference genes analyzed in the tear samples.

The expression stability of these seven candidate reference genes was evaluated using widely recognized methods: geNorm, ΔCt method, NormFinder, BestKeeper and RefFinder. Each method employs a distinct algorithm to assess gene stability, providing a comprehensive analysis of potential reference genes.

geNorm: geNorm utilizes the 2^ − ΔCt values as input and calculates an average pairwise variation (M value) for each gene, with lower M values indicating higher stability28. Furthermore, geNorm can determine the optimal number of reference genes required for accurate normalization by measuring pairwise variation between genes. In this study, geNorm was accessed through both qbase + (https://cellcarta.com/genomic-data-analysis/) and the online tool RefFinder.

ΔCt method: The ∆CT method was employed to evaluate the relative expression of gene pairs within each sample, facilitating the identification of suitable reference genes for expression analyses37. This approach involves calculating the ∆CT values, defined as the differences in quantification cycle (Cq) values between pairs of reference genes, alongside their corresponding Standard Deviation (SD) values. Minimal variability in ∆CT indicates stable expression between gene pairs, ensuring reliable normalization across diverse samples and experimental conditions. Conversely, significant fluctuations in ∆CT suggest variability in the expression of at least one gene within the pair, undermining normalization accuracy.

NormFinder: It analyzes raw Cq values, performing a grouped assessment of potential reference genes38. By combining intra- and inter-group variation, NormFinder assigns a stability value (S value), where smaller S values indicate more stable gene expression.

BestKeeper: This method also relies on raw Cq values and assumes a PCR efficiency of 2, linking gene stability to lower standard deviation (SD) and coefficient of variation (CV)27. BestKeeper indicates that a reference gene with an SD value over 1 is typically regarded as unstable and not suitable for normalization. To further confirm the stability of the chosen reference genes and to pinpoint the most dependable ones for normalization.

RefFinder (accessible at https://www.ciidirsinaloa.com.mx/RefFinder-master/): It is a robust online tool that combines results from the comparative ∆Ct method, NormFinder, geNorm, and BestKeeper to provide an overall stability ranking39. It computes a geometric mean for each candidate reference gene, where smaller values indicate greater stability. The final selection of reference genes in this study was based on the consensus across all software, improving normalization accuracy across different experimental conditions. Graphical representations of gene stability rankings and pairwise variation analyses were generated using GraphPad Prism 10.0, facilitating clear visualizations of each gene’s stability profile.

Quantitative real-time polymerase chain reaction (qRT-PCR) analysis

All qRT-PCR reactions utilized SsoAdvanced Universal Probes Supermix (Bio-Rad) alongside gene-specific Taqman assays. Each sample was conducted in duplicates in a 20 µL reaction volume on the Azure Cielo 3 Real-Time PCR System (Azure Biosystems). The PCR protocol commenced with an initial step of 2 min at 95 °C, followed by 40 cycles of denaturation at 95 °C for 10 s and annealing/extension at 60 °C for 30 s. An automatic threshold was set for each assay, using water as a negative control. The relative expression levels of target genes were calculated through the 2 − ΔΔCt method, utilizing the three most stable reference genes (18S, RER1, and ACTB), where ΔCt indicates the difference between the Cq value of the target gene and that of the reference gene. The quantification cycle (Cq) value for stable reference genes is obtained from the geometric mean of their individual Cq values. For gene expression analysis, we concentrated on the mRNA levels of the ASC and Casp-1 genes, employing TaqMan gene expression assays (ASC: Hs01547324_gH and Caspase1: Hs00354836_m1), as both play a role in inflammasome activation.

Statistical analysis

The Cq values from RT-qPCR were calculated for each sample, which included DED and glaucoma patients using BAK-preserved topical treatments, as well as controls, based on the average of duplicate technical replicates. While we did not conduct a formal power analysis prior to the study, our sample size was determined based on the availability of clinical samples and precedent in similar ocular biomarker studies2,5. Prior to analysis, data normality was verified using the Shapiro‑Wilk test. Ocular surface parameters and target gene expression were then compared among groups using a one‑way ANOVA, followed by Tukey’s multiple‑comparison test. All statistical analyses mentioned were performed using GraphPad Prism 10 software. The results are presented as mean ± standard deviation (SD), with a p value of < 0.05 considered statistically significant.