Genomic and physiological mechanisms underlying skin plasticity during water to air transition in an amphibious fish

Devonian fishes seeded the first successful vertebrate invasion of land, where various phenotypic innovations promoted a dramatic terrestrial radiation. What were the mechanistic changes that facilitated this pivotal adaptive leap to life on land? It is difficult to infer the physiological traits that enabled the original Devonian invasion of land. However, modern fishes have repeatedly evolved physiological and behavioral traits that enable exploitation of terrestrial niches (Graham, 1997). Extant amphibious fishes offer an opportunity to study the variety of mechanisms that may have facilitated the challenging transition of vertebrate lifestyles between aquatic and terrestrial habitats.

The skin of amphibious fishes plays a critical role in maintaining homeostasis during air exposure. Histological studies have shown that a number of species have capillaries close to the epidermal surface (Grizzle and Thiyagarajah, 1987; Park et al., 2000; Yokoya and Tamura, 1992), providing an effective surface for gas exchange (e.g. Bridges, 1988; Graham, 1973; Martin, 2014; Tamura et al., 1976). As well, ionocytes normally found in the gills are positioned close to the epidermal surface in some amphibious fishes (Leblanc et al., 2010; Martin, 2014; Sturla et al., 2001; Yokoya and Tamura, 1992) and appear to be important in ion regulation (Stiffler et al., 1986; Yokota et al., 1997). Cutaneous water transport, mucus production and resistance to evaporative water loss likely also contribute to water balance during prolonged air exposure (Dabruzzi et al., 2011; Heffell et al., 2018; Wilkie et al., 2007). In addition, the skin is typically a site of nitrogen excretion in fish out of water (Davenport and Sayer, 1986; Frick and Wright, 2002; Hung et al., 2009; Livingston et al., 2018; Souza-Bastos et al., 2014). But the skin must provide protection from pathogens and environmental hazards (Glover et al., 2013). Thus, the skin is multifunctional in amphibious fishes, but the mechanisms that resolve the opposing pressures of ‘exchange’ versus ‘barrier’ have received little attention.

For the skin of amphibious fishes to become a crucial organ for oxygen transport, ion and water regulation, and nitrogenous waste excretion, and to provide a mechanical barrier out of water, we hypothesized that structural, physiological and molecular changes must occur soon after air exposure. We therefore predicted that rapid and sustained regulation of molecular pathways associated with dermal/epidermal structural remodeling, gas exchange, ammonia detoxification, and water and ion transport underpin the alterations in skin structure and function that unfold during this transition from immersion to emersion.

To test these hypotheses, we used the amphibious mangrove rivulus [Kryptolebias marmoratus (Poey 1880)], which is a particularly opportune model system for studying the physiological and genomic mechanisms that support terrestrial acclimation. The species is native to the Western Atlantic mangrove forests, from southern Florida to Brazil. They frequently escape water for terrestrial microhabitats, including the mud walls of crab burrows, damp leaf litter or rotting logs, for variable durations up to 66 days (Taylor, 2012).

We compared two clonal strains of K. marmoratus during acclimation from water to air. One lineage is native to brackish water habitats in Honduras (HON) and the other to freshwater habitats in Belize (FW). We chose these lineages because of their different osmotic niches may contribute to physiological differences in water and ion regulation, and because prior studies showed that the HON strain had a consistently higher resilience to air exposure relative to the FW strain (T.S.B. and P.A.W., personal observations). Lineages were independently perpetuated by self-fertilization, resulting in a purging of within-lineage genetic variation and creation of isogenic strains (Tatarenkov et al., 2010). Strains were acclimated to their native salinities, and then over the course of 1 week of air exposure we measured changes in skin genome-wide gene expression (RNA-seq), skin structure (light microscopy, transmission electron microscopy) and survival to infer skin mechanisms that resolve the multiple competing physiological demands of terrestrial acclimation. We predicted that genotype and osmotic acclimation would interact to provoke between-strain differences in the regulation of some of these mechanisms, and thereby provide additional insight into mechanisms that underlie variable terrestrial acclimation abilities. Finally, we used models of comparative protein evolution to identify genes evolving by directional selection in K. marmoratus compared with other fish species that have limited or no terrestrial acclimation ability. We predicted that genes evolving by directional selection uniquely in the K. marmoratus lineage would be associated with the molecular pathways that regulate the species’ unusual terrestrial acclimation abilities.

