Functional and proteomic analysis of a full thickness filaggrin-deficient skin organoid model

Background: Atopic eczema is an itchy inflammatory disorder characterised by skin barrier dysfunction. Loss-of-function mutations in the gene encoding filaggrin ( FLG) are a major risk factor, but the mechanisms by which filaggrin haploinsufficiency leads to atopic inflammation remain incompletely understood. Skin as an organ that can be modelled using primary cells in vitro provides the opportunity for selected genetic effects to be investigated in detail. Methods: Primary human keratinocytes and donor-matched primary fibroblasts from healthy individuals were used to create skin organoid models with and without siRNA-mediated knockdown of FLG. Biological replicate sets of organoids were assessed using histological, functional and biochemical measurements. Results: FLG knockdown leads to subtle changes in histology and ultrastructure including a reduction in thickness of the stratum corneum and smaller, less numerous keratohyalin granules. Immature organoids showed some limited evidence of barrier impairment with FLG knockdown, but the mature organoids showed no difference in transepidermal water loss, water content or dye penetration. There was no difference in epidermal ceramide content. Mass spectrometry proteomic analysis detected >8000 proteins per sample. Gene ontology and pathway analyses identified an increase in transcriptional and translational activity but a reduction in proteins contributing to terminal differentiation, including caspase 14, dermokine, AKT1 and TGF-beta-1. Aspects of innate and adaptive immunity were represented in both the up-regulated and down-regulated protein groups, as was the term ‘axon guidance’. Conclusions: This work provides further evidence for keratinocyte-specific mechanisms contributing to immune and neurological, as well as structural, aspects of skin barrier dysfunction. Individuals with filaggrin deficiency may derive benefit from future therapies targeting keratinocyte-immune crosstalk and neurogenic pruritus.

report report report report

