text:
The JAK/STAT Pathway and Its Selective Inhibition in the Treatment of Atopic Dermatitis: A Systematic Review

In recent years, the broadening understanding of the pathogenesis of atopic dermatitis (AD) has led to the development of novel therapeutic molecules, that target core inflammatory components of the disease. The Janus kinase (JAK)/signal transducer and activation of transcription (STAT) pathway constitutes the principal signaling cascade for a large number of cytokines and growth factors and is involved in intracellular signal transduction and subsequent regulation of gene transcription. Current knowledge suggests that the robust activation of the T-helper (Th)-2 [interleukin (IL)-4, IL-5, IL-13, IL-31] and Th22 (IL-22) immune responses in both skin and serum plays a pivotal role in the immunopathogenesis of AD especially at the acute stage, followed by a variable degree of Th1 (interferon-γ, tumor necrosis factor alpha) and Th17 (IL-17) activation in chronic disease. Of note, most of the aforementioned inflammatory cytokines utilize the JAK/STAT pathway for downstream signal transduction, explaining the emerging role of JAK inhibitors in the therapeutic armamentarium of AD. The present systematic review aims to discuss the involvement of JAK/STAT pathway in the pathogenesis of AD and summarize the clinical data available on the efficacy and safety of JAK inhibitors which have been used in the treatment of AD thus far.

Introduction
Atopic dermatitis (AD) is a chronic inflammatory skin disease with multifactorial pathophysiology, involving strong genetic predisposition, skin barrier dysfunction, microbiome abnormalities, immune dysregulation and extrinsic factors [1]. Current knowledge suggests that robust activation of the T-helper (Th)-2 [interleukin (IL)-4, IL-5, IL-13, IL -31] and Th22 (IL-22) immune responses in both skin and serum play a pivotal role in the immunopathogenesis of AD especially at the acute stage, followed by a variable degree of Th1[interferon-γ (IFN-γ), tumor necrosis factor alpha (TNF-a)] and Th17 (IL-17) activation in chronic disease [1,2]. Moreover, it is known that most of the aforementioned inflammatory cytokines, including the two hallmark Th2 cytokines IL-4 and IL-13, utilize the Janus kinase (JAK)/signal transducer and activation of transcription (STAT) pathway for downstream, intracellular signal transduction [3]. Apart from immune regulation, the JAK/STAT pathway is involved in multiple biological processes, such as hematopoiesis, cell proliferation, differentiation and apoptosis [4].
The JAK family comprises four different cytoplasmatic tyrosine kinases, the JAK1, JAK2, JAK3 and tyrosine kinase-2 (TYK2), which are linked to the intracellular domains of multiple transmembrane cytokine receptors. Their function is coupled to the STAT family, which comprises seven intracellular transcription factors, the STAT1, STAT2, STAT3, STAT4, STAT5a, STAT5b, and STAT6. These proteins have an N-terminal domain, a DNA-binding domain, and a C-terminal transactivation domain, and contribute to cellular immunity, proliferation, apoptosis and differentiation [5]. Binding of the signaling cytokines to their receptors leads to the activation of the associated JAKs and subsequent phosphorylation of STATs, which, in turn, separate from the receptor, dimerize and translocate to the cell nucleus to induce gene transcription after binding to specific DNA sequences [5]. Of note, more than 50 cytokines signals via the JAK/STAT pathway have been identified, since each JAK may be associated with multiple types of cytokine receptors and, thus, exhibit a different downstream effect on cytokine signaling. Similarly, members of the STAT family may interact with different JAKs, with STAT1 and STAT4 being involved in Th1 and STAT3 and STAT6 in Th2 and Th17 immune responses, respectively [6][7][8].
The broadening understanding of AD pathogenesis in recent years has led to the development of novel therapeutic molecules, that target core inflammatory components of the disease, including the JAK inhibitors. The present review aims to discuss the involvement of the JAK/STAT pathway in the pathogenesis of AD as well as the emerging role of JAK inhibitors in the therapeutic armamentarium of AD.