Two isogenic hermaphroditic strains of K. marmoratus were maintained individually in the Hagen Aqualab at the University of Guelph, Ontario. A strain originally collected from Honduras (HON) in 1996 was held in the laboratory (15 ppt, 25°C) for 20 years (∼60 generations; HON11; Tatarenkov et al., 2010). A freshwater strain (FW) captured in December 2012 at Long Caye, Belize, in an open pool (0.3 ppt) was held in the laboratory (0.3 ppt, 25°C) for 3.5 years (∼10 generations; Platek et al., 2017). HON fish and those from the FW site in Belize are genetically distinct (Lins et al., 2018). Fish were held singly in plastic containers (120 ml Fisherbrand; Thermo Fisher Scientific, Markham, ON, Canada) in 60 ml of water under constant conditions (12 h:12 h light:dark cycle, 25°C; Blanchard et al., 2019). Reverse osmosis water and sea salt (Instant Ocean, Crystal Sea, Baltimore, MA, USA) were used to formulate water to the appropriate salinity and changed weekly. Fish were fed live Artemia three times a week. The University of Guelph Animal Care Committee approved this project (AUP 2239).

HON (0.120±0.005 g) and FW (0.130±0.006 g) fish were air-exposed for up to 7 days at 25°C. Before aerial transfer, tissues from immersed fish were sampled (control, 0 h). Upon air exposure, fish were maintained on moist filter paper (15 or 0.3 ppt) in plastic rearing containers (99% humidity) and skin samples were collected at 1, 6, 24, 72 and 168 h for transcriptomics. Fish were euthanized with tricaine methanesulfonate (1.5 mg ml⁻¹). Skin was dissected from the mid-region of the body by making a ventral incision and carefully removing extraneous tissue and muscle from the skin, and was then stored in RNAlater (Thermo Fisher Scientific) at −20°C. For histology, fish were acclimated to water or air (0 h and 168 h), euthanized by immersion in ice-cold water, and then fixed and preserved for either collagen staining or transmission electron microscopy (TEM) imaging.

Survival was measured using a different set of fish (HON, n=20; FW, n=20) that were air-exposed for up to 7 days and monitored daily. The experiment was terminated if survival reached as low as 20%.

We measured skin collagen content because it is a key contributor to the mechanical properties of skin especially in amphibious species (Craig et al., 1987; Schwinger et al., 2001), and collagen is regulated in other tissues (e.g. gill arches and filaments) during air acclimation in K. marmoratus (Turko et al., 2017). For collagen staining, whole fish were fixed and processed, as described previously (Turko et al., 2014). Tissues were then stained with Picrosirius Red, a collagen-specific dye (Electron Microscopy Sciences, Hatfield, PA, USA), photographed under circularly polarized light (Johnson et al., 2014), and the mean collagen hue for each picture was calculated (Turko et al., 2017).