Introduction
Atopic eczema (also termed atopic dermatitis or eczema 1,2 ) is an itchy inflammatory skin disorder with complex and multifactorial aetiology 3 . The phenotype is heterogeneous and pathogenic mechanisms vary between affected individuals 4 . The pathogenesis may include differential contributions of: immune activation and skin barrier dysfunction; local and systemic effects; multiple genetic and environmental factors 5 . Loss-of-function mutations in FLG, encoding filaggrin, cause the common Mendelian disorder of keratinization, ichthyosis vulgaris (OMIM # 146700). FLG null mutations are also strongly associated with increased risk of atopic eczema 6,7 and multiple other atopic traits 8-10 .
Skin is an organ that can be modelled in vitro to effectively recapitulate the multi-layered structure and gene expression patterns of human skin in vivo 11 . This offers the opportunity for selected molecular mechanisms to be investigated in relative isolation from a complex disease state, using relevant primary cells. Replicate experiments can be performed using primary cells with the same genetic background or cells from different donors 12 . Detailed functional and biochemical analyses may then be applied to define potential therapeutic targets 12,13 .
Filaggrin is a marker of keratinocyte differentiation; it is expressed in the granular layer of the upper epidermis as a polymer -profilaggrin -which undergoes post-translational modification and stepwise proteolysis to release monomeric filaggrin 14 . Profilaggrin and filaggrin have multiple functions, each of which have been described as contributing to the mechanical, biochemical, immunological and microbiological aspects of skin barrier function 15 . Clinical studies have shown that filaggrin haploinsufficiency is associated with xerosis and ichthyosis (dry and scaly skin), keratosis pilaris (increased keratin within skin follicles), palmar hyperlinearity 16 and increased transepidermal water loss 17 . Transcriptomic analysis of full thickness skin biopsies have shown evidence of an abnormal defence response in FLG haploinsufficient atopic skin 18 ; proteomic analysis to assess the effect of FLG knockdown in an epidermal organoid has also shown features of inflammation and stress protease activity 19 . However, in vitro studies have not shown consistent histological or functional effects of filaggrin deficiency in the various different skin organoid models published to date 20 and the multiple mechanisms by which filaggrin deficiency contributes to atopic disease remain incompletely understood 5 .
We have optimised a skin organoid model, with dermal and epidermal compartments cultured using donor-matched primary cells, for functional assessments and global mass spectrometry proteomic analysis. The work aims to investigate in more detail the effect of FLG siRNA-mediated knockdown on one cell type -the keratinocyte -to increase understanding of the filaggrin-deficient phenotype and to further define molecular mechanisms predisposing to atopic skin inflammation.

Source of primary human cells
Primary keratinocytes and primary fibroblasts were isolated from human skin tissue samples obtained, with written informed consent and Ethical Committee approval (East of Scotland Research Ethics Service reference 17/ES/0130 renewal 12/ ES/0083) under governance of the Tayside Biorepository. Surgical surplus samples of clinically normal skin from four adult donors (all females aged 29-65 years; one abdominal and three breast skin reductions) were used for the organoid cultures. Similar samples (n=5) were used for independent biological replicates to test for reproducibility of the functional effects of FLG knockdown.
Organoid culture methods Primary keratinocytes and dermal fibroblasts were isolated from human skin by sequential trypsin EDTA and collagenase D digestion 21 . Using our previously reported methods 12 , the keratinocytes were co-cultured with mitomycin C inactivated 3T3 feeder cells in RM media (3:1 DMEM : Hams F12, 10% FCS, 0.4µg/ml hydrocortisone, 5µg/ml insulin, 10ng/ml EGF, 5µg/ml transferrin, 8.4ng/ml cholera toxin and 13ng/ml liothyronine) (Sigma Aldrich, Gillingham, Dorset, UK) 22 . Epidermal growth factor (EGF) was omitted for the first day of culture. The fibroblasts were cultured in DMEM supplemented with 10% FCS under standard conditions.

Amendments from Version 1
The main aim of this publication is to facilitate sharing of the global mass spectrometry data generated from the skin organoid models with FLG knockdown (described in detail) and functional analyses (described and reported); this main aim has been clarified.
In response to the interest from Reviewers 1 and 3, we have provided an explanation as to why knockdown of FLG at mRNA level is demonstrated to a lesser degree than the knockdown at protein level, reflecting the prolonged duration (10 days) of organoid culture.
We have emphasised that increased dye penetration was observed in some but not all of the immature organoids, as requested by Reviewer 1.
Additional details have been provided to explain how the protein false discovery rate was determined for the removal of false positives, as requested by Reviewer 1.
We have included a fuller explanation of the analysis approach, in which we focus on consistent differences observed across the biological replicate experiments and use thresholds for increased and decreased expression to include roughly equal numbers of proteins.
We have corrected the use of the terms 'ratio' and 'fold-change' where necessary.
A new Figure 10 is provided, following advice from reviewer 2.
The colour-coding has been removed and we have labelled selected data points. During revision of this figure we discovered an error in the plotting of this figure; this detail has been corrected in the revised Methods section and in the re-drawn volcano plot.
The title has been modified to include the term 'full thickness', to clarify the complexity of our model (as recommended by reviewer 2). The list of main findings has now been removed from the title because it had become too long.
Minor typos have been corrected, including two references numbered incorrectly in the previous manuscript.
Any further responses from the reviewers can be found at the end of the article REVISED Dorset, UK) and 0.5ml thrombin (3U/ml in 2 mM CaCl 2 / 1.1% NaCl) (Sigma Aldrich, Gillingham, Dorset, UK) combined on ice, with 200,000 fibroblasts and aprotinin (0.1 U/ml) (Sigma Aldrich, Gillingham, Dorset, UK) then transferred to a 12well plate. After 30 minutes incubation at 37°C, the gels were covered in medium (DMEM, 10% FCS, 0.1 U/ml Aprotinin) and cultured overnight (day 1). On day 2, the medium was replaced with RM excluding EGF, 0.1U/ml aprotinin and 2 × 10 6 suspended keratinocytes. Culture medium was refreshed daily on days 3 and 4 using RM containing 0.1ng/ml EGF and 0.1U/ml aprotinin. On day 5 the gels were carefully removed from wells and lifted onto custom-made steel grids lined with nylon gauze (Millipore, Livingston, Scotland, UK). RM medium supplemented with 0.1ng/ml EGF and 0.1U/ml aprotinin was added up to the base of the dermal equivalent so that the epidermis remained at the air-liquid interface. Medium was refreshed on alternate days sand the cultures were used for analysis up to day 12 (representing three, five and seven days after lifting to the air-liquid interface).
A diagrammatic summary of the process of skin organoid culture is shown in Figure 1.
Epidermis was separated from the fibrin gel using hypertonic saline-induced split (4 hours, 1M NaCl, 4°C) to obtain a keratinocyte-only tissue sample for biochemical analyses.
(ii) Transmission electron microscopy (TEM) Organoid skin samples were fixed in 4% paraformaldehyde, 2.5% gluteraldehyde in 0.1M sodium cacodylate buffer (pH 7.2) for one hour then cut into small pieces, washed in buffer and post-fixed in 1% osmium tetroxide in cacodylate buffer for one hour or 1% ruthenium tetroxide in cacodylate buffer for one hour. The pieces were dehydrated through alcohol series, into propylene oxide and embedded in Durcupan resin. Ultrathin sections were stained with 3% aqueous uranyl acetate and Reynold's lead citrate and examined on a JEOL 1200EX electron microscope. TEM images were collected on a SIS Megaview III camera.
Functional assessments (i) Trans-epidermal water loss (TEWL) Organoid cultures were equilibrated at room temperature and atmospheric conditions for 30 minutes before TEWL was measured at two locations on the epidermal surface using an Aquaflux AF200 instrument (Biox Systems Ltd, London, UK) with a closed chamber and a custom (5mm diameter) probe head. TEWL measurements were taken every second for a minimum of 60 seconds until a stable reading, as determined by the software, was obtained. TEWL directly measures water evaporation from the skin surface but it may be considered as a measurement of 'inside-outside' barrier function 25 .
(ii) Capacitance To inform the interpretation of TEWL, we measured electrical capacitance; this is directly proportional to water content of the uppermost 10-20µm of tissue. Organoid cultures were equilibrated at room temperature and atmospheric conditions for 30 minutes prior to measurement of epidermal surface capacitance, using a Corneometer TM and multiprobe adapter (CM825 and MPA2, Courage and Khazaka, Cologne, Germany). Three measurements were recorded from each organoid and the mean was calculated.
(iii) Fluorescent dye penetration 50µl of 1mM lucifer yellow dye (Sigma Aldrich) was added to the epidermal surface of the organoid and incubated at 37°C for 4 hours. Metal cloning rings were used to control uniform dosing on the epidermal surface. The lucifer yellow was removed and the organoids washed in PBS before formalin fix-paraffin embedding under standard conditions. 4µM sections were deparaffinized, counterstained with DAPI (1µg/ml for 10 mins) (Life Technologies, Carlsbad, California, USA) and imaged by confocal Zeiss LSM710 microscope. Quantification of dye penetration in the upper dermis (average intensity in upper 40µm) was performed using Zeiss Zen Blue lite version 2.3 software and compared using paired t-tests. Dye penetration represents a measurement of the 'outside-inside' barrier function. The pellet was re-suspended in 200µl of 50mM ammonium bicarbonate. 100µl was taken for in-solution digest. 22µl of dithiothreitol was added and samples heated to 50°C for 15mins. 24µl of iodoacetamide was added and the samples incubated at 22°C for 30 mins. 1µl of RapiGest TM (186001860, Waters, Herts, UK) was added, followed by 2.5µl of trypsin, for a 12-hour digest. A further 1µl trypsin was added for an additional 4-hour digest. The peptide samples were fractionated into 24 fractions with high pH Reversed Phase C18 chromatography, run on a Q Exactive Classic or Plus for 160 mins and the top 15 ions selected for sequencing. Mass spectrometry resolution was 70,000 and tandem mass spectrometry (MS/MS) resolution 17,500. Data were processed using MaxQuant (v1.6.0.13) and Human UniProt Database (Dec 2017) with a protein and peptide false discovery rate of 0.01.