Materials and Methods
The MEDLINE electronic database was searched systematically via PubMed for randomized controlled trials (RCTs) assessing the efficacy and/or safety of JAK inhibitors in AD from inception to May 2022. The search strategy is summarized in Table 1. Nonhuman studies and articles reporting subgroup and post hoc analyses of data extracted from the pivotal studies were excluded. The initial search yielded 417 articles (no filters used). After title and abstract review (conducted by two independent reviewers, A.T., M.K), 40 articles were selected for full text review. Discrepancies between the two authors that performed the quality assessment of the included studies, were solved by the contribution of a third author (S.G.). Finally, 27 studies were included in the systematic review after the exclusion of irrelevant studies ( Table 2, Scheme 1). Reference lists of the included articles as well as the database of clinical trials of the US national library of medicine (https://www.clinicaltrials.gov, accessed on 10 June 2022) were also screened in order to detect any additional RCTs ( Table 2). The following data were extracted: JAK inhibitor used, study population, treatment period, disease assessment scores and safety reports. This systematic review was performed according to the Preferred Reporting Items for Systematic reviews and Meta-Analyses (PRISMA). Apart from the systematic review on the efficacy and safety of JAK inhibitors in AD, a narrative review regarding the involvement of the JAK/STAT pathway in AD was performed and is addressed below. This systematic review was performed according to the Preferred Reporting Items for Systematic reviews and Meta-Analyses (PRISMA). Apart from the systematic review on the efficacy and safety of JAK inhibitors in AD, a narrative review regarding the involvement of the JAK/STAT pathway in AD was performed and is addressed below.

Overregulation of Th2 and Th22 Immune Responses
The involvement of Th2 immune response in AD pathogenesis is well established and is associated with increased signaling through all four JAKs. Mouse models support that IL-4 is required for the production of Th2-related cytokines, including IL-4, IL-5, IL-13 and IL-31, which play a critical role in the development of AD lesions [36]. There are two types of IL-4 receptors; the type 1 receptor is comprised of an IL-4Rα and a γ chain and the type 2 receptor of an IL-4Rα and an IL-13Ra1 chain. The latter functions also as an IL-13 receptor. Cytokine binding to type 1 receptors leads to the phosphorylation of JAK1 and JAK3, resulting in activation of STAT6, while binding to type 2 receptors induces JAK1 and TYK2 expression and subsequent activation of STAT6 and STAT3 [36]. The STAT6 protein has been found to regulate T-cell proliferation and Th2 differentiation through the upregulated expression of the Th2-specific transcription factor GATA3, but also the promotion of immunoglobulin (Ig) class switching to IgE and IgG1 in B cells [37]. For that reason, various STAT6 polymorphisms have been related to high levels of IgE and greater susceptibility to AD [38]. The Interleukin-5-mediated pathway utilizes JAK1 and JAK2 as well as STAT1, STAT2 and STAT5 to transduce its signal [39][40][41]. Finally, the intracellular signaling of IL-31, as well as of thymic stromal lymphopoietin (TSLP), a highly pruritogenic protein, which also promotes the Th2 differentiation, is transmitted by JAK1 and JAK2, with subsequent involvement of STAT1, STAT3 and STAT5 [42,43].
The JAK/STAT-dependent IL-4 and IL-13 signaling is also very important in the dysregulation of keratinocyte function in AD, due to the downregulation of several barrier proteins [42]. For example, IL-4 and IL-13 downregulate loricrin and involucrin in keratinocytes, while they have been also found to reduce the expression of filaggrin in patients who did not carry the known filaggrin mutations [44]. Involucrin and loricrin are structural proteins that contribute to the protective barrier function of the stratum corneum. The superficial spinous cells express involucrin which cross-links to form the cornified envelope scaffolding. Loricrin in turn is cross-linked to involucrin and is expressed in the superficial granular cells forming composite keratohyalin granules or L granules along with profilaggrin [45].
Moreover, downstream signaling of IL-4 and IL-13 modulates the expression of innate immune response genes, including cathelicidin and β-defensins, enhancing the susceptibility to skin infections, for example with staphylococcus aureus, which, in turn, may exacerbate AD by acting on keratinocytes and immune cells [46][47][48].
IL-4 and IL-13 reduce the expression of filaggrin, loricrin, involucrin and lipid components of the skin barrier [45]. This barrier dysfunction contributes to Staphylococcus aureus colonization and subsequent inflammation. In addition, activation of the IL-13/IL-4-JAK-STAT6/STAT3 axis inhibits the aryl-hydrocarbon receptor mediated up-regulation of filaggrin, loricrin, and involucrin and also inhibits the cytoplasmic-to-nuclear translocation of the transcription factor Ovo-like-1 (OVOL1) that promotes the expression of filaggrin and loricrin [49]. Moreover, it induces keratinocytes to produce periostin that subsequently upregulates the IL-24 expression and IL-24 inhibits the expression of filaggrin via STAT3 activation [49]. The Th22 immune response is also implicated in AD pathogenesis and is particularly increased in chronic AD lesions [2]. IL-22 promotes skin modification and it has been identified as a key mediator of epidermal hyperplasia. The binding of IL-22 to its receptor leads to the phosphorylation of JAK1 and TYK2, activating STAT1/STAT3/STAT5 [50,51]. Figure 1 illustrates the involvement of the JAK/STAT signaling pathway in AD. of filaggrin, loricrin, and involucrin and also inhibits the cytoplasmic-to-nuclear translocation of the transcription factor Ovo-like-1 (OVOL1) that promotes the expression of filaggrin and loricrin [49]. Moreover, it induces keratinocytes to produce periostin that subsequently upregulates the IL-24 expression and IL-24 inhibits the expression of filaggrin via STAT3 activation [49]. The Th22 immune response is also implicated in AD pathogenesis and is particularly increased in chronic AD lesions [2]. IL-22 promotes skin modification and it has been identified as a key mediator of epidermal hyperplasia. The binding of IL-22 to its receptor leads to the phosphorylation of JAK1 and TYK2, activating STAT1/STAT3/STAT5 [50,51]. Figure 1 illustrates the involvement of the JAK/STAT signaling pathway in AD.

Regulation of Th1 and Th17 Immune Responses
Although the main cytokines that are involved in AD pathogenesis are the Th2 cytokines, Th1 and Th17 immune responses likely play a pathogenetic role in certain AD subtypes and in chronic AD lesions [52]. The hallmark cytokine of Th1 is IFN-γ, a protein that utilizes STAT1 for signal transduction. Moreover, STAT4 is necessary for IL-12 signaling, which directs the differentiation of naïve T helper cells into Th1 cells [53,54].
Th17 cytokines are elevated in blood and lesional AD and may exaggerate skin both acute and chronic AD [55]. After binding of IL17 to its receptor, the adaptor protein Act1 interacts with the TNF-receptor-associated factor 6 and the transforming growth factor-β-

Regulation of Th1 and Th17 Immune Responses
Although the main cytokines that are involved in AD pathogenesis are the Th2 cytokines, Th1 and Th17 immune responses likely play a pathogenetic role in certain AD subtypes and in chronic AD lesions [52]. The hallmark cytokine of Th1 is IFN-γ, a protein that utilizes STAT1 for signal transduction. Moreover, STAT4 is necessary for IL-12 signaling, which directs the differentiation of naïve T helper cells into Th1 cells [53,54].
Th17 cytokines are elevated in blood and lesional AD and may exaggerate skin both acute and chronic AD [55]. After binding of IL17 to its receptor, the adaptor protein Act1 interacts with the TNF-receptor-associated factor 6 and the transforming growth factorβ-activated kinase 1, leading to activation of nuclear factor (NF)-kB and p38/mitogenactivated protein kinases. This implies that IL17-mediated inflammation does not depend on JAK/STAT signaling [50]. However, literature data suggest that STAT3 is necessary for the induction of retinoic acid receptor-related orphan receptor γ (RORγt), a Th17-specific transcription factor that induces IL17 expression [53].

Suppression of Regulatory T Cells
Regulatory T cells (Tregs) in the skin are known to regulate the immune response in order to maintain immune tolerance to self-antigens and prevent excessive responses to foreign antigens. In AD, they play an important role in controlling skin inflammation and its expansion to the lungs according to a murine AD model, where the depletion of Tregs has led to disease exacerbation, due to Th2 expansion and IgE elevation [56]. Differentiation, proliferation and maintenance of Tregs depend on the transcription factor Foxp3 and mutations leading to altered expression or function of Foxp3 may lead to a severe autoimmune disease called immunodysregulation polyendocrinopathy enteropathy X-linked syndrome [56]. In a mouse model, Foxp3 has been found to have a short halflife and to be controlled by the JAK/STAT signaling [57]. Foxp3 has been found to be suppressed by GATA3, the main regulator of Th2 cells, which is upregulated by STAT6, as mentioned above. In addition, STAT5A/B and STAT3 may also be involved in the downregulation of Foxp3, all explaining why Treg function is suppressed in AD [42,58].

Stimulation of Eosinophils and Mast Cells
Interleukin-31 has been found to delay the apoptosis of eosinophils and significantly induce the secretion of pro-inflammatory cytokines and chemokines after binding on the IL-31 receptor A and subsequent activation of STAT3 in eosinophils [59]. Moreover, upregulation of IL-5 has been shown to induce eosinophilia and regulate the proliferation, survival and functions of eosinophils in a JAK2-STAT1/STAT5-dependent manner [41,60]. Interleukin-9 has been found to be increased in the lesional skin of AD patients and it has been associated with localized stimulation of VEGF from mast cells through activation of STAT3 in vitro [51,61]. Finally, the JAK/STAT pathway is involved in the maturation of B cells and their differentiation into plasma cells with the subsequent secretion of IgE, which activates cutaneous mast cells to release histamine [6,62].
Apart from the IgE-dependent mast cell activation, the JAK/STAT pathway is also involved in mast cell homeostasis and proliferation. The binding of the stem cell factor to c-kit on the human mast cell surface leads to JAK2 activation and subsequent STAT5 and STAT6 phosphorylation, which are all involved in mast cell development, survival and proliferation [63]. In addition, the same molecules (JAK2/STAT5/STAT6) may be activated by IL3 as well, a cytokine that regulates the mast cell response against antigens. Besides, the important role of the JAK/STAT pathway in mast cells is underlined by the consequences of the mutations of these molecules, including those that are found in mastocytosis [63].

JAK Inhibitors for the Treatment of AD
JAK inhibitors are small molecules, which are classified into first generation and newer JAK inhibitors. First generation JAK inhibitors bind to all three JAKs, whereas newer JAK inhibitors are more selective and demonstrate activity against specific JAKs [52]. JAK inhibitors can be administered orally due to their size and chemical structure. Some of them, including baricitinib and upadacitinib, are eliminated primarily by renal excretion and are characterized by a half-life between 9 and 14 h, making a single daily administration sufficient. On the other hand, JAK inhibitors, such as ruxolitinib and tofacitinib are metabolized in the liver and have a shorter half-life (approximately 3 h). Therefore, they should be administered twice daily [64]. All registered clinical trials on topical and systemic JAK inhibitors are listed in Table 2.

Topical JAK Inhibitors
To date, a number of trials assessed the efficacy and safety of several topical formulations in AD. Delgocitinib is a first-generation topical pan-JAK inhibitor, which was the first JAK inhibitor approved for AD in adults (0.5% ointment) and children (0.5% and 0.25% ointment) in Japan [65]. The two phase II clinical trials conducted in Japan showed the efficacy of delgocitinib in both adults and children (2-15 years of age) with a favorable safety profile [9,10]. In a double-blind, phase III RCT (n = 158, patients over 16 years of age with moderate-to-severe AD), delgocitinib 0.5% ointment met the primary endpoint, achieving improvement of percentage change from baseline modified EASI (%mEASI) at week 4 (44.3% vs. 1.7% for delgocitinib 0.5% vs. vehicle, respectively) [11]. Moreover, 95% of patients achieved at least 50% in the mEASI (mEASI-50) and 49% at least 75% in the mEASI (mEASI-75) during the 24-week extension phase. Regarding safety, the overall rate of treatment-related AEs in the delgocitinib group was 15.4%, with only eczema herpeticum occurring in three patients [11]. Similar positive results with a good safety profile were supported by a further phase III RCT with subsequent long-term extension (56 weeks) in pediatric patients between 2-15 years of age with AD [12].
Tofacitinib is another small molecule that inhibits selectively JAK1 and JAK3. In a double-blind phase IIa RCT evaluating tofacitinib 2% ointment in 69 adult patients with mild-to-moderate AD, there was an 81.7% improvement in Eczema Area and Severity Index (EASI), compared to 29.9% with the vehicle, as assessed at week 4. Tofacitinib also met several secondary endpoints, including improvement in Investigator's Global Assessment (IGA) score, Body Surface Area (BSA) and itch [13].
Ruxolitinib, a first-generation inhibitor of JAK1 and JAK2, was the first approved oral JAK inhibitor for the treatment of hematologic disorders, while it has also been investigated for the treatment of AD in the form of topical formulation. Following the positive results of a dose-ranging, double-blind, phase IIb RCT, the efficacy of ruxolitinib was evaluated in two parallel phase III trials [TRuE AD1 (n = 631), TRuE AD2 (n = 618)], where ruxolitinib 1.5% vs. 0.75% cream vs. vehicle were administered for eight weeks in adolescents and adults ≥ 12 years of age with mild-to-moderate AD [14]. Overall, the primary endpoint of IGA 0/1 with ≥2-grade improvement from baseline at week 8 was met by almost half of the patients in the active-agent groups. Moreover, both formulations met the secondary endpoints of EASI-75 and the critically relevant reduction of itch and exhibited rates of treatment-related AEs comparable to the vehicle (26.3% vs. 29.4% vs. 33.6%), with no serious AEs reported [14].

Baricitinib
Baricitinib is a selective JAK1 and JAK2 inhibitor, with less affinity for JAK3 and TYK2. The recommended dosage is 4 mg once daily but it may be reduced to 2 mg daily if the disease is well controlled in patients over 75 years of age, and in cases of renal impairment or risk of infections [70]. The results of a phase II, placebo-controlled RCT (n = 124) highlighted the efficacy of baricitinib combined with topical corticosteroids (TCS) in adults with moderate to severe AD, based on the observed achievement of EASI-50 at week 16 [15].
Subsequently, the efficacy and safety of baricitinib have been evaluated in two parallel, independent, multicentric, double-blinded, phase III RCTs [BREEZE-AD1 (n = 624) and BREEZE-AD2 (n = 615)] [16]. In both studies, adult patients with moderate to severe AD were randomized 2:1:1:1 to receive baricitinib 1 mg, 2 mg, 4 mg or placebo once daily for 16 weeks. In BREEZE-AD1, the primary endpoint of IGA 0/1 at week 16 was met by more patients on baricitinib arms compared with placebo [16.8% in baricitinib 4 mg, 11.4% in baricitinib 2 mg, 11.8% in baricitinib 1 mg, 4.8% in the placebo group, respectively]. Similarly, the same primary endpoint was met by 13.8% in baricitinib 4 mg, 10.6% in baricitinib 2 mg, 8.8% in baricitinib 1 mg, and 4.5% in the placebo group, respectively, in the BREEZE-AD2 study [16]. Moreover, baricitinib 4 mg achieved consistently all secondary endpoints in both trials, including EASI-75, improvement of skin pain and night awakenings and rapid relief of pruritus within the first week of therapy, suggesting the superiority of this dosing regimen compared to the others or placebo. After week 16, responders or partial responders from the baricitinib 4 mg and 2 mg groups (IGA ≤ 2) were further evaluated in a long-term extension study (BREEZE-AD3) and showed further improvement with a similar safety profile up to week 68, suggesting a sustained long-term benefit for these patients [17]. The efficacy of baricitinib monotherapy in adults with moderate-to-severe AD was evaluated in a further North American RCT [BREEZE-AD5 (n = 440)], showing similar efficacy and safety findings to those reported in the previous baricitinib AD studies [19].
Regarding the safety profile of baricitinib, treatment-related AEs were observed in 54% to 58% of patients, they were mild to moderate in severity and were not found to differ significantly between the four study arms in the two pivotal studies BREEZE-AD1 and BREEZE-AD2 [16]. The most frequently reported AEs were nasopharyngitis, headache and asymptomatic increased creatine phosphokinase (CPK) levels. Herpes simplex infection was observed in 7.2% vs. 3.3% vs. 5.5% vs. 1.2% of patients receiving baricitinib 4 mg vs. 2 mg vs. 1 mg vs. placebo in the BREEZE-AD1 study, but it was not so frequent in the BREEZE-AD2 study. Overall, no opportunistic infections, malignancies, thromboembolic events and major adverse cardiovascular events were reported in any study arm in both studies [16]. Pooled safety data of baricitinib were reported in a meta-analysis of six RCTs and two long-term extension studies (total n = 2531) highlighting a good overall safety profile [71]. The most common AEs were found to be upper respiratory tract and herpes simplex virus infections, exhibiting higher frequency in the baricitinib 4 mg group, while the most common serious AEs were eczema herpeticum [n = 11, adjusted incidence rate (IR) per 100 patient-years: 0.5], cellulitis (n = 6, IR: 0.3) and pneumonia (n = 3, IR: 0.1). Overall, two major adverse cardiovascular events and three thromboembolic events were identified in the total baricitinib data set [71].

Upadacitinib
Upadacitinib is a selective JAK1-inhibitor, which has been approved for the treatment of moderate to severe AD in adults and adolescents (≥12 years). The recommended dosage is 15-30 mg once daily based on disease severity, while a dosage of 15 mg once daily should be also recommended in patients ≤18 or ≥65 years of age [72]. In a multicenter, double-blinded, placebo-controlled phase IIb RCT (n = 164) patients were treated with varying doses of upadacitinib (7.5 mg, 15 mg, 30 mg) for 16 weeks [21]. The improvement of the EASI score was dose-dependent, while equal and rapid improvement of pruritus was observed in all active-treatment groups [21].
Two parallel, multicentric, double-blinded phase III RCTs [Measure Up 1 (n = 847), Measure Up 2 (n = 836)] evaluated the efficacy and safety of upadacitinib in moderate to severe AD in adults and adolescents (12-75 years) [22]. In both studies, patients were randomly assigned to receive either upadacitinib 30 mg or 15 mg or placebo once daily for 16 weeks. All upadacitinib groups met the primary endpoints of EASI-75 and IGA 0/1 at week 16. In particular, in Measure Up 1, EASI-75 was achieved by 80%. 70% and 16%, while IGA 0/1 by 48%, 62% and 8% of patients receiving upadacitinib 30 mg, 15 mg and placebo, respectively. In Measure Up 2, EASI-75 was achieved by 73%. 60% and 13%, while IGA 0/1 by 39%, 52% and 5% of patients receiving upadacitinib 30 mg, 15 mg and placebo, respectively. In addition, the secondary endpoints were also met in all active-treatment groups in both studies, including the achievement of EASI-90 and reduction of itch [22]. The analysis of the follow-up data of these two RCTs demonstrated sustained efficacy and a favorable safety profile of upadacitinib through week 52 [23].
The efficacy of upadacitinib combined with TCS in patients aged 12 to 75 years old suffering from moderate to severe AD has been assessed in a further double-blinded, placebo-controlled, phase III RCT [AD Up (n = 901)]. Eligible patients were randomized 1:1:1 to receive TCS combined with either upadacitinib 30 mg, or upadacitinib 15 mg or placebo for 52 weeks [24,25]. Both upadacitinib groups met the two primary endpoints of EASI-75 and IGA 0/1 at week 16. EASI-75 was achieved by 77.1%, 64.6% and 26.4% of patients receiving upadacitinib 30 mg, 15 mg and placebo, respectively, (p < 0.0001) while IGA 0/1 by 58.6%, 39.6% and 10.9% of those receiving upadacitinib 30 mg, 15 mg and placebo, respectively, (p < 0.001) [24,25]. Moreover, both upadacitinib arms demonstrated a rapid onset of action as soon as week 2 regarding skin clearance and week 1 regarding pruritus, which could be maintained through week 52 [24].
Upadacitinib has been also tested in a head-to-head phase III RCT, that compared upadacitinib 30 mg once daily with dupilumab 300 mg every 2 weeks for 24 weeks in adults with moderate-to-severe AD [Heads Up (n = 692)] [26]. More patients from the upadacitinib group achieved the primary endpoint of EASI-75 at week 16 compared to dupilumab [71% and 61.1%, respectively (p = 0.006)]. In addition, the secondary endpoints (EASI-75 at week 2, EASI-100 at week 16, improvement of pruritus) demonstrated the superiority of upadacitinib vs. dupilumab [26].
Regarding safety, the most commonly reported treatment-related AEs in the phase III trials were acne, nasopharyngitis, headache, upper respiratory tract infections and CPK elevation. There were isolated herpetic infections in both Measure Up 1 and Measure 2 trials (Measure Up 1: herpes zoster in six patients of the 30 mg group and five of the 15 mg group; eczema herpeticum in three cases of the 30 mg arm and four of the placebo arm; Measure Up 2: herpes zoster in three patients of the 30 mg arm, in six patients of the 15 mg arm and two patients of placebo arm; eczema herpeticum in three cases of the 15 mg arm). Apart from one venous thromboembolism in Measure Up 2, there were no opportunistic infections, active tuberculosis, or deaths [22]. The AE profile of upadacitinib was similar also in the phase III, AD up trial, in which there was no difference in AEs incidence between the treatment and placebo arms [24,25]. In the Heads Up study, the rates of herpes zoster, eczema herpeticum, other serious infections and laboratory abnormalities were higher in the upadacitinib arm and consistent with the aforementioned phase III studies, whereas rates of conjunctivitis and injection-site reactions predominated in the dupilumab arm [26]. Finally, a double-blinded, placebo-controlled phase III RCT, designed to evaluate the safety of upadacitinib combined with TCS in Japanese adult patients with moderate to severe AD [Rising up (n = 272)], demonstrated consistent results with those of prior multicentric phase II and III studies, as reported above [27].

Abrocitinib
Abrocitinib is a selective JAK1 inhibitor, that has already been evaluated in phase III RCTs in adolescents and adults. A dose-ranging, double-blinded, placebo-controlled phase IIb RCT (n = 267) evaluated the efficacy of abrocitinib 200 mg vs. 100 mg vs. 30 mg vs. 10 mg vs. placebo once daily for 12 weeks in adults with moderate to severe AD [28]. The primary endpoint of IGA0/1 was achieved by 43.8% vs. 29.6% vs. 5.8% of patients from the abrocitinib 200 mg vs. 100 mg vs. placebo, respectively. Reduction in the EASI score (secondary endpoint) was also dose-dependent, with greater reductions being observed in patients treated with abrocitinib 200 mg or 100 mg [28].
A following placebo-controlled phase III RCT evaluated the efficacy of abrocitinib compared to dupilumab in adult patients with moderate to severe AD [JADE COMPARE (n = 838)] [33]. Study participants received either abrocitinib 200 mg daily or 100 mg daily or dupilumab 300 mg every other week or a placebo for 16 weeks, while the concomitant use of topical therapy was not prohibited, unlike JADE MONO studies. The achievement rate of primary endpoints (EASI-75 and IGA0/1) was found to be statistically significant only between abrocitinib vs placebo [EASI-75: 70.3% (p < 0.001) vs. 58.7% (p < 0.001) vs. 58.1% (p > 0.05) vs. 27.1%; IGA 0/1: 48.4% (p < 0.001) vs. 36.6% (p < 0.001) vs. 36.5% vs. 14%; for abrocitinib 200 mg vs. 100 mg vs. dupilumab vs. placebo, respectively]. Statistical comparison between abrocitinib and dupilumab, regarding IGA0/1 and EASI-75 was not among the trial endpoints and was not assessed. Abrocitinib 200 mg was also found to be superior to dupilumab with respect to itch response at week 2 but both treatments were equal as regards most other secondary endpoints at week 16 [33]. Abrocitinib has been also found to be effective in patients with moderate to severe AD, who previously received dupilumab, according to the results of the JADE EXTEND (n = 203) clinical trial [34].
Safety analysis revealed a greater number of treatment-related AEs in abrocitinib groups in both JADE MONO trials (JADE MONO-1: 69-78% vs. 57%; JADE MONO-2: 63-66% vs. 42%; for abrocitinib groups vs. placebo, respectively), while similar percentages of AEs occurred also in the JADE COMPARE trial [29,30,33]. However, the rate of SAEs was similar between all groups including placebo (between 3.2% and 3.9% in JADE MONO-1 and 1.3% and 3.2% in JADE MONO-2). There was only one sudden cardiac death in the 100 mg abrocitinib group of the JADE MONO-2 trial during the follow-up period, which was not attributed to study treatment. In both studies, nasopharyngitis, nausea, headache and upper respiratory tract infections were the most commonly documented AEs. Herpes simplex infection was only observed in the abrocitinib groups in the JADE MONO trials and in the abrocitinib as well as placebo group in the JADE COMPARE trial. Thrombocythemia was reported in both trials, similarly to the previous phase II study, but it was self-limiting despite treatment continuation [29,30,33].

Other Oral JAKs
The efficacy and safety of SHR0302, a highly selective JAK1 inhibitor, in adult atopic dermatitis has been assessed in a double-blinded, placebo-controlled phase II RCT (n = 105) conducted in China between October 2019 and August 2020. Patients were randomized 1:1:1 to receive SHR0302 8 mg vs. 4 mg vs. placebo for 12 weeks. The primary endpoint of IGA 0/1 was met by 54.3% (p < 0.001) vs. 25.7% (p = 0.022) vs. 5.7% of patients in the 8 mg vs. 4 mg vs. placebo group, respectively. Regarding the safety profile, no serious AEs or infections were reported, and most AEs were mild in severity [35].

Comparative Data Regarding the Efficacy and Safety of JAK Inhibitors in AD
Despite the growing utilization of JAK inhibitors in AD based on the availability of efficacy and safety data from phase III RCTs, head-to-head comparisons of different JAK inhibitors are lacking and comparative data may yet only be derived from published meta-analyses.
A network meta-analysis of six RCTs conducted by Zhang et al. in 2021 compared the efficacy of six JAK inhibitors (abrocitinib, baricitinib, upadacitinib, tofacitinib, ruxolitinib, delgocitinib) based on the achievement of EASI-50. Among all included, oral upadacitinib 30 mg twice daily exhibited the best efficacy, while topical delgocitinib 3% twice daily was found to be superior to other topical agents [73]. A similar meta-analysis of 14 RCTs conducted by Li et al. in 2021 evaluated the efficacy of topical and oral JAK inhibitors based on the weighted mean difference (WMD) of EASI and IGA response [74]. The largest improvement of both markers was achieved at week 4 and topical agents were found to be significantly more efficacious than systemic ones. Upadacitinib was the most effective in reducing EASI while abrocitinib was most effective regarding IGA response at week 4 of treatment. Although abrocitinib showed a higher incidence of AEs, the relative risk of AEs leading to treatment discontinuation did not differ among all JAK inhibitors [74]. Finally, Silverberg et al. conducted a network meta-analysis of 11 RCTs to assess the efficacy of available targeted systemic therapies (abrocitinib, upadacitinib, baricitinib, dupilumab, tralokinumab) without topical therapy on moderate to severe AD. The study included analysis of IGA0/1, EASI-75, EASI-90 and assessment of pruritus at weeks 12 and 16 and at earlier time points. In all, upadacitinib 30 mg daily was found to be the most efficacious systemic treatment (based on all endpoints), followed by upadacitinib 15 mg daily and abrocitinib 200 mg daily [75].
In September 2021, the US Food and Drug Administration (FDA) reported that there is an increased risk of serious heart-related events, including stroke, myocardial infarction, blood clots and death with oral tofacitinib, based on the results of the ORAL surveillance study of tofacitinib versus TNFa in patients with rheumatoid arthritis. Recently, the "black box warning" was added to all currently approved oral JAK inhibitors, including baricitinib and upadacitinib, based on the fact that these agents share common pharmacological mechanisms with tofacitinib [76]. However, larger safety clinical trials must elucidate if this extrapolation of data from the pan-JAK inhibitor tofacitinib to the selective JAK inhibitors is justified, or if the differences in their specificity for the four JAKs as well as in the pathophysiology of the treated disease (Table 3) is related to different safety profiles among these agents [72,77,78].

Conclusions
Baricitinib, Upadacitinib and Abrocitinib have been approved for the treatment of moderate/severe atopic dermatitis in patient candidates for systemic therapy by both the Food and Drug Administration and the European Medicines Agency. Selective targeting of the JAK/STAT pathway allows for the control of multiple cytokines involved in AD pathogenesis, such as IL4, 13, 31, 22 and TSLP rather than blocking a single inflammatory interleukin, such as with monoclonal antibodies. Furthermore, the broad-spectrum immunosuppression caused by cyclosporin, azathioprine and mycophenolate is avoided. Literature suggests that this mode of action along with the pharmacokinetics of JAK inhibitors might account for excellent efficacy, rapid onset and long-term control of eczematous lesions and itch in AD combined with an acceptable safety profile, despite the need for screening and monitoring. Data are becoming available for both the JAK/STAT pathway and the IL blocking, along with the investigation of the other pathways involved in AD pathogenesis, such as TSLP and OX40 blockade, resulting in novel insights into the immunological components of the disorder, as well as a new era of treatment options for the patients.

Conflicts of Interest:
The authors declare no conflict of interest.  Th  T helper  TNF-a  tumor necrosis factor alpha  Treg  T regulatory  TSLP  thymic stromal lymphopoietin  TYK2  tyrosine kinase 2


Geometric Morphometrics Reveal Shape Differences in the Toes of Urban Lizards

Synopsis Urbanization, despite its destructive effects on natural habitats, offers species an opportunity to colonize novel niches. Previous research found that urban Anolis lizards in Puerto Rico had increased adhesive toepad area and more ventral toepad scales, traits that are likely adaptive and genetically based. We further investigated these phenotypic changes using geometric morphometrics to measure differences in toe shape, toepad shape, and lamellar morphology. Our results indicate that the increased toepad area of urban Anolis cristatellus lizards in Puerto Rico is not simply an isometric increase in toe size. Toes of urban populations exhibit multiple disproportional changes compared to forest lizards, with a larger proportion of the toe length covered in adhesive toepad. In addition, the toepads of urban lizards increase more in length than width. Lastly, lizards in urban populations exhibit both increased number of lamellae as well as increased spacing between individual lamellae. We also observed regional variation, with urban specimens having significantly more disparity, suggesting similar processes of urban adaptation are likely happening in parallel across the island, yet with region-specific idiosyncrasies, possibly generating more variation in toepad morphology across urban specimens as compared to forest specimens. Considering the use of geometric morphometrics, we found that specimen preparation, specifically how flat and straight toes are during imaging, to be an important factor affecting our data, more so than specimen size or any other meaningful morphological variation. In addition, we found that landmark and semilandmark data can be used to directly estimate toepad area, offering the opportunity to streamline future studies. In conclusion, our results highlight the value of considering toepad morphology in more detail beyond adhesive pad area or number of lamellae. Geometric morphometrics tools may be employed to elucidate subtle differences in shape to better allow researchers to connect changes in morphology to ecology and adhesive performance.

Introduction
As ecosystems change over time, new niches become available for organisms to exploit. This ecological opportunity can lead to speciation and novel ecological communities ( Yoder et al. 2010 ). The process often occurs when organisms colonize isolated uninhabited environments such as islands, lakes, and mountain tops. The evolution of Caribbean Anolis lizards provide an incredibly well-studied example. Within each of the Greater Antillean islands, anoles diversified along the same axes of structural and climatic habitat, with species using the same microhabitat evolving convergent morphology including limb length, tail length, body size, head shape, color, and adhesive toepad size (reviewed in Losos 2009 ). Consequently, trait-environment relationships in anoles are widely considered an evolutionary model for studying ecological speciation, functional morphology, convergence, and adaptation.
There is also suggestion that novel functional traits (i.e., key innovations), rather than island colonization, can drive expansion into previously unoccupied niches leading to increased diversification ( Burress & Muñoz 2022 ). One hypothesized key innovation is adhesive toepads in lizards ( Garcia-Porta & Ord 2013 ;Miller & Stroud 2022 ). Anole toepads are similar to those of geckos and some arthropods ( Irschick et al. 1996 ;Labonte & Federle 2016 ). Anole toepads are composed of lamellae, or scansors; specialized scales on the ventral surface covered in a dense array of setae, tiny hair-like projections that generate van der Waals forces enabling these lizards to cling to surfaces ( Ruibal & Ernst 1965 ;Autumn et al. 2002 ). Previous research has shown that toepad area is positively correlated with cling force, presumably because larger toepads have more setae to engage surfaces at the microscopic level ( Irschick et al. 1996 ;Irschick, Herrel & Vanhooydonck 2006 ). In addition, species that perch higher in arboreal habitats possess larger toepads with more lamellae, presumably conferring greater clinging ability for life in the treetops ( Glossip & Losos 1997 ;Macrini, Irschick & Losos 2003 ). Although toepads have likely played an important role in anole specialization, anoles use their toepads in conjunction with their claws to cling to rougher perches ( Macrini, Irschick & Losos 2003 ;Crandell et al. 2014 ;Naylor & Higham 2019 ;Yuan, Wake & Wang 2019 ;Falvey et al. 2020 ). It has been suggested that perch texture may affect pad performance and hence influence the evolution of toepad area ( Winchell et al. 2018b ). If a population uses smooth substrates, where claws may be less useful, it may be advantageous to have larger toepads to produce greater traction and reduce falls ( Naylor & Higham 2019 ). Given our strong background knowledge of anole ecomorphology, anoles make a great study system to investigate how a key innovation such as adhesive toepads can be shaped by adaptation to novel environments and drive ecological speciation ( Garner et al. 2019 ).
In the Anthropocene, ubiquitous human impacts have drastically changed the face of the earth. One of the clearest examples of this drastic modification is deforestation and the replacement of wilderness areas with urban environments ( Forman 2014 ). While this has led to the loss and displacement of biodiversity worldwide ( McKinney 2002 ), recent studies have found that urban environments can be a source of ecological opportunity leading to evolutionary change ( Johnson & Munshi-South 2017 ). The crested anole ( Anolis cristatellus ; Duméril and Bibron, 1837 ) is an urbanophillic species of anole lizard, occupying forest and urban habitats throughout its native range in Puerto Rico ( Winchell et al. 2016 ;Winchell et al. 2018a ). Prior research on this species has documented behavioral, ecological, morphological, and physiological shifts between adjacent urban and forest populations within its native range ( Winchell et al. 2016 ;Winchell et al. 2018aWinchell et al. , 2018bAvilés-Rodríguez & Kolbe 2019 ;Campbell-Staton et al. 2020 ).
Although responses to and tolerance of urbanization in anoles likely involve a combination of traits , two axes of adaptive morphological differentiation stand out as important in urban environments. Urban populations of A. cristatellus exhibit heritable morphological shifts towards longer limbs and augmented toepads, that is, larger with more lamellae ( Winchell et al. 2016 ;Winchell et al. 2018b ). Previous studies suggest that perch characteristics are important factors in how pad bearing lizards use urban environments ( Taylor, Daniels & Johnston 2015 ;Winchell et al. 2018a ). Increases in toepad area and lamella number in urban A. cristatellus are presumably adaptations to smoother surfaces more often encountered in urban environments. Winchell et al. (2018b) found that lizards with larger toepads and more lamellae outperformed those with smaller toepads and fewer lamellae on smooth surfaces where claws were likely ineffective. Toepads likely enhance locomotor performance on these surfaces by reducing loss of traction and aid in counteracting the backward pitching effect generated by longer limbs when climbing vertical surfaces ( Kolbe et al. 2015 ;Winchell et al. 2018b ;Battles, Irschick & Kolbe 2019 ;Naylor & Higham 2019 ). These phenotypic shifts have likely conferred a performance advantage to urban A. cristatellus , particularly on smooth vertical substrates such as painted walls and fences common to urban environments ( Winchell et al. 2018b ). In addition, these patterns of urban adaptation are likely not unique to A. cristatellus . Parallel changes have been observed in Anolis carolinensis ( Cuvier & Voigt 1832 ). Urban A. carolinensis in park-like habitats in Louisiana have larger toepads than those found in nearby forested habitat, an effect the authors hypothesize is due to the use of smoother broadleaf vegetation in anthropogenic environments ( Irschick et al. 2005 ).
We further investigated how toepads have changed in urban populations of A. cristatellus . Using geometric morphometrics, we move beyond simply measuring toepad area and lamella number, or even taking linear measurements of toepad shape, to more precisely and comprehensively quantify how toepad shape has changed in urban anoles. Specifically, we explore how increased toepad area of urban anoles is achieved, whether isometrically, with proportional changes in both pad length and width, or through disproportional changes in toepad width or in length ( Fig. 1 ) resulting in an overall shape change in addition to an increase in size. We also considered how urban anoles fit more lamellae into their large toepads ( Fig. 1 data modified from Winchell et al. 2018b ), asking how pad size and lamellae organization within the pad has changed in urban populations. As this study represents one of the first to apply geometric morphometrics to toepad morphology, we also conducted a set of analyses to describe our approach. This deeper investigation into toepad mor-  Winchell et al. (2018b) found urban A. cristatellus have larger toepads with more lamellae. Using data from Winchell et al. (2018b) , we highlight these differences (urban individuals in grey, forest lizards in green) for rear toepad area and lamella counts plotted against SVL (A). Plotting average residual toepad area and lamella counts (with 95% confidence intervals) highlights significant differences between urban and forest lizards (B). Hypothetically, toepad area may change isometrically with no change in overall pad shape (C top illustration) or nonisometrically through a disproportionate increase in toepad width or length (C center and bottom illustrations). Our analyses investigate which of these scenarios best describes previously observed toepad differences between urban and forest populations.
phology, an adaptive trait in urban and non-urban populations of anoles and geckos, will foster future functional morphology and biomechanical studies investigating relationships between toepad size, shape, ecologically relevant performance, and microhabitat use.

Materials and methods
We captured live animals from the wild and imaged their toepads for subsequent analyses. We used highdimensional analyses of variance, canonical variance analyses for visualization, descriptive principal component (PC) analyses, and linear morphometric analyses to investigate how toepad shape and lamellar morphology differ in urban anoles. We also conducted analyses generally exploring the use of geometric morphometrics to study toepad shape.

Field methods
We sampled adult male A. cristatellus from 13 sites across Puerto Rico (7 urban and 6 forest sites) in 5 municipalities (San Juan, Arecibo, Aguadilla, Mayagüez, Ponce) between 2012-2016 ( Winchell et al. 2016 ;Winchell et al. 2018b ). See table S1 for locality sampling information. Paired urban and forest sites within each municipality were within 10 km of each other. Urban sites were residential neighborhoods and university campuses, characterized by extensive impervious surfaces, minimal tree cover, and abundant anthropogenic structures (e.g., buildings). Forest sites were mature secondary forests, including tropical dry forest (Ponce), and moist forest (San Juan, Arecibo, Aguadilla, Mayagüez). Forest sites had extensive canopy cover and few anthropogenic structures except for walking paths, although some forest sites (e.g., San Juan) were immediately adjacent to urbanization (and so we cannot rule out that lizards do not encounter this habitat at some point during their lifespan).
We captured lizards as encountered using standard methods (extendable pole with floss lasso or hand capture) and did not discriminate between lizards on anthropogenic substrates versus natural substrates. We transported lizards to a field laboratory in each Fig. 2 We used geometric morphometrics to investigate changes in toe, toepad, and lamellae shape. Image (A) shows the location of our first nine landmarks and eight curves capturing the shape of the toe and pad. Image (B) illustrates the location of landmarks 10 through 19 and the five curves connecting them to illustrate the size and location of the 5th through 10th lamellae. We also extracted linear measurements of toepad shape. Image (C) illustrates line segments to quantify toepad width, toepad length, and the length of the proximal toe segment. Image (D) illustrates how we quantified average lamellae height based on the distance between consecutive lamellae landmarks. Lastly, we outlined toepad area to the nearest phalangeal joint below the pad to estimate toepad area and counted lamellae within this area (E). municipality where we imaged rear feet using a flatbed scanner (Epson V300) at 2100-2400 dpi (see Winchell et al. 2016 ). A single researcher (KMW) measured snout-vent length (SVL) to the nearest 0.5 mm using a transparent ruler. Animals were released after imaging at their location of capture.

Data collection
The above sampling represents a large dataset spanning 4 years and including over 1000 lizards ( Winchell et al. 2016 ;Winchell et al. 2018b ) from which we selected 13-20 individuals at random (using a random number generator) from each population for our analyses. We included all available individuals for sites with < 20 imaged individuals (Table S1). Focusing on the longest rear toes (fourth digit), we decided between each in-dividual's left or right rear toe by choosing the toe of higher quality imaging (i.e., straight and uniformly flattened with no damage). Images of right rear toes were mirrored so that all images were comparable. We focused on the longest rear toe since it has been shown to be important for ecologically relevant performance capabilities in other lizard groups ( Zani 2000 ;Schulte et al. 2004 ). We developed a novel geometric morphometrics landmarking scheme to capture the entire toe shape, toepad shape, and size and location of the 5th through 10th lamellae ( Fig. 2 ). It is worth noting that geometric morphometric analyses do not explicitly require developmentally homologous landmarks, so long as landmarks are consistently reproducible across specimens. In addition, recent research suggests the distal lamellae of anoles may be developmentally homologous and distinct from more proximal lamellae ( Griffing et al. 2022 ). All images were landmarked by a single researcher (BKH) using tpsDig2 ( Rohlf 2010 ). We used 19 landmarks and 13 curves for each specimen. Curves initially contained 10 semi-landmarks with the first and last semi-landmark of each curve overlapping anchoring landmarks, which were removed in the initial steps of our analyses, resulting in each curve containing eight semi-landmarks.
Landmarks one and two were placed at the base of the toe on the left and right sides where the toe joins the palm, with landmark one at the base of toe three and four and landmark two near toe five. Landmarks three and four were placed on the left and right sides of the toe, respectively, where the proximal end of the toepad begins to widen. Landmarks five and six were placed at the widest left and right points of the toepad. Landmarks seven and eight were placed on the left and right sides of the distal end of the toepad, and landmark nine was placed at the tip of the toe, on the ventral, proximal base of the claw. Landmarks were connected with curves to outline the toe (eight curves total). We also placed landmarks to capture the size and relative location of lamellae. We identified the 5th through 10th lamellae counting distally to proximally, placing landmarks on the left and right sides of each (landmarks 10-19) with a curve outlining the free edge of each lamella. We focused on these specific lamellae as they are large and can be easily traced and therefore are more reliable for comparison. It is worth noting that although the 5th-10th lamellae did not always overlap with the widest part of the toe, they often did. We calculated a millimeter/pixel scale ratio for each image from included scale bars if included. If a scale bar was absent, we used the known scan resolution to calculate a millimeter/pixel scale ratio. Since all specimens were imaged on a flatbed scanner (i.e., there was no variation in the distance from the imaging source), the number of pixels per millimeter in an image is an accurate metric to convert pixel distance to linear distance.
To complement our landmark data, we counted the number of lamellae and measured toepad area for each focal toe using ImageJ with the ObjectJ plugin ( Abràmoff, Magalhães & Ram 2004 ). Deciding how far in the proximal direction to count lamellae and define toepad area can be difficult; previous studies have defined the proximal toepad boundary at either phalanx three or four ( Glossip & Losos 1997 ;Hahn & Köhler 2010 ) whereas others have relied on lamellae width . We defined toepad area to include the dilated area of the toe and the area directly below to the nearest phalangeal toe joint, consistent with Winchell et al. (2016Winchell et al. ( , 2018b ( Fig. 2 ). We did not include area above the dilated pad area in the distal direction towards the claw. We counted the number of lamellae within this area three times per digit, retaining the consensus number of the three counts (as in Winchell et al. 2016 ;Winchell et al. 2018b ). Toepad area was measured multiple times for a subset of individuals (N = 150) to assess repeatability; we used the average of all area measures per toe in our analyses. Repeatability of toepad area measurements was high (intraclass correlation coefficient: 0.997) and within group variation was far less than among group variation (within: 0.010, among: 2.745). These data were collected by KMW.

Geometr ic mor phometr ics analyses
To process and analyze our data, we used the R studio interface (v. 3.0-2, Fox & Weisberg 2019 ) packages. We performed a generalized Procrustes analysis to align our landmark data. This initial step aligns all specimens and removes absolute size. We optimized our semi-landmark locations by minimizing bending energy, a standard approach to optimizing semi-landmark locations. To review our data for errors, we plotted each individual's landmark and semi-landmark locations both before and after alignment (see Supplemental Fig. S1). Outlying individuals were reviewed, and either re-landmarked or replaced when possible if their toes were bent or distorted. This resulted in a final dataset of 246 individuals.
We used a series of complementary analyses to test for statistically significant shape variation and to describe group variation in toepad shape. First, to investigate variation in toepad shape across individuals from urban and forest environments and test the hypothesis of isometric change, we correlated our aligned toepad shape with our categorical habitat variable (urban or forest) using a type-2 sum of squares analysis of variance including an interaction with municipality using the geomoprh "procD.lm" function. Given the Procrustes alignment process, if toepads of urban lizards are simply isometrically scaled versions of forest toepads, then there should be no significant aspects of shape that correlate with habitat. We also tested for significant differences in shape disparity between urban and forest lizards using the "morpho.dist" function in geomorph , again including municipality in our model. These analyses are designed for high-dimensionality shape data.
In addition to our analyses of variance, we used PC analyses to better understand shape variation in our dataset. To generate our PC data, we used the geomorph function "gm.prcomp." We extracted PCA axes 1, 2, and 3 for each individual and projected toepad shape at the extremes of these three PC axes. We also considered To visualize differences in toepad shape between urban and forest lizards, we plotted (A) the mean shape before Procrustes alignment to visualize absolute size and shape differences between habitat types, and (B) mean shape after Procrustes alignment to visualize only differences in shape. (C) Our CV analyses also produced projected extreme urban and forest shapes representing the maximum and minimum shapes along our CV axis. All displayed shape pairs are vertically aligned along the base of the toepad (dashed black lines). Green horizontal lines demarcate the top of forest toepads while grey horizontal lines demarcate the top of urban pads. Grid increments in A represent 0.5 mm while grid increments represent 0.02 units in Procrustes alignment space (B & C). A histogram of our single CV axis (D) highlights the difference in shape between urban and forest individuals (see Fig S3 for an evaluation of significance).
between-group PC analyses (bgPCA) using the Morpho R package to visualize shape differences between urban and forest lizards and across municipalities (see Supplement Figs. S4-S5).
We visualized shape differences between urban and forest populations and across municipalities (as well as absolute differences in size) with the average pre-alignment and aligned shape of urban and forest lizards ( Fig. 3 , see Fig. S2 for municipality specific comparisons). To visualize the relationships between shape and habitat category suggested by our above analysis, we used Canonical Variate Analysis (CVA), implemented with the "CVA" function in the Morpho R package. Similar to PC analyses, CV analyses also identify axes of variation in highly dimensional data, but instead of identifying the most variable dimensions, a CVA maximizes differences between assigned groups. We note that CV analyses can find axes that partition provided groups by chance. In fact, this is expected given the high dimensionality of geometric morphometrics data. Consequently, we did not consider our resulting one-dimensional CVA as evidence of significant differences in shape between urban and forest individuals (these conclusions were based on our analysis of variance described above). To provide context to our CVA results, we generated a set of CVA results using randomized urban and forest habitat category assignments across individuals. Using 100 randomized datasets, we compared the resulting set of CVA results to our actual CVA results to demonstrate the degree in which a CVA finds axes that separate our data by chance (Fig. S3).

Linear measurements and analyses
The Procrustes alignment process includes isometric enlarging and shrinking of individual specimens to eliminate the effect of size while maintaining shape. Because of these individual-specific transformations, if our results suggest non-isometric differences in toepad shape between urban and forest lizards, it can be difficult to infer if individuals with larger toepads had wider or longer toepads. As a result, we also extracted linear measurements from unaligned landmarks data ( Fig. 2 ) to investigate differences in toe, toepad, and lamellae proportions. We measured toepad width (distance between landmarks 5 and 6) and toepad length (distance between the midpoint between landmarks 3 and 4 and the midpoint between landmarks 7 and 8). We determined the length of the toe segment from the palm to the proximal base of the pad using the distance between the midpoint of landmarks 1 and 2 and the midpoint of landmarks 3 and 4. Lastly, we measured average lamella height (not in the vertical direction, but in the proximaldistal direction) by averaging together the distances between adjacent edge landmarks 10 and 12, 12 and 14, 14 and 16, 16 and 18, 11 and 13, 13 and 15, 15 and 17, and 17 and 19 ( Fig. 2 ).
If differences between forest and urban lizard toepads were simply isometric differences in size, then linear proportions, such as the ratio of pad width to pad length, would be similar between urban and forest populations. We calculated the natural-log-transformed ratio of toepad width to length for each individual and compared these data between urban and forest lizards using an independent two-group two-tailed t -test. This analysis evaluates if urban toepads increased more in width (higher width-to-length ratio) or length (lower width-to-length ratio), or if they are simply isometrically larger (non-significant differences between urban and forest width-to-length ratios). We also investigated how changes in toepad size relate to the rest of the toe using independent two-group two-tailed t -tests. We asked if total toe length (natural-log-transformed sum of pad length and proximal toe segment length) differs between urban and forest lizards, which we expect since urban lizards are generally larger. We also asked if the natural-log-transformed ratio of pad length to total toe length differs between urban and forest lizards and if the natural-log-transformed proximal toe length differs between urban and forest lizards to determine if the pad is covering similar proportions of the toe in urban and forest populations.
To complete our investigation using linear measurements, we next focused on lamellae and their relationship with habitat category and the rest of the toe. While previous studies considered toepad area and lamellae counts relative to lizards' body size, we did not standardize our lamellae data by body size in order to directly compare our lamellae data to other dimensions of toepad size and shape. We are interested in aspects of lamellae morphology regardless of if they are found on large bodied or small bodied individuals. We asked if our data confirm previous studies finding urban lizards have larger toepads and more lamellae via two-tailed t -tests using natural log transformed values of area and square root transformed lamellae counts ( Whitlock & Schluter 2015 ). We then asked if lamellae height (in the proximal-distal direction) differs between urban and forest populations, again using a t -test with natural log transformed data. To identify potential interactions between toepad area, toepad length, lamellae number, and lamellae height, we used the "Anova" function in the car R package to conduct type-2 analyses of variance using generalized linear models. We first asked how area is related to lamellae structure, testing if natural-logtransformed toepad area is correlated with natural-logtransformed lamellae height, square root transformed lamellae count, or the interaction of lamellae height and count. We then considered lamellae height and count, asking if natural-log-transformed lamellae height or square root transformed lamellae count is predicted by natural-log-transformed pad length, habitat category, or their interaction. Lastly, we evaluated the correlation between lamella height and count.

Methodological analyses
Given our novel approach of applying geometric morphometrics to investigate adhesive toepad shape, we also conducted a series of analyses to further elucidate the use of this method. We evaluated the utility of estimating area using geometric morphometric landmarks compared to two different measurements of toepad area and centroid size. We also evaluated the presence of toepad allometry and how size was captured by our PC axes. These methods and results are described in the supplement.

Geometr ic mor phometr ics results
Our first statistical comparisons correlating toepad shape with habitat type (urban or forest) and municipality using a Procrustes ANOVA found significant relationships between shape and habitat ( P = 0.001), municipality ( P = 0.001), as well as a significant interaction between habitat and municipality ( P = 0.001). We also found a significant difference in toepad shape variation (shape disparity) between urban and forest lizards after accounting for municipality: urban toepads were 18% more variable than forest toepads ( P = 0.026). To visualize these results, we plotted the mean shape of urban and forest lizard toepads before and after Procrustes alignment ( Fig. 3 ) as well as mean shapes from individual municipalities (Fig. S2). We observed striking difference in absolute toepad size in urban lizards ( Fig. 3 A), with larger toepads and longer toes. After specimens were scaled during our Procrustes alignment process ( Fig. 3 B), subtle differences in shape can be seen. Urban pads still appear larger, with lamellae 5-10 shifted distally, while the proximal section of the toe, from the base of the pad to the palm of the foot, appears shorter.
We further visualize these differences in our CVA projections ( Fig. 3 C). Since we only had two habitat categories, our results produced a single CV axis ( Fig. 3 ) that captured 8.6% of the variation and in which urban and forest individuals are clearly separated. This analysis separated urban and forest individuals more than expected by our null dataset, generated by 100 simulated CVAs (Fig. S3). Specifically, our randomized CV analyses found axes that separated randomly assigned categories by an average distance of 4.28 canonical variate score units (standard deviation: 0.47), with the largest observed separation among all 100 simulations of 5.81 units. Our actual urban and forest groups were separated by 5.89 units, suggesting the differences in shape between urban and forest lizards characterized by the CVA are meaningful. We projected toepad shape along our single identified CV axis, generating shapes that represent a projected urban toepad shape and a forest toepad shape. Similar to our estimates of urban and forest lizard mean toepad shapes, our projected CV axis shapes also suggest urban lizards have disproportionately longer toepads that make up a larger proportion of the total toe length with lamellae shifted distally ( Fig. 3 ).
We next investigated general patterns in our shape data using PCs analyses ( Fig. 4 ). To provide context to our PC axes, we produced shape estimates of the maximum and minimum axes extremes ( Fig. 4 D). PC 1 explained 23.8% of the variation in our dataset, PC 2 explained 20.5%, and PC 3 explained 11.5% ( Fig. 4 C). PC 1 appears to primarily capture aspects of how straight each toe was at the time of imaging, and specifically whether toes bent to the left (lower PC 1 values) or right (higher PC 1 values). In addition, PC 1 also captures aspects of relative toe and toepad propor-tions. Individuals lower on PC1 have longer proximal toe segments (distance from the base of toe to the base of the pad) and toepads shorter in length, whereas individuals with higher values on PC 1 have proportionally longer pads and/or shorter proximal toe segments. Our projections of PC 2 capture aspects of lamellae location and lamellae height. Lamellae 5-10 of individuals lower in PC 2 are more distally located and are closer together (narrower) suggesting these individuals may have more total lamellae. Conversely, lamellae 5-10 of individuals on the upper end of PC 2 sit lower on the toe and are more widely spaced, suggesting these individuals may have fewer lamellae. PC 2 also appears to capture aspects of toepad length and distal toe segment length (the distance of the claw from the pad), with specimens lower on PC 2 having a longer toepad and shorter distal segment, whereas specimens higher on PC 2 have shorter pads and longer distal segments. Lastly, PC 3 captures toepad size, including both length and width. Individuals higher on PC 3 have proportionally wider toepads that make up a larger proportion of the toe's total length. Our between-group PC analysis using Linear and area measurements help us better understand shape differences. Histogram (A) visualizes the difference in the mean ratio of width/length between urban and forest lizards. Scatterplot (B) illustrates the same data, comparing pad width versus length in millimeters. The slopes of the fitted lines through the origin represent the mean ratios for urban and forest lizards across our dataset. The grey urban line's location below the green forest line illustrates the lower width/length ratio in urban lizards, suggesting urban lizards gain toepad area via larger increases in length than width. We also include histograms of toepad width and length along the X and Y axes to better illustrate differences between urban and forest lizards and the larger increase in pad length compared to width. Mean values for each municipality are plotted with point shape indicating municipality and black lines connecting each intra-municipality urban and forest pair. Similar municipality-specific scatter plots can be seen in Fig. S6. We also compared toepad area, lamellae count, and lamellae height between urban and forest lizards (C and D). We plotted these relationships as scatter plots with embedded histograms for each axis. We included best fit lines not forced through the origin in our scatter plots for both forest and urban populations. In all plots urban individuals are in grey and forest are in green. habitat category (urban or forest) reinforced these results, suggesting urban lizards have proportionally larger pads and shorter proximal toe segments, with lamellae shifted distally (Fig. S4). Our results from the bgPCA investigating variation across municipality suggested strong differences between municipality with similar patterns as seen in our PCA ( Fig. 4 ), with our bgPC axis 1 capturing aspects of lamellae as seen in PC 2, and bgPC axis 2 capturing subtle aspects of toepad size (Fig. S5).

Linear measurements
We further investigated our hypothesis of isometric change with a series of linear measurements extracted from our unaligned landmark data. The toepads of urban lizards had an average width of 1.86 mm (sd = 0.28) and an average length of 3.61 mm (sd = 0.47). For-est lizards' toepads had an average width of 1.69 mm (sd = 0.27) and an average length of 3.20 mm (sd = 0.43). The toepads of urban lizards were on average 10.06% wider than forest lizards, but also 12.7% longer. The mean width to length ratios of toepads differed between urban and forest lizards (two-sided t -test, t = −2.05, df = 231, P = 0.04): urban lizards have a smaller mean width to length ratio (0.514, 95% CI: 0.507-0.522) compared to forest lizards (0.527, 95% CI: 0.518-0.535) ( Fig. 5 ). These results suggest that the larger toepad size of urban lizards is not isometric and that urban toepads increase in length more than width.
We also used our linear measurements to consider how the proportions of the overall toe, and not just the pad, differ between urban and forest lizards. Specifically, we asked how increases in toepad length may also be related to changes in overall toe length and the length of the proximal toe segment ( Fig. 6 and Fig. 6 To visualize differences in lizard toe proportions, we plotted the total toe lengths of urban (grey) and forest (green) lizards versus the length of only the pad (A). We highlight the significantly different ratio of pad length to total toe length between urban and forest lizards using g reen and g rey best-fit lines forced through the origin. Histograms of the pad length/total toe length ratios (B) and the proximal toe segment length (C) include vertical lines representing urban and forest population mean values. Fig. S6). We found that total toe length is longer in urban lizards ( Fig. 6 ; two-sided t -tests, t = 5.7, df = 220.7, P < 0.001) which is not surprising given that urban lizards are often larger and tend to generally have longer limbs (including the metatarsal and first phalanx on toe IV; Winchell et al. 2016 ). Forest lizards have a mean total toe length of 9.8 mm (95% CI: 9.60-10.02 mm), while urban lizards have a mean total length of 10.59 mm (95% CI: 10.42-10.78). We next compared the mean ratio of pad length over total toe length between urban and forest lizards ( Fig. 6 ) and found non-isometric differences between urban and forest lizards. The total toe lengths of urban lizards were comprised disproportionally more toepad (t = 5.6, df = 236, P < 0.001), covering 1.41% more of their toe lengthwise compared to forest lizards. Urban lizards had mean pad length to total toe length ratio of 0.338 (95% CI: 0.334-0.341) while forest lizards had a mean pad length to total toe length ratio of 0.324 (95% CI: 0.32-0.327). We observed municipalityspecific differences in this relationship (Fig. S6). Lizards in Aguadilla and San Juan exhibited minimal differences between urban and forest lizards, whereas more pronounced differences were observed in Arecibo, Mayagüez, and Ponce. Lastly, we investigated change in the proximal toe segment, the distance between the base of the pad and the base of the toe ( Fig. 6 ). We found significant differences in the mean proximal toe length (t = 4.2, df = 215.8, P < 0.001), with urban lizards having significantly longer proximal toe segments (urban = 7.01mm 95%CI 6.90-7.12, forest = 6.63 95%CI 6.49-6.76). Together these results suggest that although urban lizards have longer toes and larger toepads in an absolute sense compared to forest lizards, these differences are not isometric. Urban lizards have disproportionately longer toepads, and their toepads cover disproportionately more of their toe.
Through the above analyses, we noted many municipality-specific differences in toepad shape, linear measurements, and lamellar morphology. In Aguadilla, we did not observe substantial differences between urban and forest specimens in toe or pad proportions (Fig. S2, S6), but we did observe distal shifts in the lamellae (Fig. S2 and S5) and an absence of low-lamellae count individuals in urban populations, resulting in an increase in total average lamellae but no change in lamellae height (Fig. S6). Forest individuals from Arecibo had some of the smallest toes across the five municipalities, with urban individuals having much larger toes (Fig. S2) and pads (Fig. S6). The pads of urban individuals were disproportionately larger with shorter proximal and distal toe segments (Fig. S2,  S6). The pads of urban individuals increased in width more than length ( Fig. 5 and S6) contrary to other municipalities and the overall pattern we observed. Urban individuals also had more and wider lamellae (Fig. S6). In our Mayagüez municipality, urban lizards again had absolutely larger toes with weak relative differences between urban and forest individuals ( Fig. 3 ). We did observe a slight increase in urban lizards' pad proportion of the toe (Fig. S2, S5, S6), although the pad itself appears to be increasing in size isometrically ( Fig. 5 , S6) with equal increases in length and width (Fig. S6). Urban lizards also had more lamellae and wider lamellae (Fig. S6). The toes of forest lizards in Ponce were also very small and here we observed the biggest absolute difference in toe size between urban and forest, with urban toes being much larger. We also found strong relative differences in pad-toe proportions and a distal shift in lamellae (Fig. S2, S5, S6). Urban lizards had a clear increase in pad length more than width ( Fig. 5 ). Although forest lizards already had a high lamellae count, urban lizards increased in pad area, lamellae count, and slightly in lamellae height (Fig. S6). Lastly, forest lizards from our San Juan locality had some of the largest pads when compared with other forest lizards, with subtle differences in absolute size in urban individuals. The toes of urban lizards had disproportionately more pad and shorter proximal toe segments ( Fig. 3 ). Considering shape, urban and forest lizards were similar with urban lizards having slightly longer pads and shorter proximal toe segments (Fig. S2, S5, S6) with a subtle shift in pad length more than width ( Fig. 5 ). Urban lizards also had wider lamellae (Fig. S6).
Finally, we investigated relationships between pad area, lamellae height, and lamellae count. We found both lamellae count (type-2 analyses of variance; P < 0.001) and height ( P < 0.001) were positively correlated with toepad area, but the interaction was not significant ( P = 0.65, see Fig. 5 ). Lastly, we tested the relationship between lamellae count and height, finding it to be not significant ( P = 0.22). These results varied across municipalities (Fig. S6): Aguadilla varied little between urban and forest individuals, whereas urban lizards from Arecibo, Mayagüez, and Ponce all increased in toepad area, lamellae count, and lamellae height. Interestingly, San Juan appeared to increase primarily in lamellae height, with minor changes in area and count. Similarly, we considered how toepad length relates to lamellae characteristics (Fig. S7). We found that lamellae height was positively related with toepad length (ANOVA, P < 0.001), weakly correlated with habitat ( P = 0.066), with no significant interaction between habitat category and toepad length ( P = 0.58). Our analyses of lamellae count yielded similar results, with pad length showing a significant positive relationship ( P < 0.001), habitat category again not quite being significant ( P = 0.09) and no significant interaction ( P = 0.5).
Together these results suggest urban lizards are generally increasing all aspects of their toepads and lamellar morphology by increasing toepad area, lamellae count, and lamellae height, with longer pads also supporting more space between lamellae. Interestingly, these results do not support the idea of a trade-off between lamellae height and lamellae number, but instead suggest increases in toe pad area, specifically toepad length, likely allow increases in both lamellae count and lamellae height.

Discussion
Previous studies have found that urban populations of A. cristatellus in Puerto Rico have larger adhesive toepads and more lamellae when compared to their forest counterparts ( Winchell et al. 2016( Winchell et al. , 2018b. We asked how these populations are developing larger pads, either via isometric changes to the toe or through disproportional changes to their toe and toepad morphology. With the use of geometric morphometrics having never been applied to adhesive toepads, we also explored some of the properties of this approach. Our results reinforce previous findings that urban populations of A. cristatellus have significantly different toes and toepads than those of forest lizards ( Winchell et al. 2016 ;Winchell et al. 2018b ). We found that, in an absolute sense, the toes of urban lizards are longer with larger pads ( Figs. 3 ,5 ,6 ) and that these differences are not isometric. When only considering differences in shape, we found that urban lizards' toes are composed of relatively more pad, (i.e., longer pads and shorter proximal toe segments; Figs. 3 , S4, 6). In addition, the pads of urban lizards are disproportionately longer, relative to their increase in width ( Fig. 5 ). These results suggest adaptation to urban microhabitats has changed not only the relative size and proportions of urban lizards' toepads but aspects of the entire toe. We discuss the potential implications of these findings with respect to urban adaptation as well as more methodological implications in the following sections.

Urban adaptation
Urban lizards are exposed to different selection pressures than those in forested areas, and the observed increase in overall toe length we found in urban lizards could be driven by functional demands. The fourth toe of the hindlimbs in lizards is related to locomotion, producing propulsion forces during sprinting ( Irschick & Jayne 1999 ;Russel & Bels 2001 ). Previous studies have shown that urban areas have greater habitat openness ( Winchell et al. 2016 ;Winchell et al. 2018a ) imposing functional demands to increase sprint speed which could lead to increases in total length of the fourth toe. Similarly, changes in overall toe length, even if not directly influencing sp rint speeds, could be related to increases in hindlimb length, selected for via the use of wider perches, through indirect selection if these traits are functionally or developmentally integrated as some evidence seems to suggest ( Spezzano & Jayne 2004 ;Foster & Higham 2012 ;Poe et al. 2021 ). Related but distinct is the possibility that as lizards expand into urban areas they are presented with ecological opportunity. This explanation would be consistent with the increased disparity in the toe shape of urban lizards we observed, however, this greater disparity could also be attributed to greater heterogeneity in urban selection pressures that are imposed in different municipalities ( Yoder et al. 2010 ). Although, we were not able to explicitly address these potential alternatives, they are testable hypotheses that should be considered in future studies.
In addition to changes in absolute toe length, we also found that the toe's length is comprised of more padded area and less proximal toe length in urban lizards. The increases in pad area (via increases in pad length) are not surprising, as this presumably would allow the lizards to climb and cling more easily to the smoother perches common in urban areas as the claw becomes less effective ( Winchell et al. 2016, Falvey et al. 2020. Potentially, increasing pad area through increasing length (rather than width) could optimize both functional demands that might be placed on A. cristatellus in locomotion both on horizontal (sprinting) and vertical substrates (climbing and clinging) ( Simon et al. 2022 ). However, increases in pad length could also be the result of indirect selection on increases in hindlimb length or total toe length. Decreases in the relative proportions of proximal toe lengths might represent a reduction in phalanges that can be seen in other pad-bearing species of lizards ( Zhuang et al. 2019 ). Russell 1977 suggested that the reduction of more proximal phalanges placed the more distal phalanges in a position that better optimized usage of the claw, which has been shown to work synergistically with toepads to produce cling force in pad-bearing lizards ( Crandell et al. 2014 ;Naylor & Higham 2019 ;Yuan, Wake & Wang 2019 ;Falvey et al. 2020 ). Future studies should consider how toe anatomy changes between urban and forest lizards and if this is related to function or performance, particularly how the claw, pad, and proximal phalanges interact and how this impacts cling force as lizards adapt to using smoother substrates.
While previous studies found that urban lizards have more lamellae on their toepads ( Winchell et al. 2016 ), our results highlight that the lamellae of urban lizards tend to be more distantly interspaced (defined as lamellae height), with lamellae 5-10 more distally located on the toe. This spacing and shift in relative location of lamellae is likely how urban lizards fit more lamellae into the toepad ( Figs. 3 ,7,10). One might expect a negative relationship between lamellae height and lamellae number with the assumption that more lamellae necessitate more tightly packed lamellae, but we found no such trade off. Urban individuals had increases in both lamellae height and count along with pad area. Few studies have considered toepad morphology beyond area and number of lamellae although Macrini, Irschick, and Losos (2003) and Yuan, Wake, and Wang (2019) both found interesting patterns involving toe pad dimensions (especially width), lamellae number, claw morphology, and perch height or diameter, as well as differences between Caribbean and mainland anoles species. These studies suggest pad dimensions and intra-pad lamellae configuration likely influences performance, although the details of this relationship, and how these findings translate at the intraspecific level, remain unclear.
Our results also repeatedly found differences across municipalities, suggesting urban adaptation in toepad shape is not occurring identically across sites in Puerto Rico. This idiosyncratic pattern of parallel adaptation may explain the significant increase in morphological disparity we observed in urban specimens. Although our results confirm a general increase in toepad size and number of lamellae in urban populations, these morphological shifts appear to be occurring in slightly different ways in different populations, similar to previous studies ( Winchell et al. 2016 ). The idiosyncratic variation we saw across municipalities may be examples of plasticity, local adaptation, founder effects, or possibly our modest sample sizes. Genetic analyses suggest that across the island of Puerto Rico there is significant regional differentiation among A. cristatellus ( Reynolds et al. 2017 ), including substantial differentiation associated with environmental variation ( Wogan et al. 2020 ). Thus, the standing genetic variation in populations likely varies across municipalities, potentially influencing adaptive outcomes. Moreover, among the populations included in this study, gene flow is highly restricted across municipalities but not between urban and forest sites within each municipality ( Winchell et al. 2016, Campbell-Staton et al. 2021. In urban-forest pairs where gene flow is elevated, adaptive responses may be constrained. Notably, the urban and forest populations in Ponce were the most genetically differentiated as well as the furthest apart geographically ( Winchell et al. 2016 ) compared to the other municipalities, possibly explaining why we observed such pronounced morphological differences between urban and forest lizards in this and in previous ( Winchell et al. 2016 ) studies. It's also possible that with larger sample sizes these differences would become more pronounced and consistent across municipalities.

Methodological results
Our analyses exploring the use of geometric morphometrics to study toepad shape produced a set of important findings and highlighted some considerations for future studies (see supplement). Our PCs analyses ( Fig. 4 , S9) and the major axes of shape variation in our dataset highlight how specimen preparation can influence a dataset. We recommend future studies planning to measure toe or toepad shape using geometric morphometrics assure that their toes are straight during imaging, otherwise this variation in shape is likely to dominate their dataset. It was not until after specimen condition had been accounted for in PC1 that we saw interesting morphological differences regarding toe proportions, pad proportions, and lamellae location and height differences in our dataset.
We also found specimen size to have an influence on toepad shape, even after size had been removed by Procrustes alignment (Fig. S8 E). When specimen size is defined by toepad area, we see larger specimens having enlarged pad areas, shorter proximal toe sections, and distally shifted lamellae. When specimen size is defined by centroid size, we also see correlations with shape, with centroid size being inflated in specimens with longer toepads (Fig. S8 B). This means that centroid size is likely a poor proxy for toepad area, since short wide toepads will have smaller centroid sizes than long skinny toepads, even with similar areas. Despite this evidence of specimen size and shape being connected in our dataset, our PC analyses did not capture this relationship as expected (Fig. S9). Regardless of how size was defined (as SVL, pad area, or centroid size), it appears PC 2 and PC 3 were capturing aspects of size ( Fig. 4 ), not PC1. This again reinforces the importance of sample preparation during imaging.
Our results also highlight the importance of how toepad size and lamellae count are measured, specifically how proximally on the toe one measures area and lamellae. We found that both centroid size and landmark-based estimates of area can serve as proxies for pad area, although centroid size is affected by specimen compactness (Fig. S8), so it may not always be a great representation of size depending on the specific question being asked and the location of the landmarks for a given study. Depending on the landmark scheme a study is using, the area enclosed by a chosen set of landmarks and semi-landmarks can serve as a proxy for area, or serve as a direct measurement of area, eliminating the need to independently measure area via programs like ImageJ. Lastly, our results highlight how landmarks outlining only a few select lamellae could be informative regarding lamellae count, with distal shifts of our five focal lamellae often indicating more lamellae being present on an individual's toepad.

Conclusions
Our results highlight how microhabitat differences between urban and forest environments are likely driving the evolution of lizard toepad morphology in subtle ways, changing the toe, pad, and lamellae proportions of urban lizards. Specifically urban populations appear to be moving towards morphologies that exaggerate toepad and lamellae related traits, generally fitting more lamellae that are more distally spaced into larger, and disproportionally longer pads, with a larger proportion of the entire toe composed of pad. This process is likely happening repeatedly and independently at urban locations across Puerto Rico, with subtle differences possibly attributable to variation in the amount of time urban and forest populations have been diverging, the extent of population connectivity, or differences in the direction and/or magnitude of selection urban environments impose within each municipality, resulting in urban specimens having more disparate toepads as compared to forest specimens. Our results suggesting urban lizards having longer toes invites questions regarding how an increase in toe length may be related to phalanx kinematics during gripping and previously observed increases in overall limb length in urban populations and how these differences interact with sprint speed and locomotion kinematics. Lastly, while we did not explicitly consider differences in the distal toe segment, the length from the distal end of the toepad to the claw, we did observe differences between municipalities. Future studies may investigate how distal toe segment length above pad influences the use of the claw and how it may be related to the use of smoother urban surfaces ( Falvey et al. 2020 ). Our results also demonstrate that geometric morphometrics is a viable tool for studying adhesive toepad morphology and open a variety of exciting new research directions. In conclusion, our study furthers our knowledge of adaptive toepad morphology, moving beyond simply measuring pad area and lamellae number to gain a more detailed understanding of how adhesive toepads vary and how populations adapt to changing microhabitats with exciting implications for the biomechanics and performance of toepads and the entire locomotor system of pad bearing lizards.