Samples for TEM observation were prepared using skin sampled from both HON and FW animals that were either immersion controls or had been exposed to air for 7 days. Samples were acquired from a standardized region of skin, just posterior to the operculum and dorsal of the horizontal midline. Tissue was removed to include a small section of underlying epaxial muscle so as not to damage the skin, and to provide orientation when sectioning. Isolated tissue was placed immediately into freshly prepared, ice-cold glutaraldehyde fixative (2.5% glutaraldehyde in 0.1 mol l⁻¹ sodium cacodylate buffer, pH 7.3). Samples were fixed at 4°C overnight and then rinsed twice in 0.1 mol l⁻¹ sodium cacodylate buffer (2×10 min, pH 7.3). Samples were post-fixed in osmium teroxide (1% osmium teroxide in 0.1 mol l⁻¹ sodium cacodylate buffer, pH 7.3), rinsed twice in buffer, and dehydrated through a graded ethanol series (70%, 90% and 3×100%, 10 min each), then immersed in propylene oxide (2×10 min) and transitioned into Quetol-Spurr resin (2 h in 1:1 propylene oxide/Quetol-Spurr resin followed by 2 h in Quetol-Spurr resin). Final infiltration in Quetol-Spurr resin was carried out overnight, after which sample pieces were placed in embedding molds and polymerized at 65°C overnight. Sections (90 nm) were prepared using a Leica EM UC7 ultramicrotome, stained with uranyl acetate and lead citrate, and viewed in an FEI Tecnai 20 transmission electron microscope. Skin morphometric analysis was conducted using ImageJ software (https://imagej.nih.gov/ij/). Sub-apical desmosome 2 was the second desmosome observed in the sub-apical region of the epidermis, and its depth was determined by tracing its length. On average, for each fish, the lengths of typically six to 10 desmosomes were measured and then averaged to give an N of 1 for that animal. A total of 158 desmosomes were measured for 20 animals.

Transcriptomics data collection and analyses were described in Blanchard et al. (2019). Briefly, tissues were homogenized and RNA was extracted using RNAeasy purification kits (Agilent, Inc.). Individually indexed RNA-seq libraries were prepared using NEBnext RNA library preparation kits for Illumina, libraries were pooled, and sequence data were collected across four lanes of Illumina 4000 (PE-150). Five biological replicates were included per treatment group, except for two treatments (HON strain 72 h air exposure, and HON strain immersion control) that had four replicates. Sequence reads were mapped to the K. marmoratus gene set reported in Rhee et al. (2017) (Whole Genome Shotgun project GenBank accession LWHD00000000, GenBank assembly accession: GCA_001649575.1) using STAR v2.4.2a (Dobin et al., 2013), and read counts were generated using HTseq v0.6.1 (Anders et al., 2015). We retained genes that had greater than 10 read counts for at least four of five biological replicates within any treatment group. Read counts were log₂ transformed and normalized for gene length and total library size in EdgeR v3.20.9 (Robinson et al., 2010). Differential expression analysis was performed in Limma v3.34.9 (Ritchie et al., 2015), where strain (HON, FW) and time (0, 1, 6, 24, 72 and 168 h) were specified as fixed main effects, and strain×time as an interaction term, including false discovery rate adjustment of P-values. RNA-seq read sequences are archived at NCBI (SRA accession: SRP136920, BioProject: PRJNA448276). A matrix of per-gene per-sample raw read counts (‘final.out.txt’) and R scripts that detail all of these analyses are available in a GitHub repository (https://github.com/WhiteheadLab/mangrove_killifish). A list of all genes included in the analysis, including annotation, average expression levels and P-values associated with results from statistical tests, is in appendix S1, available from figshare (https://doi.org/10.6084/m9.figshare.13344806.v1).

Positively selected molecular changes that are lineage-specific underlie adaptive phenotypes that may be unique to a species. They leave recognizable footprints in the genome, for example in the amino acid sequence of proteins. Such footprints are detected when sequences are compared between orthologous loci among relatives. We used OrthoMCL (Li et al., 2003), which is a program that filters orthologous proteins, and GUIDANCE2 (Sela et al., 2015) to remove unreliable coding sequence (CDS) alignments and alignment positions. We then used phylogenetic analysis by maximum likelihood (PAML; Yang, 2007) to identify positively selected sites on targeted lineages, in our case those that are specific for K. marmoratus (see fig. S1, available from figshare at https://doi.org/10.6084/m9.figshare.13344806.v1). Because K. marmoratus is the only amphibious species among those included in this analysis, the genetic changes that support their unique amphibious lifestyle are likely among those showing positive selection specifically in the K. marmoratus lineage.

We downloaded protein and CDS sequences in FASTA format of genome assemblies from: (1) K. marmoratus (mangrove rivulus), (2) Austrofundulus limnaeus (bony fishes), (3) Nothobranchius furzeri (turquoise killifish), (4) Poecillia formosa (Amazon molly), (5) Fundulus heteroclitus (mummichog), (6) Oryzias latipes (Japanese medaka; ASM223467v1), (7) Gasterosteus aculeatus (three-spined stickleback), (8) Takifugu rubripes (fugu) and (9) Danio rerio (zebrafish) from NCBI and the Broad Institute (Fig. S1). We chose these species because they had well-assembled and -annotated genomes, they are all strictly aquatic (not amphibious) species and they represented a range of phylogenetic distances from K. marmoratus. To assign a database of orthologous sequences from these nine sets of genome protein FASTA sequences, we used OrthoMCL v2.0.9 (Li et al., 2003), and the following criteria were applied to achieve optimal clustering resolutions: similarity P-value ≤1×10^–5, protein percent identity ≥40% and MCL inflation of 1.5. In a preliminary clustering step, compliant FASTA files were established by applying the module orthomclAdjustFasta. The OrthoMCL module orthomclFilterFasta was used to remove proteins of poor quality to finally generate the goodProtein.fasta-file compiling all proteins (337,958 proteins; 2301 MB) with specific assigned headers for the subsequent all-versus-all BLAST.

To build orthologous groups, we first conducted an all-versus-all NCBI BLAST (ncbi-blast-2.2.24+). With the OrthoMCL module orthomclBlastParser, similarities were then summarized within a file that was loaded to mysql by the orthomclLoadBlast module. The orthomclPairs and orthomclDumpPairsFiles modules then found, sorted and created results files of all protein pairs (orthologs, in-paralogs, co-orthologs). The mcl and orthomclMclToGroups modules generated the groups.txt-file. We then used a custom Python script to extract only unique genes that were represented (1) in all species and (2) owing to their critical phylogenetic distance and poorly sequenced genomes, in all species except D. rerio and/or T. rubripes and/or G. aculeatus. We then extracted all corresponding nucleotide sequences from the CDS FASTA file by applying another custom Python script. GUIDANCE2 was then used to remove unreliably aligned positions in the multiple sequence alignments before testing for positively selected sites with PAML 4 (Yang, 2007). We conducted the analysis twice (e.g. the GUIDANCE2 filter and subsequent PAML analysis) using different aligners: (1) Muscle version 3.8.31 (Edgar, 2004) and (2) Mafft version X7.721 (Katoh and Standley, 2013). Fig. S15 (https://doi.org/10.6084/m9.figshare.13344806.v1) lists all guide trees that were used to fulfill branch-site tests for positive selection (M2-BS-H1 versus M2-BS-H0) in the K. marmoratus lineage using PAML. We inferred positive selection in K. marmoratus when the alpha-value for the likelihood ratio test was <0.01.

Treatment effects on collagen staining were tested using two-way ANOVA, with time and strain specified as main effects, including an interaction term, followed by a post hoc Holm–Šidák test. Treatment effects on desmosome length were tested using two-way ANOVA, with time and strain specified as main effects, including an interaction term. Differential expression analysis was performed in Limma v3.34.9 (Ritchie et al., 2015), where strain (HON, FW) and time (0, 1, 6, 24, 72, 168 h) were specified as fixed main effects, and strain×time as an interaction term, including false discovery rate adjustment of P-values.

We examined three sets of transcriptional responses: (1) transcripts with a significant temporal effect that was similar between HON and FW fish (shared air exposure response genes), (2) transcripts with a significant time×strain interaction (genes that differed in their air exposure response between HON and FW fish), and (3) transcripts that were not responsive to air exposure, but were constitutively differentially expressed between HON and FW fish (strain genes). Shared air exposure response genes were grouped by temporal patterns of co-expression, which were defined by the time post-air exposure at which expression was first significantly different from immersion controls. This grouping was performed separately for upregulated and downregulated (compared with immersion control) genes. For functional inference, clusters of co-expressed genes were examined for gene ontology (GO) and KEGG pathway enrichment using DAVID v6.8 (Huang et al., 2009) (tables S1–S11, available at https://doi.org/10.6084/m9.figshare.13344803.v1) using Uniprot IDs from human orthologs of K. marmoratus genes, and for canonical pathway and BioFunction enrichment using Ingenuity Pathway Analysis (Krämer et al., 2014).

We found remarkably rapid molecular changes in the skin within 1 h of air exposure. Genes associated with skin remodeling, including those involved in skin morphogenesis, cell junctions and extracellular matrix, were quickly upregulated. Histological images of the skin verified that after 7 days in air, junctional complexes and collagen content were altered. We also found a rapid upregulation of aquaporin genes and those genes involved in paracellular permeability, underlying the importance of skin osmoregulation. Large ionoytes with elaborate crypts were observed in the skin of air-exposed fish that were uniquely associated with basolateral capillaries. We also showed that genes linked to angiogenesis were quickly upregulated, consistent with an extensive capillary network close to the skin surface. Finally, our analyses of gene loci evolving by positive selection in K. marmoratus identified a number of genes involved in skin morphology, morphogenesis and angiogenesis, reinforcing the critical importance of the skin for transitioning effectively between water and air.

Air exposure induced rapid but mostly transient changes in transcription of genes associated with skin remodeling (Fig. 1). Transcription factors were enriched among the genes upregulated by 1 h in air (table S1, https://doi.org/10.6084/m9.figshare.13344803.v1), including those known to regulate skin morphogenesis (Fuchs and Raghavan, 2002), such as transcription factors P65 (NFkB or RELA), retinoic acid receptor, Krueppel-like factor 4, and transcription factors AP-1 (JUN, FOS) (fig. S2, https://doi.org/10.6084/m9.figshare.13344806.v1). Amphibian skin remodeling during metamorphosis is thyroid hormone-dependent (Brown and Cai, 2007). We observed parallel upregulation of THRA and RARG within 1 h of air exposure (fig. S2, https://doi.org/10.6084/m9.figshare.13344806.v1), which is consistent with a convergent role of thyroid hormone signaling in skin remodeling during terrestrial transitions between amphibians and amphibious fish (Brown and Cai, 2007; Inui and Miwa, 1985). This rapid and transient burst in gene expression that we observed is likely important for initiating skin remodelling, which is important for resilience to the aerial environment.

At cell junctions, multiple structures serve to attach epithelial cells together, in sheets, and to the basement membrane and extracellular matrix (ECM) (Fuchs and Raghavan, 2002), and regulate the exchange of small molecules and the passage of water. We found a rapid upregulation of genes involved in desmosomes and adherens junctions that provide strength to tissues (Fig. 1A; table S2, https://doi.org/10.6084/m9.figshare.13344803.v1). Focal adhesions are integrin-containing complexes that attach epidermal cells to the basement membrane and ECM, and here we also observed peak increases in the KEGG pathway focal adhesion 6 h after air exposure, including transcripts for integrin subunits/ligands (Fig. 1A; figs S3, S4, https://doi.org/10.6084/m9.figshare.13344806.v1). Furthermore, gap junctions (e.g. connexin-containing complexes) provide channels through which small molecules may pass between adjacent cells, and tight junction proteins (e.g. claudins, occludins) regulate paracellular permeability to water and small molecules. The molecular components of these structures, and the pathways that regulate them, were all transcriptionally regulated within hours of air exposure (Fig. 1A; figs S5, S6, https://doi.org/10.6084/m9.figshare.13344806.v1). Furthermore, TEM observations confirmed that differences in the configuration of specialized epidermal structures associated with the functional integrity of the epidermis occurred following air exposure (Fig. 2). Desmosomes serve as the intercellular glue that is required for epithelial integrity (Kowalczyk and Green, 2013), and desmosome proteins are evolving by directional selection in K. marmoratus (see below). We conclude that connections between keratinocytes that interface with the environment are, at least in part, regulated through desmosomes and that this is likely to be adaptively important for mechanical reinforcement of the skin during air exposure.

Structural constituents of the ECM were regulated during aerial acclimation. The GO terms ‘collagen’ and ‘extracellular matrix organization’ were enriched among transcripts downregulated between 1 and 7 days in air (Fig. 1A; table S6, https://doi.org/10.6084/m9.figshare.13344803.v1), including the downregulation of at least 23 collagen transcripts. This suggests that reduced collagen synthesis is an important component of skin remodeling during air exposure, which is consistent with our observed reduction in collagen staining in skin sections by 7 days in air (P<0.001; Fig. 3A). The thickness of the collagen layer or collagen content may affect exchange across the skin, but this is unknown, and the diffusion distance for gases, water and ions will be set, more importantly, by the closeness of the capillaries to the surface (Lillywhite and Maderson, 1988). Collagen also is critical for mechanical durability of the skin, and less collagen during terrestrial sojourns may be advantageous for locomotion (Szewciw and Barthelat, 2017) or other functions, but further studies are required to understand the role of collagen in K. marmoratus skin. TEM observations indicate that collagen is prominent in at least four regions: (1) the stratum compactum in the lower reaches of the dermis, (2) the scales, (3) the dermis–epidermis interface, and (4) the epidermal and sub-epidermal vasculature (Fig. 3B–E). Collagen in the sub-epidermal region of K. marmoratus skin exhibited an unusual structural arrangement relative to that reported for other fishes, as it appears in loosely ordered bundles (fig. S7, https://doi.org/10.6084/m9.figshare.13344806.v1) compared with the highly ordered and deep plywood-like lamellae localized to the upper reaches of the dermis in other fishes (e.g. Nadol et al., 1969).

In amniotes, the stratum corneum (SC) forms the outermost layer of the epidermis and functions as a barrier to protect from dehydration, crucial for the terrestrial radiation of the amniotes (Alibardi, 2016). Though fishes lack many of the key proteins that constitute the mature SC in amniotes (Strasser et al., 2014), our transcriptional data provide evidence that the molecular machinery that supports formation of the SC in amniotes is activated during air exposure. The transglutaminases (e.g. TGM1 and TGM5) that govern the process of protein cross-linking that produces the SC have been found in fish (Rodriguez Cruz et al., 2017). In K. marmoratus skin, TGM1 and TGM5 were significantly upregulated within 1 h of air exposure (fig. S8, https://doi.org/10.6084/m9.figshare.13344806.v1). The other key proteins that regulate this pathway were also upregulated after 1 and 6 h in air (e.g. EVPL, PPL, KAZN; Fig. 8). Furthermore, the sphingolipid metabolism pathway, which produces lipids integral to barrier function of the SC in amniotes (Lillywhite, 2006), is upregulated between 1 and 7 days after air exposure (Fig. 1A; table S3, https://doi.org/10.6084/m9.figshare.13344803.v1). TEM images provide no evidence for the formation of an SC (Figs 3 and 4), so the physiological manifestations of these SC-associated molecular responses merit further inquiry.

Kryptolebias marmoratus are able to effectively counteract desiccation while out of water; after several days in air they increase water influx across the skin, whereas water efflux does not change (Heffell et al., 2018), and fish maintain whole-body water levels (Leblanc et al., 2010). Enhanced water uptake may involve aquaporins or changes to paracellular tight junctions (Cerda and Finn, 2010; Kolosov et al., 2013), and upon air exposure we observed rapid upregulation of aquaporin 3, as well as claudin 10 and occludin tight junction transcripts (figs S6, S9, https://doi.org/10.6084/m9.figshare.13344806.v1).

In air-exposed fish, TEM revealed two types of sub-apical tight junctions in the epidermis. The most common were deep tight junctions (∼200 nm), which link most cells in the epidermis (fig. S10A, https://doi.org/10.6084/m9.figshare.13344806.v1). A second type of tight junction observed was a shallow (‘leaky’) tight junction, which linked ionocytes and adjacent interdigitating cells in an arrangement that was similar to the relationship between seawater fish gill ionocytes and interdigitating accessory cells. Shallow tight junctions directly interfaced with the environment within ionocyte apical crypts (fig. S10B,C, https://doi.org/10.6084/m9.figshare.13344806.v1); they were frequently observed in TEM images from HON fish acclimated to air, but none were observed in FW fish, despite the unusual presence of apical crypts in these animals (not a feature of FW gill architecture). In the seawater fish gill epithelium, ionocyte apical crypts establish a microenvironment within which paracellular movement of Na⁺ across leaky junctions can occur, down an electrochemical gradient established by transcellular Cl⁻ secretion through the Cl⁻ channel CFTR (Marshall and Grosell, 2005). Therefore, the enigmatic presence of deep apical crypts without ‘leaky’ junctions in ionocytes of air-exposed FW fish most likely reflects the need to create a microenvironment for solute transport, not for Na⁺ secretion but possibly for NH₃ trapping. In line with these observations, HON and FW animals differed in their air-induced expression of key components of fish ionocytes, including Na⁺/K⁺-ATPase, CFTR and claudins. Claudin-10 proteins are proposed to selectively facilitate Na⁺ movement (Bui and Kelly, 2014; Marshall et al., 2018) and here we observed much higher mRNA levels in HON fish (fig. S6, https://doi.org/10.6084/m9.figshare.13344806.v1). Therefore, our observations suggest that Na⁺ regulation through ‘leaky’ tight junctions that reside within apical crypts of HON epithelial ionocytes is important for aerial acclimation in this strain.

Ammonia-transporting Rhesus (Rh) proteins are localized to the apical crypt of epidermal ionocytes in K. marmoratus (Livingston et al., 2018), and are thought to facilitate cutaneous ammonia excretion on land (Litwiller et al., 2006). The genes coding for three Rh isoforms (Rh types A, B and C) and V-type H+-ATPase subunits were upregulated after 6 h of air exposure (fig. S11, https://doi.org/10.6084/m9.figshare.13344806.v1), consistent with the acid-trapping model proposed for branchial ammonia excretion (Wright and Wood, 2009). Also, glutamine synthetase (GLUL), which detoxifies ammonia (Hakvoort et al., 2017), is significantly increased at 1 h following air exposure (fig. S12, https://doi.org/10.6084/m9.figshare.13344806.v1). Interestingly, gene expression of the Na+/H+ exchanger (NHE2), also involved in ammonia excretion (Cooper et al., 2013), was 50-fold higher in HON fish relative to FW fish (fig. S13, https://doi.org/10.6084/m9.figshare.13344806.v1). This difference may reflect different strategies for ion transport, and possibly ammonia excretion, between fish acclimated/adapted to different salinities.

Morphological hallmarks of cutaneous respiration in vertebrates are augmented vascularization of the skin and the presence of blood vessels in the upper reaches of the epidermis, which reduces the blood to water diffusion distance and enhances epithelial gas transport in air (Noble, 1925). Our TEM images show that by 7 days in air, the diffusion distance from air to blood in epithelial capillaries is as short as 0.23 μm (Fig. 4C); to our knowledge, this is the shortest cutaneous diffusion distance yet measured in a fish. By comparison, diffusion distance across the skin of other amphibious fishes, as well as the gill epithelium of fishes in general, is at least one order of magnitude greater (e.g. Hughes et al., 1973). TEM observations also revealed the presence of blood vessels at the base of ionocytes (Fig. 4B), establishing a blood-to-water bridge, across which solutes could be efficiently transported.

Cutaneous angiogenesis and key angiogenesis-regulating genes were elevated within 1 day of air exposure in K. marmoratus, as reported in our companion study (Blanchard et al., 2019). Here, we show that genes involved in the KEGG pathway ‘VEGF signaling’ were enriched among the genes significantly upregulated by 6 h in air (table S2, https://doi.org/10.6084/m9.figshare.13344803.v1), and ‘angiogenesis’ and ‘vasculogenesis’ are BioFunctions that were significantly enriched (P=3.69×10^–8 and 4.65×10^–9, respectively, IPA analysis) among the genes upregulated by 24 h in air (Fig. 1A). These include transcription factors XBP1, NR4A1 and SRF (fig. S14, https://doi.org/10.6084/m9.figshare.13344806.v1), which regulate VEGF-dependent angiogenesis (e.g. Chai et al., 2004). Furthermore, bradykinin receptor promotes angiogenesis through nitric oxide signaling (Parenti et al., 2001), and we found parallel upregulation of BDKRB2 transcripts and NOS1 by 6 h in air (fig. S14, https://doi.org/10.6084/m9.figshare.13344806.v1). We conclude that angiogenesis, coupled with blood vessels in very close proximity to the apical surface and at the base of ionocytes, is important for shifting respiratory and ionoregulatory functions from the gill to the skin during aerial acclimation.

Considering 7-day survival data, the HON fish were more resilient to air exposure than FW fish (fig. S15, https://doi.org/10.6084/m9.figshare.13344806.v1); this could be due to genetic differences between strains, or to differences in acclimation environment. Though the FW strain inhabits a very different native habitat (freshwater) than HON strain animals (brackish/marine), likely resulting in evolved differences in physiology, the experiments reported here were not designed to distinguish genetic (heritable) from acclimation differences. Variable resilience could be related, in part, to skin function. Among the skin genes that were differentially expressed between strains during aerial acclimation (strain×time interaction, P<0.05), we observed three clusters of co-expressed genes (Fig. 1B). Set 1 is significantly enriched for the GO term ‘cell cycle’ (P=8.6×10^–56) (and other related GO terms: table S7, https://doi.org/10.6084/m9.figshare.13344803.v1): these genes show a transient spike in expression at 1 h for HON fish, but gradual downregulation and then recovery in FW fish. Set 2 is enriched for the GO terms ‘proteasome’ (P=5.3×10^–18) and ‘lysosome’ (P=9.4×10^–4) (table S8, https://doi.org/10.6084/m9.figshare.13344803.v1); these genes show gradual upregulation only in FW fish, with little temporal response in HON fish. Set 3 is enriched for the GO terms ‘translational initiation’ (P=8.7×10^–6), ‘cell junction’ (P=7.7×10^–5) and ‘sodium transport’ (P=3.1×10^–2) (table S9, https://doi.org/10.6084/m9.figshare.13344803.v1); these genes show gradual downregulation only in FW fish, with little temporal response in HON fish. Immediate differences in cell cycle regulatory response may indicate differences between HON and FW animals in their ability to sense and respond to cellular stress incurred upon air exposure (set 1). Cellular damage may have accumulated in the less-resilient FW strain, such that proteolytic activity is upregulated primarily in that strain (set 2). HON and FW animals differ in their regulation of structural and functional components of the osmoregulatory apparatus in ionocytes (set 3); indeed, in addition to the sodium transport genes (including Na⁺/K⁺-ATPase, Na⁺-bicarbonate cotransporter, Na⁺/H⁺ exchanger), a number of other key gill osmoregulation genes are differentially expressed between strains, including V-type H⁺ ATPase, CFTR, carbonic anhydrase, claudins and aquaporin (fig. S16, https://doi.org/10.6084/m9.figshare.13344806.v1). We conclude that several differences in the regulation of skin pathways from acclimation and/or adaptation to freshwater may constrain the ability of K. marmoratus to transition between water and air. Considering that Devonian tetrapods were likely derived from marine or brackish lineages (Bray, 1985; Goedert et al., 2018; Niedźwiedzki et al., 2010), and that brackish and terrestrial environments share challenges posed by dehydration, perhaps salty aquatic environments favor osmoregulatory physiologies that are more easily co-opted for terrestriality than freshwater environments. Our data are consistent with this hypothesis and suggest that physiological differences in ionoregulation are mechanistically associated with terrestrial acclimation abilities.

Our initial test for positive selection in K. marmoratus protein sequences included 1073 single-copy orthologous gene loci in all nine fish species included (table S12, https://doi.org/10.6084/m9.figshare.13344803.v1; fig. S1, https://doi.org/10.6084/m9.figshare.13344806.v1), and revealed 24 genes showing positive selection in the K. marmoratus lineage only (table S13, https://doi.org/10.6084/m9.figshare.13344803.v1; fig. S1, https://doi.org/10.6084/m9.figshare.13344806.v1). An extended analysis expanded the query set to include orthologs present in all species except one of either D. rerio, T. rubripes or G. aculeatus; this included 461 further genes for analysis (table S14, https://doi.org/10.6084/m9.figshare.13344803.v1), and revealed 16 more gene loci evolving by positive selection in the K. marmoratus lineage (table S15, https://doi.org/10.6084/m9.figshare.13344803.v1; fig. S1, https://doi.org/10.6084/m9.figshare.13344806.v1). Because of the constraints for confidently identifying single-copy orthologs across multiple taxa, this analysis queried only ∼4% of the genes in the K. marmoratus genome (table S16, https://doi.org/10.6084/m9.figshare.13344803.v1).

We detected positive selection in genes involved in skin morphology and morphogenesis, including two key genes that contribute to desmosome structure (desmoglein and plakophilin) and in transcription factor Kruppel-like factor 9 (KLF9), which regulates skin morphogenesis (Spörl et al., 2012). Desmosomes are an evolutionary innovation of vertebrates (Green et al., 2020) and share a conserved role in maintaining epidermal integrity to mechanical stress including in amphibians (Bharathan and Dickinson, 2019). Desmocollin and desmoglein genes have expanded in the tetrapods (Green et al., 2020) and desmosome genes are subject to selection during re-invasion of aquatic environments in mammals (Chikina et al., 2016), consistent with a role for desmosomes in terrestrial-aquatic adaptations. Positive selection was also detected in heart- and neural crest derivatives-expressed protein 1 (HAND1), which is required for angiogenesis, and regulates the expression of many angiogenic pathways including VEGF, angiopoietin and ephrin signaling (Morikawa and Cserjesi, 2004). Pathways related to skin remodeling, desmosome function and angiogenesis are also among those showing patterns of differential regulation during terrestrial acclimation (Fig. 1A). We conclude that some of the proteins evolving by positive selection within the K. marmoratus lineage are adaptive for supporting skin plasticity, which enables the amphibious lifestyle of this species.

Amphibious abilities have evolved 87 times independently among the extant fishes (Damsgaard et al., 2020), many species of amphibians transition between aquatic and terrestrial forms during metamorphosis, and three lineages of terrestrial mammals have transitioned back to fully aquatic forms (cetaceans, pinnipeds and sirenians). This diversity provides ample opportunity to study the conserved and divergent mechanisms that support transitions between aquatic and terrestrial physiologies that are so important in the history of vertebrate diversification. So far, few studies have detailed the mechanisms that underlie skin plasticity and evolution during air exposure. Our integrated molecular, structural and evolutionary analyses revealed the importance of rapid regulation of cell–cell adhesion structures, molecular pathways that govern skin morphogenesis, and evolution of desmosome proteins to mechanically reinforce and buttress barrier functions of this tissue. The challenge of supporting respiratory and ion exchange functions is associated with ionocyte structures and rapid regulation of ionoregulatory pathways, angiogenesis, evolution of angiogenesis-regulating proteins, and an unusual blood vessel architecture that links very short air–blood diffusion distance with delivery to energetically demanding ionocytes. Studies in marine mammals and amphibious gobies provide evidence of parallel importance of skin morphogenesis (Chikina et al., 2016) and ion transport mechanisms (You et al., 2014), respectively, in aquatic–terrestrial transitions. Overall, our data reveal that when K. marmoratus transitions from water to air, the skin quickly alters its barrier and exchange functions to delicately balance these dual purposes on land.

We thank Dr Lisa Johnson for bioinformatics advice and training, Jennifer Roach for assistance with RNA isolation and RNA-seq library preparation, Norbert Grundmann for assistance in using the PALMA cluster Muenster, Chun Chih Chen for assistance with morphometric work, Drs Nick Bernier, Beren Robinson and Andy Turko for advice on experimental design, and Caitlin Brookbanks for clerical assistance.