Mass spectrometry proteomic analysis
Proteins identified as contaminants (trypsin and other lab originating proteins) and contained within the MaxQuant database were removed. To account for the possibility of keratin as a contaminant, we analysed the blank runs for keratin peptide intensity and subtracted this from the measured keratin intensities. A reversed database was used to identify false positives. The peptide and protein false discovery rate was set at 1%, and these were removed from further analysis.
Total intensities of the samples were normalized to total protein abundance, with adjustment to the lowest total protein yield. All samples were run in quadruplicate and reproducibility was assessed by Pearson correlation.

Mass spectrometry lipidomic analysis
Organoid samples were sonicated in phosphate buffered saline (pH7.4) and then extracted according to the method of Folch et al. 26 . The samples were analysed by liquid chromatographymass spectrometry (LC-MS) on a Thermo Exactive Orbitrap mass spectrometer (Thermo Scientific, Hemel Hempsted, UK) coupled to a Thermo Accela 1250 ultra high pressure liquid chromatography (UHPLC) system. Samples were injected on to C18 column (Thermo Hypersil Gold, 2.1 mm x 100 mm, 1.9 µm) and separated using a water/acetonitrile/isopropanol gradient 27 . Ceramide 17:0 (Avanti Polar Lipids, Alabaster, AL, USA) was included as an internal standard. Ion signals corresponding to individual ceramide molecular species were extracted from raw LC-MS data sets. Concentration of ceramide was expressed as pmol/mg after normalisation to mg of wet weight tissue.
Data analysis of qPCR, immunoblotting and functional assessments Measurements were compared between non-targeting controltreated organoid models and FLG-siRNA-treated models, using paired t-test.
Histological sections of 10 replicate skin organoid experiments were measured using a standardized technique to minimize bias, as follows: three points were measured on each slide, at the right side, middle and left of each sample, across viable cell layers and stratum corneum; measurements were generated using Fuji ImageJ (May 2017) and the mean of the three measurements calculated.

Proteomic data analysis
The total intensities of the samples were normalised to total protein abundance, with adjustment to the lowest total protein yield. Log2 ratio of protein expression (log2 FLG knockdown -log2 NT control) was visualised as a volcano plot, using GraphPad Prism (version 5) for human proteins detected in (three or all four) of the replicate experiments. P values were calculated by paired t-test applied to log 10 transformed data. Because of the burden of multiple testing in this large dataset, rather than defining arbitrary thresholds for 'statistical significance', our analyses focussed on changes that showed consistency across the biological replicate experiments. Approximately 1000 of the >8000 proteins were selected as showing consistently increased or decreased abundance, defined as follows: Proteins were considered to be reduced in abundance if they showed a reduction in all four biological replicate experiments comparing FLG-knockdown organoid samples to matched non-targeting controls and ratio ≤0.5 in three or four out of four of the replicates and reduced to undetectable levels of protein in a maximum of two replicates. Proteins were considered to be increased in abundance if they showed an increase in all four biological replicate experiments comparing FLGknockdown organoid samples to matched controls treated with non-targeting controls and ratio ≥1.2 in three or four out of four of the replicates. These magnitudes of change were used to define roughly equal proportions of up-and down-regulated proteins for pathway analysis. The groups were then assessed using gene ontology and network analysis in STRING, version 11 28 and pathway prediction using the Reactome Knowledgebase, version 68 29 . This analysis, based on consistency of findings across multiple biological replicate experiments was used in preference to a defined threshold of statistical significance,  with the aim to identify proteins and pathways of importance in the context of a common complex but heterogenous trait.

Lipidomic data analysis
The ratio of ceramide and omega-hydroxy ceramide species were compared between the FLG knockdown organoids and matched control samples in five biological replicate experiments.

Results
All samples are wild-type for the prevalent FLG null mutations and filaggrin shows knockdown at mRNA and protein level Genotyping of DNA from tissue samples from the four skin donors did not detect any of the prevalent FLG null mutations 30 ; this does not, however, exclude the possibility of rare (<0.01 allele frequency) null mutations.
FLG knockdown in the skin organoids persisted to day 10 after siRNA transfection, as demonstrated by qPCR ( Figure 2) and western blotting ( Figure 3). We note the greater magnitude of knockdown quantified at mRNA than protein level in the organoid epidermis 10-days after transfection. This is likely to reflect the inherent differences in mRNA and protein stabilities: FLG mRNA has a half-life of ~2.2 hours 31 whereas the half-lives of profilaggrin and filaggrin are approximately 6 and 24 hours respectively 32 .
FLG knockdown results in morphological changes in the epidermis Histological examination confirmed a reduction of the granular layer, in keeping with a reduction in profilaggrin ( Figure 4A). There was also a slight reduction in thickness of the stratum corneum (ratio 0.83 +/-0.05 (mean +/-SEM), but no difference in thickness of the keratinocyte cell layers ( Figure 4B).
Transmission electron microscopy (TEM) of the skin organoid model shows an effective recapitulation of the physiological layers of keratinocyte differentiation within skin, including a mature, multi-layered stratum corneum with corneodesmosomes and a granular layer with keratohyalin granules and desmosomes ( Figure 5A). TEM of the FLG-knockdown organoid shows a similarly mature epidermis, but the stratum corneum has fewer layers and the granular layer contains keratohyalin     granules that are generally smaller and less numerous than in the control organoid ( Figure 5B).
FLG knockdown shows some evidence of a barrier defect in the immature organoids but no functional effect on skin barrier in the mature organoid model The organoid model described here shows progressive development of skin barrier function after lifting to the air-liquid interface: there is a progressive reduction in TEWL (Figure 6), progressive reduction in capacitance (Figure 7) and the hydrophilic dye, lucifer yellow, is excluded from the epidermis by day 7 in both the control and FLG-knockdown organoids (Figure 8).
Results from five biological replicate experiments showed there was no difference in the mean TEWL (Figure 6), or mean capacitance measurements (Figure 7) between the mature organoid models with FLG knockdown and donor-matched non-targeting siRNA-treated controls. Capacitance was increased in the immature FLG-knockdown model (day 5 after air-exposure, Figure 7) and lucifer yellow dye was not fully excluded by the stratum corneum in some of the the FLG-knockdown replicate experiments until day 7 (Figure 8). However, this delay in maturation of the 'outside-to-inside' barrier function was not consistently demonstrated in all biological replicates.
Mass spectrometry proteomic analysis identifies over 8000 protein species per sample After removal of contaminants and false positive results, >8000 proteins were identified in each epidermal sample. The proteomic data have been deposited to the ProteomeXchange Consortium via the PRIDE 33 partner repository with dataset identifier PXD014875. Quality control analysis confirmed that similar distributions of proteins were identified in each experiment ( Figure 9A); protein extracts were analysed in quadruplicate to test for technical reproducibility and the Pearson correlations were 0.869-0.938 ( Figure 9B). The protein expression changes were visualised using a volcano plot ( Figure 10).  645 of the proteins are identified in STRING and this group showed significant enrichment for interactions (p<1.0e-16) computed by combining the probabilities from the different evidence channels and corrected for the probability of randomly observing an interaction 21 . Gene ontology (GO) analysis in STRING showed evidence for functional enrichment of GOterms including RNA processing (false discovery rate (FDR) 4.40e-11), catalytic activity (FDR 8.59e-08) and intracellular component (FDR 3.88e-26). A list of the 10 most significant GO analysis results in each of the GO domains is shown in Table 1 (see Underlying data) 34 .
Pathway analysis of down-regulated proteins identifies themes including RNA metabolism, adaptive immunity and axon guidance Reactome detected 39 pathways in this set of proteins ( Table 2, Underlying data) 35 . Descriptions include: metabolism of RNA (including 67 out of the 652 genes in this pathway, FDR 3.36e-12); class I MHC mediated antigen processing and presentation (30/365 genes in pathway, FDR 0.0031); and axon guidance (34/541 genes, FDR 0.0207).
Down-regulated proteins indicate molecular mechanisms of relevance to atopic eczema Manual searching of the list of down-regulated proteins identified several of specific interest based on a priori knowledge of eczema pathomechanisms. These highlight several putative therapeutic targets: the aryl hydrocarbon receptor plays a role in signalling pathways for skin barrier repair 22 and has been proposed as a therapeutic target for the treatment of skin inflammation 36 ; caspase 14 plays a role in filaggrin processing 37 and is required for epidermal protection against water loss and UVB-induced damage 38 ; dermokine is a soluble regulator of keratinocyte differentiation 39 ; sequence similarity to a murine protein (UniProtKB:Q8R1R3) indicates that StAR-related lipid transfer protein 7 may play a role in protecting mucosal tissues from exaggerated allergic responses; a reduction in the cytokine TGF-beta 1 may predispose to atopic eczema 40 ; and loss-of-function mutations in CYP4F22 cause lamellar ichthyosis 41 , suggesting that a reduction in CYP4F22 protein expression may contribute to the ichthyotic phenotype observed in ichthyosis vulgaris and atopic skin. 374 proteins are identified in STRING and the group showed significant enrichment for interactions (p<1.0e-16). GO analysis revealed enrichment of terms including translation (FDR 4.12e-19), RNA binding (FDR 9.77e-12) and cytoplasmic part (FDR 9.64-36). A list of the 10 most significant GO analysis results in each of the GO domains is shown in Table 3 (see Underlying data) 42 .
Up-regulated proteins indicate molecular mechanisms of relevance to atopic eczema Manual searching of the list of up-regulated proteins identified several of relevance to atopic eczema, including: STAT1, an intracellular signalling molecule that shows increased expression in atopic eczema, and as part of the JAK-STAT pathway is targeted by several traditional herbal remedies 44 as well as novel small molecule therapies 45 ; carbonic anhydrase 2, which is increased in skin in a variety of forms of eczema 46 ; and FK binding protein, which binds to tacrolimus, an established topical treatment for atopic eczema.
The type I keratins 10, 16 and 17 each show increased expression in the FLG-siRNA-treated organoid models with mean ratio 1.7, 1.3 and 1.3, respectively. Keratin 10 is an intermediate filament protein of structural importance in the skin; keratin 16 is an epidermis-specific keratin that plays a role in innate immunity in response to skin barrier breach 47 ; keratin 17 is expressed in epidermal appendages, including the hair follicle.
Finally, it is noteworthy that GAPDH, which is often considered to be a stable 'housekeeping' gene, is in the dataset of consistently up-regulated proteins. However, in the context of this work, filaggrin knockdown is clearly apparent at a greater magnitude than the more subtle change in GAPDH (Figure 3).
Epidermal ceramide content shows no significant difference in FLG-knockdown organoids compared with controls Ceramide and omega-hydroxy ceramide species in five biological replicate experiments showed no consistent difference between the FLG knockdown organoids and matched control samples (Figure 11).

Discussion
The functional and biochemical analyses of a filaggrin-deficient skin model have provided findings that are consistent with our and others' previous analyses in vitro 19 and in vivo 48 . We have included both epidermal and dermal compartments in this model (to allow for dermal-epidermal cross-talk 49,50 ) and our more detailed proteomic analysis has revealed additional evidence of molecular mechanisms with relevance to atopic skin. These datasets are shared as a resource for the research community.
Analysis of the immature organoid identified a delay in epidermal maturation following FLG knockdown, but the lack of a replicable functional effect in our model, when mature, reflects previous reports in which epidermal equivalents created from FLG-null keratinocytes do not show impairment in barrier function 20 . This is in keeping with the observation that FLGnull mutations are not fully penetrant with respect to ichthyosis vulgaris or atopic eczema 16,51 .
Because of the nature of our data, generated using biological replicate samples, and the focus on modelling a common complex trait with considerable inter-individual variation, we elected to use the criterion of consistency across the biological replicates to define proteins of interest, followed by network and pathway analyses. We recognise that other, more complex data analytical approaches may be applied, and/or more focussed interrogation of specific pathways; we share the proteomics data and welcome such further analysis.
Our analysis has highlighted aspects of keratinocyte-immune signalling. This is an important component of the immunological barrier that may be overlooked in tissue with multiple different cell types. The predicted pathways entitled 'innate immunity' and 'adaptive immunity' included proteins showing both increased and decreased expression (Table 2 and Table 4, Underlying data) 35,43 . This immune dysregulation is likely to contribute to the predisposition to infection and inflammation observed in filaggrin-deficient skin in vivo: filaggrin-deficient patients show increased viral infection of the skin 52 and increased incidence of irritant 53,54 and allergic contact dermatitis 55 . The regulation of innate and adaptive immune pathways within keratinocytes (in addition to cells of the haematopoietic lineage) may therefore represent a worthwhile target for future therapeutic interventions.
The lack of effect of FLG knockdown on the ceramide content in the model epidermis is also consistent with previously reported findings in vitro 56,57 and in vivo 58 . We did not analyse in detail the other epidermal lipids and therefore we cannot comment on whether there were relevant changes. Differences in lamellar body structure 59 have been reported in FLG-mutant skin biopsies but this was not clearly observed in the organoid model. However, ultrastructural examination of the model stratum corneum did reveal abnormalities in the lamellar bilayer structure and less prominent corneodesmosomes ( Figure 5), as previously reported in vivo 59 . It is also noteworthy that pathway analysis in Reactome identified the terms 'metabolism of lipids' and 'phospholipid metabolism' from the group of consistently down-regulated proteins (Table 4, Underlying data) 43 .
The most significant themes identified in the groups of proteins showing increased expression upon FLG knockdown relate to RNA binding and translation; this is in contrast to RNA metabolism and RNA processing, which are demonstrated in the proteins showing decreased expression. The functional significance of this increase in transcriptional and translational activity cannot be determined from the gene ontology and pathway analysis. We may hypothesise that it reflects the increased 'stress' response, as we have previously observed in FLG genotype-stratified transcriptome analysis of atopic skin biopsies 18 and/or the activation of innate and adaptive immune mechanisms.
It is tempting to speculate that the increased expression of keratins 10 and 16 may contribute to palmar hyperlinearity and keratin 17 may contribute to keratosis pilaris. These are characteristic features of ichthyosis vulgaris but their pathogenesis has not yet been explained as a direct result of filaggrin haploinsufficiency 16,60 .
Several of the proteins showing reduced expression relate to aspects of differentiation, including filaggrin itself, and caspase 14 which plays a role in the terminal degradation of filaggrin 37 . The aryl hydrocarbon receptor (AhR) is a ligandactivated transcription factor that responds to multiple different environmental and metabolic stimuli to control transcriptional pathways of relevance to atopic eczema, including immunity and differentiation 61 . Activation of AhR reduces skin inflammation 62 and enhances barrier repair 63 ; therefore, a reduction in AhR expression in the FLG knockdown organoid would be anticipated to increase the eczema diathesis by opposing these effects. AKT1 activity is known to be reduced in atopic skin, leading to an alteration in protease expression and reduced filaggrin expression and processing 64 . This mechanism may further exacerbate the filaggrin deficiency in our FLG knockdown organoid. Dermokine acts as a soluble regulator of keratinocyte differentiation and mice deficient in the epidermal dermokines (beta and gamma isoforms) show defective cornification 65 . Transforming growth factor (TGF) beta-1 is a multifunctional protein: it regulates the growth and differentiation of multiple cell types, including keratinocytes and it can promote either Th17 or T-regulatory cell lineage differentiation in a concentration-dependent manner; it appears to have an immunosuppressive effect in atopic disease since a geneticallydetermined lower production of TGF-beta-1 is associated with increased risk of atopic dermatitis 40 .
The Reactome term 'axon guidance' was detected in proteins showing increased and decreased expression. It has long been recognised that there is an increase in the density of sensory neurons in eczematous skin 66 and nerve growth factor expression is upregulated in the keratinocytes of patients with atopic eczema 67 . Therefore, dysregulation of axon guidance may be of relevance to the pruritus that is so characteristic of atopic eczema and contributes substantially to the morbidity of this condition 3 . Antagonists of transient receptor potential vanilloid subfamily member 1 (TRPV1), expressed on sensory nerves and keratinocytes, are now being investigated for the control of pruritus 68,69 .
Taken together, our data demonstrate that filaggrin haploinsufficiency leads to abnormalities in both structural and immune features within the keratinocyte compartment of the epidermis. Keratinocytes display immune activity, producing multiple cytokines and chemokines in addition to their roles in innate immunity 70,71 . Together with the evidence of dysregulation in axon guidance, these findings indicate that patients with filaggrin deficiency -either haploinsufficiency or a reduction in filaggrin expression that occurs secondary to atopic inflammation 72 -may derive clinical benefit from future therapies targeting keratinocyte-immune crosstalk mechanisms and neurogenic pruritus. The skin organoid model described here represents a valuable opportunity to study genetic effects on keratinocytes in relative isolation. However, the lack of other cell types (notably T cells, B cells, dermal dendritic cells and innate lymphoid cells), the absence of complex dermal structures (including neurons and blood vessels) and skin appendages (for example hair follicles, sweat and sebaceous glands) is also a limitation in that it precludes assessment of the complex interactions that are likely to occur between different tissue compartments in skin.
In conclusion, we share these experimental results in detailed form, because of the additional insight provided by this in vitro analysis of filaggrin deficiency in a skin organoid model; our findings compliment and extend previous work in vitro and in vivo. Our analyses have re-emphasised the role of keratinocyte-specific mechanisms in contributing to immune and neurological as well as structural aspects of skin barrier function and offer new insight into molecular mechanisms for the ichthyosis vulgaris phenotype. Elias . investigated the effect of filaggrin deficiency in an optimized skin organoid model. Next to the et al inside/out and outside/in barrier function assessed respectively by TEWL and penetration of a hydrophilic dye, a global mass proteomic analysis has been done. As discussed by authors, the used model has important advantage that the changes related to filaggrin deficiency could be studied in specific cells (keratinocytes) and both epidermal and epidermal compartments were present. On the other side, the model misses other relevant cell types. The study is well designed and presented.

Skin barrier function assessment
Although Lucifer yellow is a commonly used model penetrant it has to be realized that owing to its hydrophilicity and a large molar mass (444 Da) it has quite low penetration rate and a long lag time. The incubation time of 4 hours before microscopic examination is likely too short to detect penetration across the SC unless the barrier damage is substantial. In general, it would be good to include in functional studies both, hydrophilic and lipophilic penetrants and measure their penetration at different time points. This might provide more detailed insight into alterations in the structure and composition of the SC and relate them to the structural changes. As noted by the authors, the Flg-knock down ''reveal various abnormalities: in the lamellar bilayer structure and less prominent corneodesmosomes'' and furthermore down-regulation of proteins involved in metabolism of lipids' and 'phospholipid metabolism'' and thinner stratum corneum. Though, no changes in either TEWL or dye penetration have been found. This model would be of great value to further study functional consequences of filaggrin deficiency, including previous damage of the skin e.g. by skin irritants.

Other points
The authors note that ''FLG knockdown in the skin organoids persisted to day 10 after siRNA transfection, as demonstrated by qPCR ( Figure 2) and western blotting (Figure 3)''. In Fig. 2 it is obvious that the knock-down at the RNA level was not complete and lower than achieved in similar studies (

Are sufficient details of methods and analysis provided to allow replication by others? Yes
If applicable, is the statistical analysis and its interpretation appropriate? I cannot comment. A qualified statistician is required.
Are all the source data underlying the results available to ensure full reproducibility? Yes

Are the conclusions drawn adequately supported by the results? Yes
No competing interests were disclosed.

Competing Interests:
Reviewer Expertise: Skin barrier in inflammatory skin diseases, dermatotoxicology.

I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.
Author Response 19 Nov 2019 , Ninewells Hospital, Dundee, UK Sara J. Brown

Skin barrier function assessment
We agree that Lucifer yellow has a slow penetration rate hence our long incubation (4 hours, compared to some other publications reporting 1 hour or less). We will, in future work, consider use of a hydrophobic dye as an additional, complimentary assessment -thank you for this useful suggestion.

Other points
The mRNA quantification was carried out 10 days after siRNA transfection and we agree it is interesting to note this decay in knockdown at transcriptome level, whilst the protein knockdown persists to a greater degree. This is likely to reflect the different half-lives of mRNA FLG (approximately 2.2 hours ) whereas the half-lives of profilaggrin and filaggrin are approximately 6 and 24 hours respectively. We have added this additional explanation to the updated manuscript.
Thank you for the knowledgeable advice about ceramide analysis. We agree that the stratum corneum is the most relevant tissue layer in which to perform lipidomic analysis, but an accurate separation of the whole stratum corneum (without upper epidermal cells) proved to be technically challenging in the organoid model. We therefore chose to analyse full thickness epidermis to 1 2 1.

4.
the prior proteomic analysis published on filaggrin knockdown keratinocytes (Elias ., 2017 ), et al as that model did not include matched donor fibroblasts. Table 2 presents Reactome data and highlights Axon guidance, and it would be helpful to see this in the context of the other pathways described. I would favour this being included in the PDF, rather than Fig 9, which may be moved to supplementary data.
The volcano plot (Fig 10) is of interest, but more data could be presented here, for example, the 7 points which are log2 fold change >+5 could be identified to make this more informative for the reader. I am uncertain if the colour scheme with p value cut offs is helpful; the authors carefully justify using selected fold change thresholds in the methods section over arbitrary thresholds that are p-value driven and this depiction is at odds with this.
Could the title be improved to reflect the more complete nature of this skin equivalent model that incorporates fibroblasts? Also, how were the three elements "increased transcriptional-translational activity, keratinocyte-immune crosstalk and disordered axon guidance" selected or prioritised from the list of Reactome findings? There is a minor typographical error on page 3. "Medium was refreshed on alternate day the cultures sand were used for analysis up to day 12 (representing three, five and seven days after lifting to the air-liquid interface)".

If applicable, is the statistical analysis and its interpretation appropriate? Yes
Are all the source data underlying the results available to ensure full reproducibility? Yes

Are the conclusions drawn adequately supported by the results? Yes
No competing interests were disclosed.

Competing Interests:
Reviewer Expertise: Dermatology, Genetics, Primary cell culture models.

I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.
Author Response 19 Nov 2019 , Ninewells Hospital, Dundee, UK Sara J. Brown

Comments:
The authors highlight the proteins that have relevance to the known pathogenetic mechanisms of eczema. It would be helpful to emphasise the novel findings detected in this model compared to the prior proteomic analysis published on filaggrin knockdown keratinocytes (Elias et al., 20171), as that model did not include matched donor fibroblasts.
Thank you for this insightful comment. The current work differs in several ways from the previous work (performed by the same first author in different research institutions): (i) The current work used FLG wild-type primary keratinocytes isolated from female breast skin and co-cultured with 3T3 feeder cells in the presence of serum, whereas the previous publication used keratinocytes isolated from foreskin donors and cultured in serum-free conditions. (ii) The current work used a full-thickness skin organoid model ((Keratinocytes seeded on top of donor matched fibroblasts embedded in a fibrin dermal matrix) whereas the previous work was anepidermal equivalent model (Keratinocyes grown in a transwell) .
(iii) The current work benefited from advancements in proteomic technology and generated a more extensive, detailed proteomic dataset.
The previous work provided, for the first time, evidence that filaggrin deficiency can alter the expression levels of proteins relevant to the pathogenesis of AE, even without an external inflammatory stimulus and the current work supports this.
Mechanistic themes emerging from this study, including inflammatory, proteolytic and cytoskeletal pathways, broadly overlap with findings from this previous work. However, novel findings detected in the current model compared to the prior proteomic analysis (Elias et al., 2017) may be summarised as follows: Evidence that filaggrin-deficient keratinocytes are a source of immune signalling. Filaggrin deficiency leads to an increase in transcriptional and translational activity but a reduction in multiple markers of terminal differentiation.
There is a dysregulation of axon guidance (with elements of this pathway represented in both the up-and down-regulated protein sets). Thank you for this suggestion. We agree that Table 2 is a valuable summary of the findings from Reactome, but it is too large (including detailed lists of protein IDs, which we are keen to include) to fit into the .pdf. We hope that readers will easily access the full table via the direct link provided to FigShare. We have considered alternative methods to present our Reactome pathway analyses, but the visualisation of pathways within a global over-view does not lend itself to our data where relatively small numbers of pathways (39 and 113) show significant dysregulation.
The volcano plot (Fig 10) is of interest, but more data could be presented here, for example, the 7 points which are log2 fold change >+5 could be identified to make this more informative for the reader. I am uncertain if the colour scheme with p value cut offs is helpful; the authors carefully justify using selected fold change thresholds in the methods section over arbitrary thresholds that are p-value driven and this depiction is at odds with this.
Thank you for this advice. On further consideration we agree that the colour-coding is not helpful and the addition of labels to the most extreme data points is informative for readers. These changes have been made in a corrected version of Fig 10 in the updated manuscript.
Could the title be improved to reflect the more complete nature of this skin equivalent model that incorporates fibroblasts? Also, how were the three elements "increased transcriptional-translational activity, keratinocyte-immune crosstalk and disordered axon guidance" selected or prioritised from the list of Reactome findings?
In response to your helpful suggestion we have added the words 'full thickness' to the title to emphasise the nature of our model. Because the title was already very long we have truncated it, since the main findings ('increased transcriptional-translational activity, keratinocyte-immune crosstalk and disordered axon guidance') are clearly stated in the abstract.
These three elements were prioritised as noteworthy findings based on a combination of GO and pathway analyses as follows: (a) Aspects of transcriptional and translational activity are represented in all 10 of the most significant GO terms for Biological process in the consistently down-regulated (Table 1) and up-regulated (Table 3) proteins; 'Metabolism of RNA' is the top Reactome term in the down-regulated protein dataset (Table 2) and the Reactome term 'Translation' is top in the upregulated proteins (Table 4).
(b) Keratinocyte-immune crosstalk is highlighted to emphasise that these experiments, in which epidermal keratinocytes are isolated from other cells of the immune system, still generate signals classified within the terms 'innate' and 'acquired immunity'. In the Reactome analyses of consistently down-regulated proteins 7/39 enriched pathways are related to the immune or antiviral systems and 15/113 enriched pathways in the up-regulated proteins relate to immune systems or antiviral systems.
(c) Disordered axon guidance was represented by only one Reactome term in each of the up-and down-regulated protein groups, but together these account for 74/541 (14%) of the proteins in this pathway. This was prioritised as a novel signal of interest, likely to be relevant to AD pathophysiology but not currently targeted for therapy.
There is a minor typographical error on page 3. "Medium was refreshed on alternate day sand the cultures were used for analysis up to day 12 (representing three, five and seven days after lifting to the air-liquid interface)".
Thank you, we have corrected this typo in the updated version of our manuscript.
No competing interests were disclosed. Competing Interests:

6.
Benjamini-Hochberg? Did the authors try e.g. a limma moderated t-test that has proven useful for proteomics data? Even without multiple testing correction, a statistical test appears more suitable than selecting differentially abundant proteins only based on consistency between comparisons/experiments. Moreover, the authors apply a cut-off of 2-fold (0.5) to identify proteins with lower abundance upon filaggrin knockdown and of 1.2-fold (1.2) for those with higher abundance in this condition. It appears this skewed effect size has been selected aiming at equalizing hit numbers for pathway enrichment. This should be reconsidered. The terms 'ratio' and 'fold-change' are not used properly, since the authors report actual ratios as fold-changes, e.g. a ratio of 0.5 would correspond to a fold-change of -2 and of 2 to a fold change of +2. Since identifying proteins with statistically significantly differential abundance between conditions is crucial, also for follow-up data analysis, the authors should re-evaluate this part of the study with help of an expert in proteomics statistics.

If applicable, is the statistical analysis and its interpretation appropriate? Partly
Are all the source data underlying the results available to ensure full reproducibility? Yes Thank you for your interest in the skin organoid model and our analyses. We believe that a focus on the keratinocyte phenotype is valuable, particularly in AD where the keratinocyte's contribution to barrier function is of particular interest. Whilst we agree that there are many other cell types contributing to AD pathogenesis, we consider that the isolation of keratinocytes may be viewed as an advantage in this experimental system because it provides an opportunity to investigate the specific contribution of keratinocyte biology to features of atopic inflammation and other pathomechanisms in this highly complex trait. pathomechanisms in this highly complex trait.
We agree that functional follow-up may be an important step and we have conducted extensive validation in other applications of the skin organoid and a simpler epidermal model. The detailed dataset and our analyses reported here in are shared as a resource for Wellcome Open Research the community, to allow others to review the differentially expressed proteins and follow up selected candidates as indicated. We have further clarified this as the principle aim of our publication in Version 2.

Specific comments:
The filaggrin knock-down at the RNA level was only approximately 50%, whereas the knock-down at the protein levels seems more robust. How do the authors explain this inconsistency?
We agree, there was a 0.5-fold reduction in mRNA compared to 0.36-and 0.25-fold reduction FLG in pro-filaggrin and filaggrin proteins respectively. These differing levels of reduction in expression reflect the use of siRNA which creates a transient knockdown. siRNA is our preferred methodology because we aim to maintain the primary cells at low passage and to minimise non-physiological experimental influences. The different degrees of knockdown reflect the long experimental time points (10 days post transfection for a mature organoid) and the inherent differences in mRNA and protein stabilities.
mRNA has a reported half-life of ~2.2 hours whereas the half-lives of FLG profilaggrin and filaggrin are approximately 6 and 24 hours respectively. We have added this additional explanation to the updated manuscript. Since our main experimental aim was to reduce filaggrin at the protein level, we are satisfied that knockdown of the functional monomer is maintained for the duration of our experiment and this degree of reduction is in keeping with the reduction in expression observed in -null heterozygotes or atopic skin. FLG Why are there no error bars in Fig. 2 NT-siRNA?
The filaggrin expression presented in Figure 2 has been calculated in a pair-wise manner (FLG -siRNA versus -siRNA control) for each independent biological donor (N=8) and the replicate NT data were then combined for statistical analysis. Since only one technical replicate was performed per donor, the -siRNA condition is calculated as 1 in each biological replicate and hence NT displays no error in the combined dataset. This internally controlled method of analysis accounts for differences in baseline filaggrin expression and technical variations between replicate experiments.

The data should be confirmed by immunostaining of the cultures for filaggrin.
Thank you for this suggestion. We agree that immunostaining can represent a valuable confirmatory approach, however the H&E stained section show a clear reduction in the granular layer and the filaggrin knockdown at protein level was confirmed by mass spec proteomics (mean 0.27-fold, N=4) which is likely to be more sensitive for the detection of quantitative differences than immunostaining.
Staining of the cultures for loricrin would be helpful to confirm abnormalities in the granular layer.
The electron microscopy ( Figure 5) shows the granular layer in some detail, with a clear reduction in the number and size of granules. We therefore did not elect to stain the organoid cultures for loricrin specifically. It is impossible to see obvious differences in Lucifer yellow penetration and there is also no statistically significant difference. This should be expressed more carefully.
Thank you for bringing this to our attention. We did observe a difference in Lucifer yellow dye penetration in 3/5 of the biological replicates at day 3 but we agree this is not apparent in the mature model at day 7 (example images shown in Figure 8A). We aimed to explain that our quantification (in Figure 8B) shows all of the individual data points and, whilst some experiments showed an increase in dye penetration in the immature organoid, we state "However, this delay in maturation of the 'outside-to-inside' barrier function was not consistently demonstrated in all biological replicates." In light of your comment we have added further clarification to the text.
The proteomics data provide a very important resource. It would have been helpful to show the validation of some of the data using independent samples and another method, e.g. Western blot and/or immunohistochemistry. This would be easy for selected proteins, such as the different keratins.
Thank you, we agree the proteomics data are a valuable resource and we are delighted to share this resource. We also agree it will be relatively easy for other groups to validate proteins of interest in their work; therefore, rather than carrying out validation ourselves for selected keratins (or multiple different proteins as we have for our parallel proteomics analyses, recently published ) we prefer to share the data for further independent validation.

A comparison of the data with available transcriptomics data from AD patients would be interesting.
We agree that this would be a very interesting analysis; it is not within the remit of this study but we or others may pursue this in the future.
For the proteomics analyses, it is not clear how the false discovery rate was exactly determined (scrambled or reversed database?) and how false positives were removed. This requires clarification.
We apologise for the lack of clarity in this description. The false discovery rate was determined using a reversed database. The peptide and protein false discovery rate was set at 1%, and these identified peptides were denoted with REV_ and removed from further analysis. These details have been added to the revised manuscript.
Statistical analysis to identify proteins with differential abundance between conditions lacks stringency and appears asymmetric with regard to positive and negative effect sizes. Lack of statistically significant hits upon multiple testing correction is a common problem with proteomics datasets. Which multiple testing correction algorithm was applied? Bonferroni or Benjamini-Hochberg? Did the authors try e.g. a limma moderated t-test that has proven useful for proteomics data? Even without multiple testing correction, a statistical test appears more suitable than selecting differentially abundant proteins only based on consistency between comparisons/experiments.
We appreciate that there are no clear guidelines on methodology for the analysis of large proteomics datasets (or indeed many other large datasets) and different expert opinions exist regarding the most appropriate tools and thresholds to define 'significance'. We have benefited from the considerable expertise within our multidisciplinary collaboration, leading to the rather 1 from the considerable expertise within our multidisciplinary collaboration, leading to the rather simple data analysis that we have performed, whilst providing the full dataset for others to analyse as they may wish.
A full Bonferroni correction does not identify any differentially expressed proteins because of the large number detected in this global analysis. Because of the nature of our data, generated using biological replicate samples (primary cells from 4 separate human donors), and the focus on a common complex trait, we prefer to use the criterion of *consistency* across the biological replicates, with filtering for fold change to provide similar numbers of up-and down-regulated proteins, followed by network and pathway analyses.
Moreover, the authors apply a cut-off of 2-fold (0.5) to identify proteins with lower abundance upon filaggrin knockdown and of 1.2-fold (1.2) for those with higher abundance in this condition. It appears this skewed effect size has been selected aiming at equalizing hit numbers for pathway enrichment. This should be reconsidered.
Thank you for highlighting this question. We have stated clearly and openly that we select did these cut-offs aiming to equalize protein numbers in this right-skewed dataset. Using the thresholds of £0.5 and ³1.2 defines 1422 (~25%) and 1871 (~34%) of the total number of proteins. This is an established approach in proteomic analysis; we appreciate other valid approaches can be applied and we welcome further interrogation of our experimental data. We have further clarified this in the discussion.
The terms 'ratio' and 'fold-change' are not used properly, since the authors report actual ratios as fold-changes, e.g. a ratio of 0.5 would correspond to a fold-change of -2 and of 2 to a fold change of +2.