# Mechanical positional information guides the self-organized development of a polygonal network of creases in the skin of mammalian noses

Paule Dagenais, Ebrahim Jahanbakhsh, Aurélien Capitan, Hélène Jammes, Karine Reynaud, Camino de Juan Romero, Victor Borrell, Michel C Milinkovitch

<figure><img src="image_1.png" /></figure>

## To cite this version:

Paule Dagenais, Ebrahim Jahanbakhsh, Aurélien Capitan, Hélène Jammes, Karine Reynaud, et al. "Mechanical positional information guides the self-organized development of a polygonal network of creases in the skin of mammalian noses." Current Biology, 2024, 34, pp.1-16.e1-e4.
10.1016/j.cub.2024.09.055. hal-04765961

HAL Id: hal-04765961

https://hal.inrae.fr/hal-04765961v1

Submitted on 4 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Current Biology

## Mechanical positional information guides the self-organized development of a polygonal network of creases in the skin of mammalian noses

### Graphical abstract

<figure><img src="image_3.png" /></figure>

### Authors

Paule Dagenais,
Ebrahim Jahanbakhsh,
Aurélien Capitan, ...,
Camino De Juan Romero,
Victor Borrell, Michel.C. Milinkovitch

### Correspondence

michel.milinkovitch@unige.ch

### In brief

Dagenais et al. demonstrate that the characteristic polygonal folding patterns on the noses of many mammalian species, such as dogs, ferrets, and cows, emerge via a peculiar self-organizational process involving mechanical positional information encoded by the spatial distribution of stiffness in adjacent tissues.

### Highlights

* Rhinarial polygonal domains are not placode-derived but mechanically self-organized

* The rhinarial skin pattern is mechanically coupled with underlying blood vessels

* Epidermal domes rise against stiffer blood vessels—like arches against pillars

* Rhinarial folds are as variable among cow clones as among unrelated individuals
<figure><img src="image_4.png" /></figure>

Article

# Mechanical positional information guides the self-organized development of a polygonal network of creases in the skin of mammalian noses

Paule Dagenais,<sup>1,2</sup> Ebrahim Jahanbakhsh,<sup>1,2</sup> Aurélien Capitan,<sup>3</sup> Hélène Jammes,<sup>4</sup> Karine Reynaud,<sup>5,6</sup>

Camino De Juan Romero,<sup>7</sup> Victor Borrell,<sup>7</sup> and Michel.C. Milinkovitch<sup>1,2,8,9,*</sup>

<sup>1</sup>Laboratory of Artificial & Natural Evolution (LANE), Department of Genetics & Evolution, University of Geneva, 1211 Geneva, Switzerland

<sup>2</sup>SIB Swiss Institute of Bioinformatics, 1211 Geneva, Switzerland

<sup>3</sup>Université Paris-Saclay, INRAE, AgroParisTech, GABI, 78350 Jouy-en-Josas, France

<sup>4</sup>BREED INRAE, UVSQ, Université Paris-Saclay, 78350 Jouy-en-Josas, France

<sup>5</sup>CNRS, IFCE, INRAE, Université de Tours, PRC, 37380 Nouzilly, France

<sup>6</sup>École Nationale Vétérinaire d'Alfort, EnvA, 94700 Maisons-Alfort, France

<sup>7</sup>Instituto de Neurociencias, CSIC-UMH, 03540 San Juan de Alicante, Spain

<sup>8</sup>X (formerly Twitter): @LANEVOL

<sup>9</sup>Lead contact

*Correspondence: michel.milinkovitch@unige.ch

https://doi.org/10.1016/j.cub.2024.09.055

## SUMMARY

The glabrous skin of the rhinarium (naked nose) of many mammalian species exhibits a polygonal pattern of grooves that retain physiological fluid, thereby keeping their nose wet and, among other effects, facilitating the collection of chemosensory molecules. Here, we perform volumetric imaging of whole-mount rhinaria from sequences of embryonic and juvenile cows, dogs, and ferrets. We demonstrate that rhinarial polygonal domains are not placode-derived skin appendages but arise through a self-organized mechanical process consisting of the constrained growth and buckling of the epidermal basal layer, followed by the formation of sharp epidermal creases exactly facing an underlying network of stiff blood vessels. Our numerical simulations show that the mechanical stress generated by excessive epidermal growth concentrates at the positions of vessels that form rigid base points, causing the epidermal layers to move outward and shape domes—akin to arches rising against stiff pillars. Remarkably, this gives rise to a larger length scale (the distance between the vessels) in the surface folding pattern than would otherwise occur in the absence of vessels. These results hint at a concept of “mechanical positional information” by which material properties of anatomical elements can impose local constraints on an otherwise globally self-organized mechanical pattern. In addition, our analyses of the rhinarial patterns in cow clones highlight a substantial level of stochasticity in the pre-pattern of vessels, while our numerical simulations also recapitulate the disruption of the folding pattern in cows affected by a hereditary disorder that causes hyperextensibility of the skin.

## INTRODUCTION

Many morphological patterns in biology emerge from interactions of signaling molecules (morphogens) that act as regulators of cell dynamics and whose distributions in time and space can be described by Turing's reaction-diffusion self-organizational process.<sup>1,2</sup> For example, many skin appendages (avian feathers, mammalian hairs, and most reptilian scales) develop from placodes, i.e., anatomical and biochemical signaling centers organized as reaction-diffusion-taxis systems.<sup>5-6</sup> Nevertheless, this framework alone is insufficient to account for the full diversity of patterns found in living organisms—in particular, folded morphologies that arise as mechanical instabilities generated by differential growth between adherent tissues.<sup>7</sup> Hence, tissue patterning dominated by mechanics offers an essential and complementary viewpoint to understand morphogenesis

beyond the classic paradigms of chemical positional information<sup>8</sup> and Turing's patterns.<sup>1</sup> In fact, it is most probable that many biological systems acquire their complex shape through an interaction of the latter two processes with mechanical forces.<sup>9-11</sup>

Constrained growth in multilayered soft matter is a common means through which elaborate patterns can emerge spontaneously without a prior stencil of molecular positional information, resulting in typical networks of sharp creases or smooth folds. Growth-induced folding patterns have been studied analytically<sup>12,13</sup> and in numerical models.<sup>7,14,15</sup> The purely mechanical aspects of constrained growth can be experimentally investigated with gel analogs.<sup>15-22</sup> Indeed, gel polymers swell when immersed in specific solvents, causing mechanical instabilities to develop between the static core and the growing outer layer as the liquid penetrates from the surface. For example, Figure S1
illustrates how polygonal mosaics can form spontaneously in a stratified tissue (here a bilayer hydrogel) due to constrained growth. Global regimes of surface folds are predictable and can be classified into phase diagrams based on variables distinguishing the two layers, such as strain, compressive stress, stiffness, shear modulus, thickness, and curvature.⁷,¹⁵,²³⁻²⁵ Hence, layered soft materials can self-organize into complex folded morphologies through a purely mechanical process, with the potential to form convoluted (“brainy”) shapes, polygonal tessellations, or wavy ridges.

Multiple examples of folded architectures in biology have been successfully described in terms of growth-induced mechanical patterning, from (sub-)micro-scale geometries such as the membrane of white blood cells,²⁶ the surface nano-ridges in plant petal epidermal cells,²⁷,²⁸ and the convoluted wall of respiratory airways,²⁹ to macroscopic structures like fingerprints³⁰ (but see Glover et al.³¹), the inner mucosa of tubular organs such as the stomach and the gut,²⁴,³² the surface of tumors,³³,³⁴ the cerebral cortex,³⁵,³⁶ the buckled geometry of grass blades and leaves,³⁷ and the wrinkled integument of various animals and plants.²⁵ Note that these folded architectures often correspond to somewhat smooth undulations, whereas examples of mechanically patterned polygonal networks with sharp creases—as during the development of crocodile head scales (E.J. and M.C.M., unpublished data, and Milinkovitch et al.³⁸)—seem substantially less frequent in living systems. Natural folded shapes provide great inspiration for biomimetic engineering, sparking the design of artificial organs and sensors in which multi-scale patterning occurs spontaneously under controlled mechanical conditions without having to resort to traditional and tedious techniques such as photolithography for template design.²⁵,³⁹

A beautiful example of folding patterns is found in the glabrous skin of the rhinarium (naked nose) of many mammalian species, including dogs, cows, and ferrets, in which polygonal domes are surrounded by sharp epidermal creases (Figure 1). These patterns were first cataloged in dozens of species by Hill,⁴⁰ who named them rhinoglyphics by analogy with dermoglyphics (fingerprints). The rich inventory of rhinoglyphics, from the oblique elongated ridges of cows to the small polygonal domes of dogs (Figure 1), might reflect a diversity of functions.⁴¹,⁴² The innervation of the epidermal domes suggests that they form individual sensory units.⁴² One aspect common to all species is that the rhinarial skin grooves keep the nose wet by collecting and retaining physiological fluid that is either spread on the nose by the tongue,⁴³ or released by gland ducts—opening in the nostrils, as in dogs⁴⁴ or directly on the surface of the epidermal domes, as in cows.⁴⁵ A moist surface facilitates the collection of airborne or substrate-bound molecules, which can subsequently travel into the middle cleft of the nose (philtrum) and then into the naso-palatine ducts in the anterior palate, leading to the vomeronasal organs for chemosensation.⁴⁶⁻⁵⁰ In addition, mammals tend to keep their rhinaria at specific temperatures,⁴¹,⁴⁵,⁵¹ and the fluid collected in skin creases might contribute to evaporation-based thermoregulation, similarly to the retention and subsequent evaporation of water in the network of crevices that adorns the skin of African elephants.⁵² In rodents, thermoregulation via the rhinarium participates to selective brain cooling,⁴³ while in carnivores, a cold rhinarium is thought to favor heat and infrared radiation sensing for more efficient preying.⁴¹,⁴²,⁵³

Other studies have suggested a tactile role for rhinoglyphics, similar to fingerprints, which could be facilitated by a warm rhinarium as in herbivores.⁴⁰⁻⁴²,⁵⁴⁻⁵⁶

Here, using opportunistic sampling of embryos, we study the ontogenic development of the rhinarial skin in cows, dogs, and ferrets. To this end, we use light-sheet microscopy to image whole-mount rhinaria from sequences of embryonic and juvenile stages, from which we extract the corresponding 3D multilayered skin geometries and quantify the spatial distribution of proliferating cells in the dermis and epidermis. These analyses indicate that the polygonal networks of creases of these three mammalian species arise during embryogenesis as a self-organized folding process with the particularity of being mechanically coupled with an underlying network of stiff blood vessels. In parallel, we successfully recapitulate the polygonal patterning of rhinarial skin by performing numerical simulations of growth in a 3D domain combining the multilayered skin and underlying vasculature—first in a simplified flat-slab geometry with biomimetic mechanical properties, then in a realistic model of a dog embryonic rhinarium reconstructed from light-sheet data. Moreover, our growth simulations allow us to explore the mechanical parameter space (growth rate and stiffness of individual skin layers) and pinpoint transitions in folding regimes corresponding to biologically relevant versus unnatural surface patterns.

Through these complementary approaches, we identify the following common mechanism that governs the morphogenesis of rhinoglyphics across different species. At the onset of ontogenic patterning, excessive growth of the epidermis basal layer causes its buckling. However, the spatial distribution of the resulting compressive field is mechanically constrained by the presence of a sub-epidermal template of stiffer blood vessels that form rigid base points such that, at later stages, the sharp creases that form polygonal domes in the epidermis develop along the underlying network of vessels. These results prompt us to propose a mechanical equivalent to chemical positional information.⁸,⁵⁷ Indeed, the development of rhinoglyphics in cows, dogs, and ferrets indicates that mechanical morphogenetic patterns based on growth-induced self-organization can be strongly constrained and guided by a pre-existing template of mechanical positional information. Other biological systems have been shown to rely on mechanical interactions with the surrounding anatomical features to attain their shape, for example, the branching morphogenesis in mouse lungs, where epithelial growth is constrained by smooth muscles.⁵⁸ However, the patterning process uncovered here is strikingly different in the sense that the network of folds can be predicted exactly from a pre-existent anatomical pattern formed by the vessels, i.e., the stochastic component of motif positions (that usually characterizes self-organized folding patterns) is erased by the effect of mechanical coupling with a pre-existing pattern of vessels. In other words, the presence of underlying stiff vessels transforms rhinarial skin folding into a deterministic process, akin to the formation of hairs, feathers, or scales that follow a pre-existent map of local chemical cues (the signaling placodes).⁴,⁶

## RESULTS

The three species studied here (cow, dog, and ferret) exhibit distinct adult rhinarial skin phenotypes (Figures 1A-1C). To
Current Biology

Article

<figure><img src="image_5.png" /></figure>

CellPress

OPEN ACCESS

<figure><img src="image_6.png" /><figcaption>**Figure 1.** Patterns of polygonal domes in the folded rhinarial skin of three mammalian species</figcaption></figure>

(A) Cows, (B) dogs, and (C) ferrets. Left column: adult phenotypes; insets show magnified portions of the rhinarium with one polygon marked in blue (edges) and red (nodes). White arrowheads: shallower sulci within the main epidermal domes. White arrow: gland duct opening at the surface of the cow rhinarium. Central column: 3D surface reconstructions (right side) of adult individuals with maximum curvature (left side) highlighting the epidermal creases surrounding the polygonal domes. Right column: statistical distributions of polygonal order and aspect ratio of epidermal domes, averaged over all samples for each species (13 cows, 2 dogs, and 2 ferrets), with corresponding color-coded spatial distributions on rhinarial skin surfaces for one sample of each species. Color bars for the polygon order are the same for all three species.

See also Figures S1–S3 and S6.

quantify their differences, we performed statistical analyses of
the corresponding lattices of polygonal epidermal domes after
their segmentation on 3D nose geometries using surface curva-
ture (Figures 1, central column of panels, and S2A). The lattices
of the three species exhibit different statistical attributes (Fig-
ure 1, right column of panels). For example, cows exhibit elon-
gated units (Figure 1A) diagonally oriented with respect to the
axis of the philtrum, whereas dogs and ferrets display polygons
with lower aspect ratios (ARs) and no preferred orientation
(Figures 1B and 1C). The rhinaria of cow, dog, and ferret can
also be distinguished using the mean relative area of the poly-
gons (i.e., relative to the size of the nose), which is much smaller
in dogs than in the other two species. On the other hand, the

mode of the polygonal order distributions is identical in all three
species, indicating that the most common shape for epidermal
domes is the pentagon. Interestingly, this contrasts from the
polygonal lattices that form spontaneously in hydrogel models
with a majority of hexagons¹⁶ (see also Figure S1). Moreover, nu-
merical simulations of growth instabilities in highly compressible
solids (i.e., with low Poisson's ratio) also result in mostly hexag-
onal patterns. ¹⁴ Note that secondary, shallower sulci are
observed in all three species, creating incomplete subdivisons
or dimples within the main polygons.

Histological sections of the rhinaria at different developmental
stages in cow, dog, and ferret embryos (Figure S3) show the un-
dulated basal layer of the epidermis and a primary vascular
<figure><img src="image_7.png" /><figcaption>**Figure 2. Timeline of hierarchical folding in rhinarial skin during embryogenesis**</figcaption></figure>

First row: schematic timeline representation of three successive key morphological events; blood vessels (white circles), the basal layer (dark gray), the dermis (orange), the epidermis (light gray), and secretory glands (purple) are shown. Second to fourth rows: histological sections (hematoxylin and eosin staining) perpendicular to the skin surface at selected embryonic stages (corresponding to the four stages of the timeline) in cows, dogs, and ferrets, respectively. Embryonic days (E) are indicated.

(A) Stage 0: a pre-pattern of blood vessels is located beneath the epidermis, slightly protruding against the basal layer.

(B) Stage 1: folding of the basal layer, forming outward (cows) or inward (dogs and ferrets) cups between the vessels.

(C) Stage 2: folding of the epidermis, forming surface creases facing the vessels.

(D) Stage 3: deepening of the creases, with secretory gland ducts emerging within the epidermal domes (cows), and remodeling of the basal layer (dogs and ferrets). Black arrowheads: blood vessels; black arrows: epidermal creases; black outline in cow section: glandular tissue; white dotted lines: basal layer of epidermis (main cups between the vessels); white solid lines: secondary folds of basal layer over the vessels. All scale bars are 0.1 mm.
See also Figure S3.

plexus⁵⁹ consisting of a polygonal network of sub-epidermal blood vessels. At embryonic day 150 (E150), sections in the cow rhinarium show the presence of dermal glands whose ducts (Figures S3Ai and S3Aii) form rows at the surface of the polygonal domes (see inset of Figure 1A). Conversely, in dogs and ferrets, fluid is secreted within the nostrils by the lateral nasal glands,⁴⁴ and molecules trapped by the fluid in the network of epidermal creases are drained into the philtrum (blue arrowheads, Figures S3B and S3C), then into the naso-palatine ducts leading to the chemosensory vomeronasal organs (red arrowheads, Figure S3), as illustrated by the sketch in the inset of Figure S3B.

## Ontogeny of the rhinarial skin pattern

Using histological sections and volumetric reconstructions from light-sheet data of whole-mount rhinarium samples at selected embryonic stages, we identify the three main morphogenetic

events (hereafter called stages 1, 2, and 3) leading to the final rhinarial skin pattern observed in the three species (Figure 2). First, at stage 0, a plexus of large blood vessels is observed below the epidermis before the onset of the folding process (Figure 2, column A). The morphogenesis of rhinoglyphics is then initiated at stage 1 with undulations of the epidermal basal layer (Figure 2, column B), forming outward (in cows) or inward (in dogs and ferrets) cups between the vessels. The surface of the epidermis then folds at stage 2 (Figure 2, column C), forming outward bumps surrounded by skin creases tracing lines positioned exactly above the sub-epidermal blood vessels. Hence, the epidermal creases forming the rhinoglyphics observed after birth mirror the network of sub-epidermal embryonic blood vessels. During the final stage 3 (Figure 2, column D), the epidermal domes continue to grow, deepening the adjacent creases while the skin surface keratinizes. In cows, gland ducts also emerge at the skin surface at that stage.
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 3.** Volumetric imaging (light-sheet microscopy) of the sequential patterning of skin folds in the embryonic rhinaria of cows, dogs, and ferrets</figcaption></figure>

(A, D, and G) Stage 1: folding of the basal layer around the pre-pattern of vessels.

(B, E, and H) Stage 2: onset of epidermal folding, with surface creases mirroring the underlying vessels.

(C, F, and I) Stage 3: deepening of epidermal creases, with rows of gland ducts opening within each dome for cows, and remodeling of the basal layer for dogs and ferrets. The images are maximum intensity projections (MIPs) of whole rhinaria (A, D, and G and left panels of B, E, H, F, and I) or slabs close to the skin surface (right panels of B, E, and H, both panels of C and panel (ii) of I). Panels (i) of (F) and of (l) show surface normal shading. Nuclear staining signal (TOPRO3) is shown in gray; autofluorescence signal (keratinized epidermal layer) is shown in green for (E) (right), (F), and (l) (left); autofluorescence of red blood cells is shown in green for (C) (left) and (F) (panel (ii)). Red dotted lines: epidermal surface; red arrows: epidermal creases; green lines and green shaded areas: main folds of the epidermal basal layer between the vessels; yellow dotted lines: secondary folds of the basal layer above the vessels; white dotted lines: blood vessels; white arrows: gland ducts. Scale bars: 0.4 mm. Anatomical axes indicated in (A) (l-r, left-right; r-c, rostral-caudal; d-v, dorsal-ventral) are valid for all panels except where otherwise stated.

See also Figures S4 and S7 and Videos S1, S2, S3, S4, S5, S6, and S9.

At the same time, the basal layer is being remodeled from its initial folding template to more elaborate motifs specific to each species: in dogs and ferrets, the primary inward fold of the basal layer (white dotted lines in Figure 2, column C) expands and subdivides, while a secondary network of smaller inward folds (white solid lines in Figure 2, column D) appears over the vessels (black arrow-heads), right below the epidermal creases (black arrows). As embryos approach their full development, the sub-epidermal blood vessels become much smaller relative to the size of the rhinarium, such that the correlation between their distribution and the rhinarium skin folding pattern in juveniles and adults can easily go unnoticed. However, this correlation is indeed maintained in adults, supporting a thermoregulatory function of the rhinoglyphics: evaporation of the fluid in the skin creases might cool down the underlying blood circuit.

To better visualize the full 3D topology of individual skin layers at each key step of rhinoglyphics morphogenesis, we produced light-sheet volumetric microscopy imaging of rhinaria at selected embryonic stages in each species (Figure 3; Videos S1, S2, S3, S4, S5, and S6). Maximum intensity projections (MIPs) for the whole rhinaria at the first stage of rhinoglyphics development (Figures 3A, 3D, and 3G) demonstrate that (1) blood vessels (dashed white lines) form a template that presets the epidermal surface pattern emerging at stages 2 and 3, and (2) the epidermal basal layer (green lines and green shaded areas) forms a complementary motif by folding into cups between the vessels. In cows, rows of developing gland ducts are already visible within each epidermal unit at stage 1 (Figure 3A, white arrows), although they do not emerge at the skin surface until stage 3. At stage 2, the upper epidermis is still barely keratinized, and
<figure><img src="image_10.png" /><figcaption>**Figure 4.** 3D segmentation of rhinarial tissues from volumetric light-sheet microscopy images of dog rhinaria</figcaption></figure>

(A–C) Segmented dermis (pink), epidermis (beige), basal layer (dark gray), blood vessels (red), and nasal cartilage (purple) are shown at three embryonic stages: E38 (A), E40 (B), and E44 (C). Anatomical axes are indicated: left-right, rostral-caudal, and dorsal-ventral.

(A) Top: transversal slice of the labeled image close to the base of the nose. Bottom: 3D surface reconstruction showing the network of blood vessels (left) protruding against the epidermal basal layer (right) viewed from inside the nose.

(B) Top: transversal slice of a light-sheet image close to the base of the nose (TOPRO3 nuclear staining with overlaid segmentation of the basal layer in blue). Bottom: 3D surface reconstruction showing the onset of epidermal creases formation viewed from outside the nose (left) as well as inward cups in the basal layer viewed from inside the nose (right). Note how these two anatomical elements mirror one another. Dotted red lines highlight the corresponding positions of the blood vessels.

(C) Top: sagittal slice intersecting the nostril, with top and bottom regions of the rhinarial skin magnified in the left insets; the right insets (i) and (ii) show the 3D-segmented blood vessels and basal layer (viewed from inside the nose) for the corresponding regions. Bottom: deepening of skin creases surrounding the polygonal domes in the epidermal surface (left), highlighted by the color-coded maximum curvature of the surface (right).
See also Figure S5 and Videos S2, S3, S7, and S8.

MIPs of the whole rhinaria reveal the polygonal patterns of the deeper basal layer and blood vessels (left panels of Figures 3B, 3E, and 3H; Video S3). In contrast, MIPs of transversal slabs taken close to the nasal surface (right panels of Figures 3B, 3E, and 3H; Videos S1 and S2) intersect all skin layers and thus unveil their mirrored topologies: at the surface, nascent epidermal creases (red arrows) appear over the template of sub-epidermal vessels (dashed white lines). At stage 3 (Figures 3C, 3F, and 3I; Videos S4, S5, and S6), the upper epidermis further keratinizes and can therefore be imaged with autofluorescence, revealing surface creases surrounding the well-defined polygonal domes. At this point, the secretory glands in cows are fully developed and connected to the surface by ducts (white arrow in Figure 3C). Autofluorescence of red blood cells confirms that the sub-epidermal vessels are indeed blood vessels (dashed white lines in Figures 3C and 3F). Note that, at this stage, the basal layer of the epidermis in dogs and ferrets is remodeling, as the main cups are subdivided by shallower sulci and dimples. At the same time, a secondary network

of inward folds (dashed yellow lines) develops right below the epidermal creases, therefore also mirroring the lines of the vessels.

To facilitate the 3D visualization of the processes described above, we use volumetric imaging data of the developing dog rhinarium (e.g., Videos S3, S4, and S5) to segment the main layers of the skin (upper epidermis, basal epidermal layer, dermis, and blood vessels). This analysis (Figure 4; Videos S7 and S8) allows us to highlight subtle anatomical features that contribute to the patterning of the polygonal folds but are otherwise hard to illustrate in the 2D MIPs of Figure 3. At stage 0 (E38), a network of blood vessels has formed in the upper dermis, with the vessels slightly protruding against the basal layer of the epidermis (Figure 4A). At E40, the basal layer has formed inward folds between the vessels (stage 1). Concurrently, very shallow epidermal creases begin to form directly above the vessels (stage 2), i.e., around the emerging outward domes (Figure 4B, lower). As stage 3 (E44) is approached, the surface creases deepen, making the polygonal domes more defined (Figure 4C,
# Current Biology

## Article

<figure><img src="image_11.png" /></figure>

CellPress

OPEN ACCESS

lower, Videos S2 and S3), while the epidermal basal layer begins to remodel into a more complex hybrid regime (Figure 4C, upper), with some of the folds facing outward (inset i) and others facing inward (inset ii) in-between the blood vessels.

The two observed folding events—the initial folding of the epidermal basal layer over the embryonic sub-epidermal blood vessels and the subsequent creasing of the epidermal surface on top of these same vessels—likely correspond to mechanical self-organization dominated by homogeneous growth. Alternatively, it could be argued that folding and creasing occur because of successive local proliferation (i.e., at the positions of, or in-between, vessels) of the epidermal basal layer and then of the upper epidermis, respectively. In principle, these two alternative morphogenetic processes (homogeneous versus localized proliferation) can be distinguished with assays identifying the spatial distribution of proliferating cells. Hence, we performed light-sheet microscopy imaging of full rhinaria after whole-mount pH3 (phospho-histone H3) labeling in each species at the three key developmental stages identified above (Video S1). At each of these stages, the distributions of pH3-marked cells appear homogeneous within each skin layer and do not show points of focalized proliferation (Figures S4 and S5; Video S9), i.e., we do not observe any spatial correlation between the distribution of proliferating cells and the distribution of skin folds. Figure S5 shows the relative spatial density of proliferating cells, i.e., the number of pH3-positive cells divided by the total number of cells within a sphere centered on each pH3-positive cell. By varying the diameter of the sphere (Figure S5B), we observe that proliferation patterns only emerge at length scales larger than the characteristic scale of the folding patterns, except for proliferation foci at the bottom of the noses, which correspond to hair follicles. Therefore, the spatial distribution of cell proliferation is much smoother than what would be required for patterning due to local proliferation. Note, however, that this analysis is limited by the fact that processes other than cell proliferation—such as extracellular matrix (ECM) deposition, cell hypertrophy, or contractility, which were not measured here—might be modulating local tissue growth locally. On the other hand, the relative growth rates vary among the dermis, basal layer, and epidermis during the course of ontogeny. For example, at the stage of epidermal basal layer folding, just before the onset of epidermal surface patterning, we observe that the relative densities of proliferating cells are uniformly higher in the epidermal basal layer compared with these densities in the rest of the epidermis and in the dermis. More specifically, in a dog embryo at E38, the normalized densities of proliferating cells (numbers of proliferating cells per unit volume) are 1.0, 6.4, and 1.5 for the dermis, basal layer, and epidermis, respectively, while the normalized total cell densities (number of proliferative and non-proliferative cells per unit volume) are 1.0, 1.3, and 1.1 between these layers, yielding relative proportions of proliferating cells of 1.0, 4.8, and 1.4. At E40, the normalized total cell densities become 1.0, 1.6, and 1.6 (i.e., the dermis remains less dense than the epidermal layers), and the relative proportions of proliferating cells become 1.0, 4.0, and 1.3 (i.e., the relative growth remains higher in the epidermal basal layer than in the dermis and upper epidermis). Note that these measurements rely on the 3D segmentation of skin layers from volumetric

light-sheet data, which was done for one sample per stage (E32, E38, and E40), as shown in Figure S5.

Thus, the combined observations of homogeneous proliferation within different layers and substantial differences in growth rates between layers strongly support the hypothesis that the initial buckling of the basal layer and the subsequent creasing of the epidermal surface that generates the final rhinarial folding patterns both emerge from mechanical self-organization rather than being driven by local cell proliferation at the positions of (or in-between) vessels.

## Growth simulations in multilayered skin models

To further explore the mechanical coupling between the multilayered skin and the underlying vascular network, we perform numerical simulations of growth in an idealized flat slab of skin with a square grid of blood vessels at the dermo-epidermal interface. We adjust the relative amount of tangential growth within each tissue layer to values estimated from light-sheet imaging of pH3 labeling: higher growth in the epidermal basal layer, intermediate growth in the upper epidermis, and lower growth in the dermis. Due to the technical difficulty of evaluating the elastic (Young's) moduli of individual skin layers, we use the cellular densities obtained from light-sheet data as a simple proxy for stiffnesses—translating to a stiffer basal layer compared with the dermis and epidermis. This approximation is reasonable for epithelial cells with little ECM but possibly problematic in the case of connective tissues where a complex ECM has a high influence on elastic properties (e.g., tendons and ligaments are hypocellular but orders of magnitude stiffer than tissues like the liver with high cell density⁶⁰). Note also that direct measurements made by Glover et al.³¹ on embryonic human skin have revealed a higher stiffness for the upper epidermis compared with the basal layer. Therefore, we test here both situations, i.e., stiffer upper epidermis or stiffer epidermal basal layer (Figure 5). Importantly, we follow an effective modeling approach where we do not seek to replicate the exact absolute parameter values of the system. Instead, we are testing the combined effects of growth ratios and relative stiffnesses between the different skin layers, whereas it would be irrelevant to insist on implementing exact numerical values for each mechanical property.

First, we show that, given the right set of parameters, such simulations on the flat-slab model recapitulate both the initial buckling of the basal layer forming inward folds between the underlying vessels (yellow arrowheads in Figure 5A) and the non-intuitive positioning of the epidermal surface creases over a template of vessels (black arrows in Figure 5A). The relative stiffnesses and growth values of the modeled skin layers in Figure 5A are inspired by the cellular and proliferation densities obtained from light-sheet microscopy data (higher density of cells in the epidermal basal layer and higher growth in the epidermal layers, especially in the basal layer).

In the absence of vessels, the epidermis folds into a brainy pattern (Figure 5B) whose characteristic length scale (indicated by a black scale bar) is proportional to the thickness of the upper layer.¹²,⁶¹ In contrast, the template of vessels introduces a new length scale in the surface folding pattern (the distance between the vessels, indicated by a white scale bar in Figure 5B), which, remarkably, supplants the natural length scale of the system. An essential condition to obtain the surprising, but correct,
<figure><img src="image_12.png" /><figcaption>Figure 5. Numerical simulations of growth in a flat-slab multilayered skin model with square network of vessels</figcaption></figure>

$E_b$, $E_d$, $E_e$, and $E_v$ are the Young's moduli of the basal layer (red), the dermis (green), the epidermis (orange), and the vessels (white), respectively; $G_b$, $G_d$, $G_e$, and $G_v$ are the growth factors of the corresponding layers. No displacement is imposed at the sides and bottom boundaries.

(A) Left: simulation with bio-inspired skin properties ($E_b = 0.05$, $E_d = E_e = 0.01$, $E_v = 10$, $G_b = 1$, $G_e = 0.75$, $G_d = G_v = 0$, Poisson's ratio = 0.35) recapitulating the key events of the natural hierarchical folding in the rhinarial skin. Stage 0: initial configuration with smooth layers and a square lattice of cylindrical vessels right below the basal layer. Stage 1: the basal layer folds and forms inward cups (yellow arrowhead) between the stiffer vessels. Stage 2: the epidermal surface folds to form outward domes, surrounded by inward creases (black arrows) facing the vessels, and smaller dimples in the middle of the domes (black arrowheads). In stages 1 and 2, mechanical stress concentrates at the vessels (dashed white circles), and higher velocity occurs in the inward cups of the basal layer and the outward domes of the epidermis (white arrows).

(B) Simulation with parameters identical to (A), but without the vessels, generates an unnatural brainy pattern. The natural length scale (black scale bar) of the system simulated in (B) is smaller than the distance between the vessels (white scale bar).

(C–H) Transitions in folding regimes as crucial mechanical parameters are varied. Layers properties, except when otherwise stated, are $E_b = E_e = 1$, $E_d = 0.1$, $E_v = 10$, $G_b = G_e = 1$, $G_d = G_v = 0$, Poisson's ratio = 0.35. White arrows indicate transitions from a non-biological folding regime (epidermal creases between the vessels) to a biological folding regime (epidermal creases facing the vessels). The effects of varying the relative values of the parameters between layers (C–H) are discussed in the text.

localization of epidermal creases exactly over the underlying blood vessels is that the latter must be substantially stiffer (roughly >5 times) than the epidermal basal layer and upper epidermis (Figure 5C). This condition is very likely met in the rhinarium by the active regulation of blood pressure within the skin vasculature. Indeed, the resistance of vessels to deformation is imposed not only by the passive structural stiffness of the vessel walls but also by internal fluid forces and by the active contraction of smooth muscles wrapped inside their walls.⁶² The distributions of mechanical stress and velocity in these simulations show that the stress is focalized at the positions of vessels, forming rigid base points between which both the basal and upper epidermal layers move outward to shape domes, akin to arches rising against rigid pillars.

Strikingly, the outward epidermal bumps further develop central shallower dimples (black arrowheads in Figure 5A), reminiscent of the true geometry of epidermal domes in adult dogs (see inset of Figure 1B). In addition, we find that unnatural brainy folding patterns emerge as the stiffness of the basal layer is increased over 10 to 20 times that of the overlaying upper epidermis (Figure 5D). Note that this transition occurs irrespective of the absolute stiffness values of the two layers (Figure 5E), confirming that the folding regime is in fact determined by their relative stiffnesses. As for the dermis, its stiffness must be similar

to, or lower than, those of the epidermal layers in order to generate realistic folding patterns (Figure 5F), although the system transitions to a non-realistic serpentine folding regime for excessively soft dermis.

Unsurprisingly, natural folding of the basal layer requires that its tangential growth be at least as large as that of the upper epidermis (Figure 5G). Remarkably, the exact positioning of the epidermal surface creases over the blood vessels remains robust even under the period-doubling regime, which is observed when the growth ratio between basal and upper epidermis exceeds about three times (Figure 5G). Finally, growth of the dermis must be lower than those of the epidermal layers (Figure 5H); otherwise, a network of creases with free branches forms instead of individual domes surrounded by closed loops.

We investigate further the deterministic coupling between the pre-pattern of vessels and the epidermal creasing pattern by varying the geometrical distribution of vessels (Figure 6A) from an ideal hexagonal grid to a broken hexagonal grid (comprising free branches) or a noisy quasi-hexagonal grid (comprising pentagons but no free branches). The latter resembles more closely the real vascular networks of mammalian noses, allowing us to assess how spatial variations could influence the epidermal folding patterns. These analyses indicate that the key features of the folding process described above
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 6. Numerical simulations of growth in a flat-slab multilayered skin model with (near-) hexagonal networks of vessels**</figcaption></figure>

Epidermal basal layer is shown in red, upper epidermis in orange, dermis in green, and vessels in white. No displacement is imposed at the sides and bottom boundaries.

(A) Simulations with (i) an ideal hexagonal network of sub-epidermal vessels, (ii) a broken near-hexagonal network with free branches (open black circles), and (iii) and (iv) noisy hexagonal networks including pentagons. Mechanical parameters are identical to those of Figure 5A for panels (i), (ii), and (iii) and identical to the default parameters of Figures 5C–5H for panel (iv). The left panels show the pattern of vessels, whereas the right panels show the epidermal surface color-coded by the velocity magnitude at the surface. In all cases, the spatial distribution of epidermal creases corresponds to that of the vessels. Black arrowheads indicate smaller dimples in the middle of epidermal domes. Note that the simulation of panel (iv) stops at 65% because multiple polygons merge at steady state; polygons remain stable if vessels are made even stiffer.

(B) Chronological stages of hierarchical folding in the multilayered model from (A), panel (iii), illustrated in XY planes (top) and XZ planes (bottom). The same key events of the folding process occur regardless of the vessel network geometry (see also Figure 5A): the basal layer first forms inward cups between the stiffer vessels (stage 1), followed by the formation of outward domes at the surface (stage 2), surrounded by epidermal creases (black arrow) overlaying the vessels. Mechanical stress concentrates at the vessels (dashed white circles), and higher velocity occurs in the inward cups of the basal layer and the outward domes of the epidermis (white arrows).

(C) Light-sheet microscopy images of the rhinarial skin pattern of a juvenile ferret (P30) shown for comparison, in transversal and sagittal planes (top and bottom panels, respectively).

See also Figure S1 and Video S6.

(Figure 5A) are robust against variations of the geometry of the underlying pattern of vessels. As illustrated in Figure 6B, epidermal creases robustly form on top of the stiffer vessels, such that variations in the spatial distribution of the vascular pre-pattern are reflected by the distribution of surface folds, i.e., even in the presence of free branches and/or irregular polygons. The quasi-hexagonal model closely resembles the natural geometry of rhinarial skin folds, as shown by comparison with a juvenile ferret (Figure 6C). The correspondence between the simulated and real pattern is clear, with sharp epidermal creases facing the vessels, main inward cups in the basal layer (between the vessels), and secondary folds of the basal layer (over the vessels). Hence, the polygonal order of the folding pattern is imposed by the topology of the vessels.

The presence of secondary features of the folding pattern, such as central dimples in the epidermal domes (black arrowheads in Figures 5A and 6A), or the inward/outward direction of the folds in the basal layer (Figure 2B), can be regulated in simulations by varying mechanical properties. For example, the central dimples disappear when increasing skin stiffness (panel iv of Figure 6A), and the basal layer folds inward instead of outward when it grows faster than the upper epidermis (Figures 5A and 5G).

**Growth simulations in realistic volumetric segmentation models**

We also perform numerical simulations with a more realistic model of an E38 dog rhinarium reconstructed from the
<figure><img src="image_15.png" /><figcaption>**Figure 7.** Growth simulations in a multilayered skin model reconstructed from volumetric imaging of a dog rhinarium at E38</figcaption></figure>

(A) Initial configuration (dermis, epidermis, basal layer, blood vessels, and nasal cartilage in pink, beige, dark gray, red, and purple, respectively) after smoothing of the epidermal surface. Panel (i) shows a transversal section with vessels protruding from under in the epidermal basal layer. Panel (ii) shows the adaptive-size tetrahedral mesh of the 3D model. No displacement is imposed at the truncated boundaries.

(B) Final folding pattern obtained by setting the mechanical properties to biologically relevant values from the glabrous skin of 14-week-old human embryos³¹ (soft dermis, stiffer epidermis, and intermediately stiff epidermal basal layer: here set to $E_d = 0.1$, $E_e = 0.25$, $E_b = 0.15$), and other mechanical properties taken from our flat-slab simulations (stiffer vessels and cartilage: $E_v = 10$, $E_c = 50$; Poisson's ratio = 0.3). Relative growth rates follow our measurements of cellular proliferation ($G_b = G_e = 1.5$, $G_d = 0$). Panel (i) shows a close-up, color-coded by maximum surface curvature, that can be compared with the rhinarial pattern of a juvenile dog (panel (ii), MIP from light-sheet data; autofluorescence of keratinized epidermis). Panel (iii) shows that stress is focalized in the segmented vessels (dashed outlines), whereas velocities are high in the domes (white arrows).

(C) Final folding patterns obtained with the same parameters as in (B), with (left) or without (right) stiffer vessels. In the absence of positional information from stiffer vessels, the skin forms a non-biological serpentine pattern.

(D) Statistical distributions of geometrical attributes of the epidermal domes from the real patterns of two juvenile dogs (red histograms) versus the simulated patterns of the model with stiffer vessels (blue) and without vessels (yellow).
See Figures S2, S5, and S7.

volumetric segmentation of light-sheet microscopy data (Figure 7). Although the surface of the epidermis seems smooth at this embryonic stage, we nevertheless apply Gaussian smoothing to the epidermal envelope (left panel of Figure 7A) to erase any topological information from the surface that could influence the final simulation results. At this stage, the pre-pattern of blood vessels is present beneath the basal layer, and the latter has barely started to form inward folds, while the epidermal creases

have not yet appeared (top panel of Figure 4A). Using this realistic nose geometry, as well as a set of relative elastic properties measured by Glover et al.³¹ in the fingertip skin of human embryos (i.e., soft dermis, stiffer epidermis, and intermediate basal layer), our growth simulations mimic the formation of polygonal epidermal domes: panels (i) and (ii) in Figure 7B compare the morphology of the simulated folds (color-coded by maximum surface curvature) with that of a real dog embryonic rhinarium.
# Current Biology
Article

<figure><img src="image_16.png" /></figure>

Note that these simulations recapitulate the positioning of the surface creases facing the underlying template of stiffer vessels (Figure 7C, left). Panel (iii) of Figure 7B shows that, as in the simulations with the flat-slab models (Figures 5A and 6B), the stress builds up around the stiff vessels, while velocities are largest in the center of the growing epidermal domes. To demonstrate that the sub-epidermal stiffer vessels are essential to the folding process, we compare the results of folding simulations with stiff vessels (Figure 7C, left) with those generated when the stiffness of the vessels is set to the same value as the dermis (which is equivalent to ignoring the presence of vessels altogether). In the latter case, the resulting pattern is clearly unnatural, with serpentine units, while the vessels are deformed along with the surrounding dermis (Figure 7C, right).

We then quantify the topology of the modeled morphologies following the same approach as for the real rhinarial patterns (Figures 1 and S2). Using a curvature-based segmentation of the polygonal domes, we generate statistical distributions of selected geometrical attributes, namely their area, AR, polygon order, and mean height (Figure 7D). This quantification was restricted to the upper half of the nose model, where the 3D reconstruction of the vessels (Figure 4A) is most accurate. Note that all lengths are normalized by the distance between the nostrils. In Figure 7D, we compare the distributions of these geometrical features between natural patterns (from two juvenile dog samples) and the simulated phenotypes (model with stiffer vessels and model without vessels, as shown in Figure 7C). We find that the natural distribution of ARs is fairly well reproduced by the model with stiffer vessels, whereas the tail of the distribution increases toward larger values for the serpentine pattern that develops in the absence of vessels. Both models result in larger unit areas compared with the natural pattern, with an even more drastic increase without the vessels. One plausible explanation for this discrepancy is that the global geometry of the model is static during the simulated skin growth, whereas the natural skin patterning occurs concurrently with (1) the overall growth of the rhinarium and (2) the remodeling of the vascular plexus. Both of these processes are likely to yield further unit subdivisions (explaining the smaller relative sizes of the units in real rhinaria).

The distributions of polygonal order in the real and simulated patterns indicate that the model with stiffer vessels contains a slight bias toward hexagonal units, even though it also contains a large number of pentagons. This is most likely due to an incomplete segmentation of the network of vessels (Figure 4A). Indeed, smaller branches could hardly be detected by the segmentation algorithm (see STAR Methods), and the local mechanical coupling between the folding skin and the underlying vessels is thus missing in those places. Finally, the epidermal domes in the natural organ appear more flat than the simulated ones (Figures 7B and 7C). This difference is reflected in the height distributions (Figure 7D), leaning toward much larger values for the simulated phenotypes (even more so for the model without vessels). This result is not surprising for two reasons. First, the keratinization of the epidermal surface during embryogenesis likely contributes to flattening the polygonal domains in real embryos. Second, the samples were dehydrated during their preparation for light-sheet microscopy and histology (see STAR Methods), which cause tissue shrinkage and additional surface flattening.

Overall, our simulations with both models (flat slab with simplified grid of blood vessels and realistic morphology segmented from light-sheet imaging data) indicate that the lattice of polygonal skin domains of the rhinarium forms through a self-organized process dominated by mechanical folding under compressive stress induced by differential growth rates among skin layers. Importantly, the localization of the skin creases is strongly constrained by the presence of a circuit of stiffer vessels, i.e., the growing multilayered skin harnesses the pre-existing positional information encoded in the vascular network to constrain the positions of the creases. These results pertain to a concept of *mechanical positional information*, whereby material properties of anatomical elements mechanically encode positional information without requiring chemical signaling.

## Stochasticity of rhinarial skin patterns

Given that vasculogenesis (the assembly of an early polygonal vascular network) likely emerges from the self-assembly of endothelial cells,⁵⁹ the corresponding positions of rhinarial skin creases are more likely to be stochastically distributed than genetically determined. To test this hypothesis, we statistically analyze the spatial distributions of rhinoglyphics in multiple individuals of genetically identical cows (i.e., genetic clones) versus unrelated controls (as well as dog and ferret individuals for comparison). To do so, we first generate 3D reconstructions of the rhinarial surfaces of all selected individuals using structure-from-motion,⁶³ structured light, or light-sheet microscopy, depending on the size of the samples (see STAR Methods). We then segment the acquired surfaces into individual polygonal domes using a curvature-based approach (Figures 1, S2A, and S6A) and quantify the statistical distributions of various geometric attributes, such as polygon length, height, orientation, and number of nodes (Figures 1 and S2B). Finally, we measure the Jensen-Shannon divergence ($D_{JS}$) between these statistical distributions to quantify the level of disparity between the patterns of different samples (see STAR Methods). The $D_{JS}$ maps of all geometrical attributes are averaged together and over groups of equivalent samples (e.g., bovine clones or controls) to obtain a single parameter normalized between 0 (statistically identical patterns) and 1 (maximally different patterns), as shown in Figures S6B and S6C. For example, Figure S6C shows that the AR of polygonal domes distinguishes well between cows and the other two species (high $D_{JS}$ values), but not between ferrets and dogs (low $D_{JS}$ values). In contrast, the number of branch nodes associated with each polygon hardly discriminates the lattices of the different sample groups. The $D_{JS}$ maps also demonstrate that ferrets display the most divergent patterns within the three compared species, especially in terms of the relative dimensions of epidermal domes (width, length, and height).

Most importantly, these analyses indicate that the rhinarial pattern of creases is not genetically determined as mean $D_{JS}$ values (Figure S6B) among bovine clones and among unrelated individuals are similarly small, i.e., the rhinoglyphics of genetically identical animals exhibit the same overall level of stochasticity as unrelated individuals. In addition, Kolmogorov-Smirnov tests comparing distributions of geometric attributes confirm that the patterns of cow clones are not significantly different from those of unrelated individuals ($p$ values >3.5% for all geometrical attributes, except for the height of epidermal
domes; some differences in sampling and preservation of the samples might have caused a post-mortem deformation of the surfaces due to dehydration, therefore introducing a bias in the height distributions for each sample group). Conversely, a strong statistical divergence is observed between the different species, while the patterns are homogenous within each species (Figure S6B), suggesting that the parameters controlling self-organized vasculogenesis/angiogenesis are substantially different among these divergent species.

Finally, we analyzed the rhinarial patterns of cows affected by the Ehlers-Danlos syndrome, a hereditary collagenopathy that impacts the mechanical properties of connective tissues, resulting in hyperextensibility of the skin. These samples further demonstrate that mechanical properties are determinant in shaping the network of skin folds. Indeed, their patterns clearly differ from those of healthy cows (clones and controls), as indicated by high $D_{JS}$ values; actually, they are more similar to the patterns of dogs (Figure S6B) than to those of healthy cows. Specifically, cows with Ehlers-Danlos syndrome differ from healthy animals mainly in terms of the length, width, and AR of their polygonal units (Figure S6C). This is due to the fact that the network of creases presents many unterminated branches (pink dots in Figure S6A), although the global orientation of the epidermal units seems to be preserved in these COL5A1 mutants. By comparing Ehlers-Danlos mutants with cow controls using a Kolmogorov-Smirnov test, we obtain low *p* values (<1%) for the height, length, width, area, and AR of the epidermal domes, confirming that the patterns are statistically distinct. However, the *p* values obtained for the number of branch and end nodes are high (around 30%) because the disrupted patterns of those two individuals are characterized by the fusion of several domes into a single, very large unit containing many nodes (Figure S6B), therefore not contributing much to the statistical significance for that geometrical attribute. Based on our growth simulations on the flat-slab model, which showed that too soft dermis results in the incomplete formation of bordering creases around the polygonal domes (right-most panel in Figure 5F), we suggest that the defective pattern of Ehlers-Danlos cows is caused by improper mechanical folding of the epidermal layers over the sub-epidermal vascular template, due to insufficient dermal stiffness.

DISCUSSION

Our analyses indicate that the ontogenic patterning of rhinarial skin creases is established through a growth-induced self-organized mechanical folding process. However, the dynamics and the resulting steady state are strongly constrained by pre-existing positional information bore by a sub-epidermal vascular network and transferred to the creasing epidermis through mechanical (instead of chemical) interactions, i.e., in the form of a stress field. Hence, akin to the principle of (chemical) positional information<sup>8,57</sup> involved in multiple morphogenetic processes, here we present a case of *mechanical positional information* determined by specific anatomical components. Instead of involving chemical signaling, this principle relies on positional information being mechanically encoded by the specific material properties (e.g., stiffness) of the corresponding anatomical elements.

An outstanding feature of the mechanical patterning of the rhinarium is the hierarchical sequence in which the basal layer of the epidermis first folds between the vessels due to excessive proliferation, followed by the folding of the epidermis into outward domes surrounded by sharp creases facing the lines of the vessels. Our numerical simulations confirm that these basic principles of morphogenesis, dominated by mechanical parameters, can be effectively modeled without involving any biochemical interactions. Here we use simple numerical models with isotropic mechanical properties to recapitulate the essential folding steps leading to rhinarial skin patterning. In reality, the complex architecture of the ECM could induce anisotropic responses to growth stresses. In particular, the network of collagen fibers displays a preferred orientation perpendicularly to the skin plane in embryonic dog rhinaria, as revealed by light-sheet microscopy images (Figure S7). We predict that encompassing this type of anisotropy in the elastic model would generate yet more realistic folding patterns, for example, by limiting the strain perpendicularly to the skin plane and thus reducing the curvature of the modeled epidermal domes, which appear more flat in natural phenotypes. Moreover, here we considered the epidermal basal layer (red layer in Figure 5) to effectively represent the basal layer plus its basement membrane. The latter has been shown to play an important mechanical role in the development of skin cancer.<sup>64-68</sup> Therefore, refining the anatomical resolution of our models by distinguishing basal epidermal layer and basement membrane could also improve the accuracy of the simulated folding process.

It is likely that the pre-pattern of blood vessels mechanically transferred to the folding epidermis provides an additional level of robustness to the morphogenesis of *functional rhinoglyphics*. Indeed, unconstrained self-organized skin folding (i.e., relying solely on differential growth rates between different layers of the skin) would generate a network of creases that would not coincide with the spatial distribution of vessels, hence making rhinoglyphics less efficient for functions such as evaporative cooling or nose temperature control. In other words, evolution harnessed a mechanical coupling between subcutaneous blood vessels and the embryonic rhinarium skin, thereby ensuring that skin folds coincide (spatially) with the self-organized vasculogenesis process. This point is spectacularly illustrated by our observations in cows: genetically identical individuals do not exhibit more similar nose folds than unrelated individuals because the positions of the vessels are not genetically controlled but result from self-organized vasculogenesis. Note that self-organized vasculogenesis must precede mechanical self-organized skin folding to ensure the robust emergence of skin creases at the exact position of vessels. On the other hand, divergent evolution of the vasculogenesis parameters explains why different species (such as cow, dog, and ferret) exhibit different rhinoglyphics phenotypes. Note that the relative growth rates and material properties of the skin layers (dermis, epidermal basal layer, and upper epidermis) must have coevolved to maintain an effective mechanical coupling of the folding process with the subcutaneous vascular network.

Early on, morphological resemblance between rhinarial folds and the papillary ridges (fingerprints) found in the glabrous skin of human hands has been noted.<sup>40</sup> Similarly to the rhinoglyphics whose ontogeny is revealed here, fingerprints have been
suggested to initially develop as undulations in the epidermal basal layer, presetting the subsequent formation of epidermal ridges.⁶⁹,⁷⁰ Fingerprints have long been proposed to originate from growth-induced compressive stress,³⁰,⁷¹⁻⁷³ an hypothesis compatible with the fact that the final pattern of swirls and loops correlates with the surface curvature of the embryonic fingertip.⁷⁴ Moreover, it has been suggested that sub-epidermal nerves or blood vessels in fingertips could somehow accompany the formation of folds in the basal layer and upper epidermis.³⁰,⁷⁰,⁷⁵⁻⁸⁰ These observations are compatible with fingerprints emerging from a self-organized folding process similar to that of rhinoglyphics uncovered here, i.e., constrained by mechanical positional information originating from underlying stiff vessels.

Alternatively, fingerprints might be produced by a fundamentally different ontogenic patterning process. Indeed, more recent analyses of fingerprint development propose that Turing's reaction-diffusion—involving WNT as an activator and BMPs as inhibitors—drives focalized epidermal proliferation, thereby triggering the formation of early fingerprints as U-shaped epidermal buds.³¹ Thus, fingerprint patterning might emerge from a chemical self-organized process (reaction-diffusion), possibly constrained by chemical positional information provided by anatomical elements such as blood vessels or nerves, whereas we show that rhinoglyphics develop through a mechanical self-organized process constrained by mechanical positional information provided by an early network of subcutaneous blood vessels.

Although the biochemical and mechanical paradigms are not mutually exclusive (patterning of skin folds could involve dynamic interactions between a growth-induced mechanical stress field and templates of morphogen concentrations), we show here that the development of a polygonal network of skin creases in the rhinarium of multiple mammalian species is dominated by a mechanical self-organized folding process. Additionally, we show for the first time that a pre-existing network of vessels can modulate the growth-induced stress field and thereby regulate robust patterning of epidermal skin folds. The principle of mechanical positional information that we develop here opens new perspectives on the emergence of multiple intricate geometries whose morphogenesis is dominated by mechanics.

## RESOURCE AVAILABILITY

### Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Michel C. Milinkovitch (michel.milinkovitch@unige.ch).

### Materials availability

This study did not generate new materials.

### Data and code availability

* All data reported in this paper will be shared by the lead contact upon request.

* Original code for the differential growth simulations has been deposited at CodeOcean (capsule with the source code, examples of input meshes, and simulations readily executable) and is publicly available as of the date of publication. The link is provided in the key resources table.

* Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

## ACKNOWLEDGMENTS

We thank Grigorii Timin, Rory Cooper (LANE), and Ivana Gantar (Wyss Center) for advice on light-sheet microscopy protocols; Valérie Haechler (LANE) for technical assistance; Szabolcs Zakany (LANE) for advice on structure-from-motion 3D reconstructions; Adrien Debry and Florent Montange (LANE) for sample collection; and Athanasia Tzika, Carine Langrez, and Ingrid Burgelin (LANE) for advice on histology. We thank Yves Millemann (ENVA), Sarah Barbey, Mekki Boussaha, Cécile Grohs, Frédéric Launay (INRAE), François Lejuste (Eliance), the staff of the INRAE experimental units of Bressonvilliers and Le-Pin-au-Haras, the Abattoir SOCOPA (Gacé, France), and the Abattoir des Marais (Estavayer-le-Lac, Switzerland) for their help in obtaining and sampling embryonic and adult bovine muzzle specimens. We thank Sandra Thoumire, the CERCA (Centre d'Etude en Reproduction des Carnivores), and the clinicians at the Alfort Veterinary School for their help in the dog oestrus follow-up, inseminations, and collection of prenatal tissues. Funding: this work was supported by grants to M.C.M. from the Georges & Antoine CLARAZ Foundation, the Swiss National Science Foundation (FNSNF; grant 31003A_179431), the International Human Frontier Science Program Organization (HFSP RGP0019/2017), and the European Research Council (ERC; advanced grant EVOMORPHYS) under the European Union's Horizon 2020 research and innovation program. V.B. was supported by an ERC Synergy Grant (UNFOLD 101118729). The funding bodies played no role in the design of the study; collection, analysis, and interpretation of data; or writing of the manuscript.

## AUTHOR CONTRIBUTIONS

M.C.M. conceived the study and supervised all aspects of the research. A.C. and H.J. sampled and provided cow material. K.R. sampled and provided dog material. C.D.J.R. and V.B. sampled and provided ferret material. P.D. performed gel analog experiments, rhinarium tissue preparation, cell proliferation assays, histology, and volumetric imaging. P.D. and E.J. performed numerical simulations and statistical analyses. P.D. and M.C.M. analyzed and interpreted all data. M.C.M. and P.D. wrote the manuscript. All authors read and approved the final manuscript.

## DECLARATION OF INTERESTS

The authors declare no competing interests.

## STAR★METHODS

Detailed methods are provided in the online version of this paper and include the following:

* KEY RESOURCES TABLE

* EXPERIMENTAL MODEL AND SUBJECT DETAILS

* METHOD DETAILS

* Biological samples

* Histology

* Light-sheet microscopy

* Hydrogel analogues

* QUANTIFICATION AND STATISTICAL ANALYSIS

* 3D segmentation of tissues

* 3D surface reconstruction of large samples

* Differential growth numerical simulations

* Statistical analysis of polygonal networks

## SUPPLEMENTAL INFORMATION

Supplemental information can be found online at https://doi.org/10.1016/j.cub.2024.09.055.

Received: January 24, 2024

Revised: July 5, 2024

Accepted: September 20, 2024

Published: October 22, 2024
REFERENCES

1. Turing, A.M. (1952). The Chemical Basis of Morphogenesis. Philos. Trans. R. Soc. Lond. B 237, 37-72.

2. Milinkovitch, M.C., Jahanbakhsh, E., and Zakany, S. (2023). The Unreasonable Effectiveness of Reaction Diffusion in Vertebrate Skin Color Patterning. Annu. Rev. Cell Dev. Biol. 39, 145-174. https://doi.org/10.1146/annurev-cellbio-120319-024414.

3. Musser, J.M., Wagner, G.P., and Prum, R.O. (2015). Nuclear -catenin localization supports homology of feathers, avian scutate scales, and alligator scales in early development. Evol. Dev. 17, 185-194. https://doi.org/10.1111/ede.12123.

4. Di-Poï, N., and Milinkovitch, M.C. (2016). The anatomical placode in reptile scale morphogenesis indicates shared ancestry among skin appendages in amniotes. Sci. Adv. 2, e1600708. https://doi.org/10.1126/sciadv.1600708.

5. Cooper, R.L., Martin, K.J., Rasch, L.J., and Fraser, G.J. (2017). Developing an ancient epithelial appendage: FGF signalling regulates early tail denticle formation in sharks. EvoDevo 8, 8. https://doi.org/10.1186/s13227-017-0071-0.

6. Ho, W.K.W., Freem, L., Zhao, D., Painter, K.J., Woolley, T.E., Gaffney, E.A., McGrew, M.J., Tzika, A., Milinkovitch, M.C., Schneider, P., et al. (2019). Feather arrays are patterned by interacting signalling and cell density waves. PLoS Biol. 17, e3000132. https://doi.org/10.1371/journal.pbio.3000132.

7. Tallinen, T., and Biggins, J.S. (2015). Mechanics of invagination and folding: Hybridized instabilities when one soft tissue grows on another. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 92, 022720. https://doi.org/10.1103/PhysRevE.92.022720.

8. Wolpert, L. (1969). Positional information and the spatial pattern of cellular differentiation. J. Theor. Biol. 25, 1-47. https://doi.org/10.1016/s0022-5193(69)80016-0.

9. Nelson, C.M., Jean, R.P., Tan, J.L., Liu, W.F., Sniadecki, N.J., Spector, A.A., and Chen, C.S. (2005). Emergent patterns of growth controlled by multicellular form and mechanics. Proc. Natl. Acad. Sci. USA 102, 11594-11599. https://doi.org/10.1073/pnas.0502575102.

10. Brinkmann, F., Mercker, M., Richter, T., and Marciniak-Czochra, A. (2018). Post-Turing tissue pattern formation: Advent of mechanochemistry. Plos Comput. Biol. 14, e1006259. https://doi.org/10.1371/journal.pcbi.1006259.

11. Bailles, A., Gehrels, E.W., and Lecuit, T. (2022). Mechanochemical Principles of Spatial and Temporal Patterns in Cells and Tissues. Annu. Rev. Cell Dev. Biol. 38, 321-347. https://doi.org/10.1146/annurev-cell-bio-120420-095337.

12. Ben Amar, M., and Ciarletta, P. (2010). Swelling instability of surface-attached gels as a model of soft tissue growth under geometric constraints. J. Mech. Phys. Solids 58, 935-954. https://doi.org/10.1016/j.jmps.2010.05.002.

13. Cao, Y.P., and Hutchinson, J.W. (2012). From wrinkles to creases in elastomers: the instability and imperfection-sensitivity of wrinkling. Proc. R. Soc. A. 468, 94-115. https://doi.org/10.1098/rspa.2011.0384.

14. Tallinen, T., Biggins, J.S., and Mahadevan, L. (2013). Surface Sulci in Squeezed Soft Solids. Phys. Rev. Lett. 110, 024302. https://doi.org/10.1103/PhysRevLett.110.024302.

15. Wang, Q., and Zhao, X. (2015). A three-dimensional phase diagram of growth-induced surface instabilities. Sci. Rep. 5, 8887. https://doi.org/10.1038/srep08887.

16. Tanaka, T., Sun, S.T., Hirokawa, Y., Katayama, S., Kucera, J., Hirose, Y., and Amiya, T. (1987). Mechanical Instability of Gels at the Phase-Transition. Nature 325, 796-798. https://doi.org/10.1038/325796a0.

17. Komori, T., Takahashi, H., and Okamoto, N. (1988). Swelling of Sodium Acrylate Gels. Colloid Polym. Sci. 266, 1181-1186. https://doi.org/10.1007/BF01414408.

18. Mora, T., and Boudaoud, A. (2006). Buckling of swelling gels. Eur. Phys. J. E Soft Matter 20, 119-124. https://doi.org/10.1140/epje/i2005-10124-5.

19. Sultan, E., and Boudaoud, A. (2008). The buckling of a swollen thin gel layer bound to a compliant substrate. Journal of Applied Mechanics 75, https://doi.org/10.1115/1.2936922.

20. Trujillo, V., Kim, J., and Hayward, R.C. (2008). Creasing instability of surface-attached hydrogels. Soft Matter 4, 564-569. https://doi.org/10.1039/b713263h.

21. Dervaux, J., and Ben Amar, M. (2011). Buckling condensation in constrained growth. J. Mech. Phys. Solids 59, 538-560. https://doi.org/10.1016/j.jmps.2010.12.015.

22. Dervaux, J., and Ben Amar, M. (2012). Mechanical Instabilities of Gels. Annu. Rev. Condens. Matter Phys. 3, 311-332. https://doi.org/10.1146/annurev-conmatphys-062910-140436.

23. Stoop, N., Lagrange, R., Terwagne, D., Reis, P.M., and Dunkel, J. (2015). Curvature-induced symmetry breaking determines elastic surface patterns. Nat. Mater. 14, 337-342. https://doi.org/10.1038/Nnat4202.

24. Chan, H.F., Zhao, R.K., Parada, G.A., Meng, H., Leong, K.W., Griffith, L.G., and Zhao, X.H. (2018). Folding artificial mucosa with cell-laden hydrogels guided by mechanics models. Proc. Natl. Acad. Sci. USA 115, 7503-7508. https://doi.org/10.1073/pnas.1802361115.

25. Tan, Y.L., Hu, B.R., Song, J., Chu, Z.Y., and Wu, W.J. (2020). Bioinspired Multiscale Wrinkling Patterns on Curved Substrates: An Overview. Nanomicro. Lett. 12, 101. https://doi.org/10.1007/s40820-020-00436-y.

26. Wang, L.F., Castro, C.E., and Boyce, M.C. (2011). Growth strain-induced wrinkled membrane morphology of white blood cells. Soft Matter 7, 11319-11324. https://doi.org/10.1039/c1sm06637d.

27. Feng, L., Zhang, Y.A., Xi, J.M., Zhu, Y., Wang, N., Xia, F., and Jiang, L. (2008). Petal effect: A superhydrophobic state with high adhesive force. Langmuir 24, 4114-4119. https://doi.org/10.1021/la703821h.

28. Lugo, C.A., Airoldi, C., Chen, C., Crosby, A.J., and Glover, B.J. (2023). Morphoelastic modelling of pattern development in the petal epidermal cell cuticle. J. R. Soc. Interface 20, 20230001. https://doi.org/10.1098/rsif.2023.0001.

29. Huang, C.J., Quinn, D., Suresh, S., and Hsia, K.J. (2018). Controlled molecular self-assembly of complex three-dimensional structures in soft materials. Proc. Natl. Acad. Sci. USA 115, 70-74. https://doi.org/10.1073/pnas.1717912115.

30. Kücken, M., and Newell, A.C. (2004). A model for fingerprint formation. Europhys. Lett. 68, 141-146. https://doi.org/10.1209/epl/i2004-10161-2.

31. Glover, J.D., Sudderick, Z.R., Shih, B.B.J., Batho-Samblas, C., Charlton, L., Krause, A.L., Anderson, C., Riddell, J., Balic, A., Li, J., et al. (2023). The developmental basis of fingerprint pattern formation and variation. Cell 186, 940-956.e20. https://doi.org/10.1016/j.cell.2023.01.015.

32. Shyer, A.E., Tallinen, T., Nerurkar, N.L., Wei, Z.Y., Gil, E.S., Kaplan, D.L., Tabin, C.J., and Mahadevan, L. (2013). Villification: How the Gut Gets Its Villi. Science 342, 212-218. https://doi.org/10.1126/science.1238842.

33. Dervaux, J., Couder, Y., Guedeau-Boudeville, M.A., and Ben Amar, M. (2011). Shape Transition in Artificial Tumors: From Smooth Buckles to Singular Creases. Phys. Rev. Lett. 107, 018103. https://doi.org/10.1103/PhysRevLett.107.018103.

34. Risler, T., and Basan, M. (2013). Morphological instabilities of stratified epithelia: a mechanical instability in tumour formation. New J. Phys. 15, https://doi.org/10.1088/1367-2630/15/6/065011.

35. Tallinen, T., Chung, J.Y., Biggins, J.S., and Mahadevan, L. (2014). Gyrfication from constrained cortical expansion. Proc. Natl. Acad. Sci. USA 111, 12667-12672. https://doi.org/10.1073/pnas.1406015111.

36. Tallinen, T., Chung, J.Y., Rousseau, F., Girard, N., Lefèvre, J., and Mahadevan, L. (2016). On the growth and form of cortical convolutions. Nature Phys. 12, 588-593. https://doi.org/10.1038/Nphys3632.

37. Dumais, J. (2007). Can mechanics control pattern formation in plants? Curr. Opin. Plant Biol. 10, 58-62. https://doi.org/10.1016/j.pbi.2006.11.014.
38. Milinkovitch, M.C., Manukyan, L., Debry, A., Di-Poï, N., Martin, S., Singh, D., Lambert, D., and Zwicker, M. (2013). Crocodile Head Scales Are Not Developmental Units But Emerge from Physical Cracking. *Science* 339, 78–81. https://doi.org/10.1126/science.1226265.

39. Rodríguez-Hernández, J. (2015). Wrinkled interfaces: Taking advantage of surface instabilities to pattern polymer surfaces. *Prog. Polym. Sci.* 42, 1–41. https://doi.org/10.1016/j.progpolymsci.2014.07.008.

40. Hill, B.W.C.O. (1948). Rhinoglyphics: Epithelial sculpture of the Mammalian Rhiarium. *Proc. Zool. Soc. Lond.* 118, 1–35. https://doi.org/10.1111/j.1096-3642.1948.tb00364.x.

41. Gläser, N., and Kröger, R.H.H. (2017). Variation in rhinarium temperature indicates sensory specializations in placental mammals. *J. Therm. Biol.* 67, 30–34. https://doi.org/10.1016/j.jtherbio.2017.04.010.

42. Tuminaite, I., and Kröger, R.H.H. (2021). Rhinarium skin structure and epidermal innervation in selected mammals. *J. Morphol.* 282, 419–426. https://doi.org/10.1002/jmor.21313.

43. Caputa, M. (2004). Selective brain cooling: a multiple regulatory mechanism. *J. Therm. Biol.* 29, 691–702. https://doi.org/10.1016/j.jtherbio.2004.08.079.

44. Blatt, C.M., Taylor, C.R., and Habal, M.B. (1972). Thermal Panting in Dogs - Lateral Nasal Gland, A Source of Water for Evaporative Cooling. *Science* 177, 804–805. https://doi.org/10.1126/science.177.4051.804.

45. Tuminaite, I. (2020). Functional morphology of mammalian rhinarium skin. *Doctoral Thesis* (Lunds Universitet).

46. Ramser, R. (1935). Zur Anatomie des Jakobson'schen Organs Beim Hunde (Tierarztliche Hochschule zu Berlin).

47. Kelliher, K.R., Baum, M.J., and Meredith, M. (2001). The ferret's vomeronasal organ and accessory olfactory bulb: Effect of hormone manipulation in adult males and females. *Anat. Rec.* 263, 280–288. https://doi.org/10.1002/ar.1097.

48. Hillenius, W.J., and Rehorek, S.J. (2005). From the Eye to the Nose: Ancient Orbital to Vomeronasal Communication in Tetrapods? Held in Boston, MA (Springer), pp. 228–241.

49. Rossie, J.B., and Smith, T.D. (2007). Ontogeny of the nasolacrimal duct in primates: functional and phylogenetic implications. *J. Anat.* 210, 195–208. https://doi.org/10.1111/j.1469-7580.2006.00682.x.

50. Ortiz-Leal, I., Torres, M.V., Villamayor, P.R., López-Beceiro, A., and Sanchez-Quintiro, P. (2020). The vomeronasal organ of wild canids: the fox (Vulpes vulpes) as a model. *J. Anat.* 237, 890–906. https://doi.org/10.1111/joa.13254.

51. Kröger, R.H.H., and Goiricelaya, A.B. (2017). Rhinarium temperature dynamics in domestic dogs. *J. Therm. Biol.* 70, 15–19. https://doi.org/10.1016/j.jtherbio.2017.10.013.

52. Martins, A.F., Bennett, N.C., Clavel, S., Groenewald, H., Hensman, S., Hoby, S., Joris, A., Manger, P.R., and Milinkovitch, M.C. (2018). Locally-curved geometry generates bending cracks in the African elephant skin. *Nat. Commun.* 9, 3865. https://doi.org/10.1038/s41467-018-06257-3.

53. Bálint, A., Andics, A., Gácsi, M., Gábor, A., Czeibert, K., Luce, C.M., Miklósi, Á., and Kröger, R.H.H. (2020). Dogs can sense weak thermal radiation. *Sci. Rep.* 10, 3736. https://doi.org/10.1038/s41598-020-60439-y.

54. Barker, D.J., and Welker, W.I. (1969). Receptive fields of first-order somatic sensory neurons innervating rhinarium in coati and raccoon. *Brain Res.* 14, 367–386. https://doi.org/10.1016/0006-8993(69)90116-4.

55. Catania, K.C. (2000). Epidermal sensory organs of moles, shrew-moles, and desmans: A study of the family talpidae with comments on the function and evolution of Eimer's organ. *Brain Behav. Evol.* 56, 146–174. https://doi.org/10.1159/000047201.

56. Elofsson, R., Tuminaite, I., and Kröger, R.H.H. (2015). A complex sensory organ in the nose skin of the prosimian primate Lemur catta. *J. Morphol.* 276, 649–656. https://doi.org/10.1002/jmor.20363.

57. Wolpert, L. (1971). Positional information and pattern formation. *Curr. Top. Dev. Biol.* 6, 183–224. https://doi.org/10.1016/s0070-2153(08)60641-9.

58. Goodwin, K., Mao, S., Guyomar, T., Miller, E., Radisky, D.C., Kosmrlj, A., and Nelson, C.M. (2019). Smooth muscle differentiation shapes domain

branches during mouse lung development. *Development* 146, dev181172. https://doi.org/10.1242/dev.181172.

59. Czirok, A., and Little, C.D. (2012). Pattern formation during vasculogenesis. *Birth Defects Res. C Embryo Today* 96, 153–162. https://doi.org/10.1002/bdrc.21010.

60. Guimarães, C.F., Gasperini, L., Marques, A.P., and Reis, R.L. (2020). The stiffness of living tissues and its implications for tissue engineering. *Nat. Rev. Mater.* 5, 351–370. https://doi.org/10.1038/s41578-019-0169-1.

61. Bowden, N., Brittain, S., Evans, A.G., Hutchinson, J.W., and Whitesides, G.M. (1998). Spontaneous formation of ordered structures in thin films of metals supported on an elastomeric polymer. *Nature* 393, 146–149. https://doi.org/10.1038/30193.

62. Camasão, D.B., and Mantovani, D. (2021). The mechanical characterization of blood vessels and their substitutes in the continuous quest for physiological-relevant performances. A critical review. *Mater. Today Bio* 10, 100106. https://doi.org/10.1016/j.mtbio.2021.100106.

63. Martins, A.F., Bessant, M., Manukyan, L., and Milinkovitch, M.C. (2015). R2OBBIE-3D, a Fast Robotic High-Resolution System for Quantitative Phenotyping of Surface Geometry and Colour-Texture. *PLoS ONE* 10, e0126740. https://doi.org/10.1371/journal.pone.0126740.

64. Chang, J., and Chaudhuri, O. (2019). Beyond proteases: Basement membrane mechanics and cancer invasion. *J. Cell Biol.* 218, 2456–2469. https://doi.org/10.1083/jcb.201903066.

65. Fiore, V.F., Krajnc, M., Quiroz, F.G., Levorse, J., Pasolli, H.A., Shvartsman, S.Y., and Fuchs, E. (2020). Mechanics of a multilayer epithelium instruct tumour architecture and function. *Nature* 585, 433–439. https://doi.org/10.1038/s41586-020-2695-9.

66. Winkler, J., Abisoye-Ogunniyan, A., Metcalf, K.J., and Werb, Z. (2020). Concepts of extracellular matrix remodelling in tumour progression and metastasis. *Nat. Commun.* 11, 5120. https://doi.org/10.1038/s41467-020-18794-x.

67. Reuten, R., Zendehroud, S., Nicolau, M., Fleischhauer, L., Laitala, A., Kiderlen, S., Nikodemus, D., Wullkopf, L., Nielsen, S.R., McNeilly, S., et al. (2021). Basement membrane stiffness determines metastases formation. *Nat. Mater.* 20, 892–903. https://doi.org/10.1038/s41563-020-00894-0.

68. Banerjee, S., Lo, W.C., Majumder, P., Roy, D., Ghorai, M., Shaikh, N.K., Kant, N., Shekhawat, M.S., Gadekar, V.S., Ghosh, S., et al. (2022). Multiple roles for basement membrane proteins in cancer progression and EMT. *Eur. J. Cell Biol.* 101, 151220. https://doi.org/10.1016/j.ejcb.2022.151220.

69. Penrose, L.S., and Ohara, P.T. (1973). Development of Epidermal Ridges. *J. Med. Genet.* 10, 201–208. https://doi.org/10.1136/jmg.10.3.201.

70. Seidenberg-Kajabova, H., Pospisilova, V., Vranakova, V., and Varga, I. (2010). An Original Histological Method for Studying the Volar Skin of the Fetal Hands and Feet. *Biomed. Pap. Med. Fac. Univ. Palacky Olomouc Czech. Repub.* 154, 211–218. https://doi.org/10.5507/bp.2010.032.

71. Bonnevie, K. (1929). Zur Mechanik der Papillarmusterbildung I. Wilhelm Roux' Archiv für Entwicklungsmechanik der Organismen 117, 384–420. https://doi.org/10.1007/BF02110970.

72. Bonnevie, K. (1932). Zur Mechanik der Papillarmusterbildung II. Wilhelm Roux' Archiv für Entwicklungsmechanik der Organismen 126, 348–372. https://doi.org/10.1007/BF00576269.

73. Kücken, M., and Newell, A.C. (2005). Fingerprint formation. *J. Theor. Biol.* 235, 71–83. https://doi.org/10.1016/j.tbi.2004.12.020.

74. Cummins, H. (1926). Epidermal-ridge configurations in developmental defects, with particular reference to the ontogenetic factors which condition ridge direction. *Am. J. Anat.* 38, 89–151. https://doi.org/10.1002/aja.1000380105.

75. Bonnevie, K. (1927). Die ersten Entwicklungsstadien der Papillarmuster der menschlichen Fingerballen. *Nyt Mag Natur* 65, 1956.

76. Hirsch, W., and Schweichel, J.U. (1973). Morphological Evidence Concerning Problem of Skin Ridge Formation. *J. Ment. Defic. Res.* 17, 58–72. https://doi.org/10.1111/j.1365-2788.1973.tb01186.x.
77. Okajima, M. (1982). A methodological approach to the development of epidermal ridges viewed on the dermal surface of fetuses. Prog. Clin. Biol. Res. 84, 175-188.

78. Dell, D.A., and Munger, B.L. (1986). The Early Embryogenesis of Papillary (Sweat-Duct) Ridges in Primate Glabrous Skin - the Dermatotopic Map of Cutaneous Mechanoreceptors and Dermatoglyphics. J. Comp. Neurol. 244, 511-532. https://doi.org/10.1002/cne.902440408.

79. Moore, S.J., and Munger, B.L. (1989). The Early Ontogeny of the Afferent Nerves and Papillary Ridges in Human Digital Glabrous Skin. Brain Res. Dev. Brain Res. 48, 119-141. https://doi.org/10.1016/0165-3806(89) 90097-7.

80. Wertheim, K. (2011). Embryology, Physiology, and Morphology of Friction Ridge Skin. In Fingerprint Sourcebook, N.I.o. Justice, ed. (NCJRS Virtual Library).

81. Bauer, A., Bateman, J.F., Lamandé, S.R., Hanssen, E., Kirejczyk, S.G.M., Yee, M., Ramiche, A., Jagannathan, V., Welle, M., Leeb, T., and Bateman, F.L. (2019). Identification of Two Independent COL5A1 Variants in Dogs with Ehlers-Danlos Syndrome. Genes-Basel 10, https://doi.org/10.3390/ genes10100731.

82. Jacinto, J.G.P., Häfliger, I.M., Veiga, I.M.B., Letko, A., Benazzi, C., Bolcato, M., and Drögemüller, C. (2020). A Heterozygous Missense Variant in the COL5A2 in Holstein Cattle Resembling the Classical Ehlers-Danlos Syndrome. Animals-Basel 10, 2002. https://doi.org/10.3390/ani10112002.

83. Timin, G., and Milinkovitch, M.C. (2023). High-resolution confocal and light-sheet imaging of collagen 3D network architecture in very large samples. *iScience* 26, 106452. https://doi.org/10.1016/j.isci.2023.106452.

84. Arzt, M., Deschamps, J., Schmied, C., Pietzsch, T., Schmidt, D., Tomancak, P., Haase, R., and Jug, F. (2022). LABKIT: Labeling and Segmentation Toolkit for Big Image Data. Front. Comput. Sci. 4, https://doi.org/10.3389/fcomp.2022.777728.

85. Lösel, P.D., van de Kamp, T., Jayme, A., Ershov, A., Faragó, T., Pichler, O., Jerome, N.T., Aadepu, N., Bremer, S., Chilingaryan, S.A., et al. (2020). Introducing Biomedisa as an open-source online platform for biomedical image segmentation. Nat. Commun. 11, 5577. https://doi.org/10.1038/s41467-020-19303-w.

86. Alliez, P., and Fabri, A. (2009). Computational Geometry Algorithms Library. ACM SIGGRAPH Asia 2009 Courses (Association for Computing Machinery).

87. Schönberger, J.L., Zheng, E., Frahm, J.-M., and Pollefeys, M. (2016). Pixelwise View Selection for Unstructured Multi-View Stereo. Held in Cham (Springer International Publishing), pp. 501-518.

88. Ben Amar, M., and Jia, F. (2013). Anisotropic growth shapes intestinal tissues during embryogenesis. Proc. Natl. Acad. Sci. USA 110, 10525-10530. https://doi.org/10.1073/pnas.1217391110.

89. Cignoni, P., Callieri, M., Corsini, M., Dellepiane, M., Ganovelli, F., and Ranzuglia, G. (2008). MeshLab: an Open-Source Mesh Processing Tool. Sixth Eurographics Italian Chapter Conference, 129-136.
<figure><img src="image_17.png" /></figure>

## STAR★METHODS

### KEY RESOURCES TABLE

<table><thead><tr><th>REAGENT or RESOURCE</th><th>SOURCE</th><th>IDENTIFIER</th></tr></thead><tbody><tr><td colspan="3"><strong>Biological samples</strong></td></tr><tr><td>Dogs (embryos and juveniles) of the Beagle breed</td><td>École Nationale Vétérinaire d'Alfort, (facility 947-046-2, ENVA, Maisons-Alfort, France)</td><td>N/A</td></tr><tr><td>Ferrets Mustela putorius furo (embryos, juveniles and adults)</td><td>Euroferret (Denmark) and Animal Facilities of the Universidad Miguel Hernández (Alicante, Spain)</td><td>N/A</td></tr><tr><td>Bovine embryos</td><td>INRAE experimental unit of Le-Pin-au-Haras (Gouffern-en-Auge, France)</td><td>N/A</td></tr><tr><td>Bovine embryos</td><td>Abattoir SOCOPA (Gacé, France)</td><td>N/A</td></tr><tr><td>Bovine embryos</td><td>Abattoir des Marais (Estavayer-le-Lac, Switzerland)</td><td>N/A</td></tr><tr><td>Muzzles from 15 adult Holstein cows (9 genetically identical clones, 4 unrelated controls, and 2 paternal half-sisters affected by Ehlers-Danlos syndrome)</td><td>ENVA and INRAE experimental unit of Bressonvilliers (Leudeville, France)</td><td>N/A</td></tr><tr><td colspan="3"><strong>Chemicals</strong></td></tr><tr><td>Acrylamide/Bis-acrylamide (37.5 :1 solution)</td><td>Sigma-Aldrich</td><td>#A3699</td></tr><tr><td>Fast Green FCF</td><td>Thermo Scientific</td><td>#A16520.22</td></tr><tr><td>TOPRO3</td><td>ThermoFisher</td><td>#T3605</td></tr><tr><td colspan="3"><strong>Antibodies</strong></td></tr><tr><td>Phosphorylated Histone H3 (pH3)</td><td>Santa Cruz Biotechnology</td><td>Cat#sc-374669; RRID: AB_11150094</td></tr><tr><td colspan="3"><strong>Software and algorithms</strong></td></tr><tr><td>ImageJ</td><td>https://imagej.net/ij/</td><td>N/A</td></tr><tr><td>Matlab</td><td>https://www.mathworks.com/products/matlab.html</td><td>N/A</td></tr><tr><td>Meshlab</td><td>https://www meshlab.net/</td><td>N/A</td></tr><tr><td>Amira</td><td>https://www.thermofisher.com/ch/en/home/industrial/ electron-microscopy/electron-microscopy-instruments- workflow-solutions/3d-visualization-analysis-software/ amira-life-sciences-biomedical.html</td><td>N/A</td></tr><tr><td>Imaris</td><td>https://imaris.oxinst.com/imaris-viewer</td><td>N/A</td></tr><tr><td>Biomedisa</td><td>https://biomedisa.info/</td><td>N/A</td></tr><tr><td>CGAL</td><td>https://www.cgal.org/</td><td>N/A</td></tr><tr><td>COLMAP</td><td>https://colmap.github.io/</td><td>N/A</td></tr><tr><td>Dear ImGui</td><td>https://github.com/ocornut/imgui</td><td>N/A</td></tr><tr><td>CUDA® C++ Core Libraries</td><td>https://github.com/NVIDIA/cccl</td><td>N/A</td></tr><tr><td>Eigen</td><td>https://gitlab.com/libeigen/eigen</td><td>N/A</td></tr><tr><td>libigl</td><td>https://github.com/libigl/libigl</td><td>N/A</td></tr><tr><td colspan="3"><strong>Deposited data</strong></td></tr><tr><td>Source code for the growth simulations</td><td>This paper</td><td>https://codeocean.com/capsule/ 6939939/tree/v1; DOI: https://doi.org/ 10.24433/CO.5622496.v1</td></tr></tbody></table>

## EXPERIMENTAL MODEL AND SUBJECT DETAILS

Sampling of dog embryos were carried out in accordance with animal welfare ethical principles and approved by the ethics committees of the French national review board CNREEA #16. The procreation of bovine embryos was approved by the French Ministry of Higher Education, Research and Innovation (permit number APAFIS #32140-2021062809539809). The bovine clones were generated as part of another study approved by the local Institutional Animal Care and Use Committee of AgroParisTech/INRAE (COMETHEA, permit number 2012/160). No ethical approval was required for the sampling of Ehler-Danlos affected cows, control cows and embryos obtained opportunistically at the Abattoir des Marais, as they were not deliberately bred for this study and since
sampling was performed after euthanasia or slaughter of cows not known to be pregnant. Ferrets were treated according to Spanish (RD 53/2013) and European regulations, and experimental protocols were approved by the CSIC (Consejo Superior de Investigaciones Científicas) Ethics Committee. All experiments reported in this article comply with the ethical guidelines of the research institutions involved in this study and with the EU guidelines and regulations (Directive 2010/63/UE).

## METHOD DETAILS

### Biological samples

All samples were collected post-mortem. Dog foetuses of the Beagle breed were opportunistically collected at the École Nationale Vétérinaire d'Alfort, (facility 947-046-2, ENVA, Maisons-Alfort, France). In dogs, embryonic stages are counted from the day of ovulation (E0) which was precisely monitored by ultrasound and progesterone assays. Fertilization and beginning of embryonic development in dogs occurs 3 days after ovulation. The samples were fixed and stored in formalin. Timed-pregnant sable ferrets (*Mustela putorius furo*) were obtained from Euroferret (Denmark) and maintained on a 16-hour light/8-hour dark cycle at the Animal Facilities of the Universidad Miguel Hernández (Alicante, Spain). The day after mating was counted as embryonic day E0, and the day of birth counted as postnatal day P0. Embryos and kits were transcardially perfused with 4% paraformaldehyde and stored. Fixed samples were stored in methanol, ethanol or formalin. Some of the bovine embryos were produced by artificial insemination of cull cows after oestrus synchronization at the INRAE experimental unit of Le-Pin-au-Haras (Gouffern-en-Auge, France) and collected at the Abattoir SOCOPA (Gacé, France). In bovine samples, embryonic stages are given in days after insemination (E0), which occurred a few hours after the observation of oestrus. The remainder were opportunistically obtained from the abattoir des Marais (Estavayer-le-Lac, Switzerland). The muzzles from 15 adult Holstein cows were collected at the ENVA and the INRAE experimental unit of Bressonvilliers (Leudeville, France). They were fixed in formalin and stored in this medium, which was changed regularly, until the samples were used for analysis. The specimens included nine genetically identical clones, four unrelated controls, and two paternal half-sisters affected by Ehlers-Danlos syndrome — a disease caused by structurally abnormal dermal collagen fibrils<sup>81,82</sup> — which were found to be heterozygous for a *de-novo* mutation truncating the collagen type V alpha 1 chain (COL5A1 p.E1565X).

### Histology

2D sections of embryonic skin pieces or complete rhinaria (small enough to fit into a regular microtome, typically less than 2 cm³) were prepared using a standard histology protocol (fixed in 4% PFA, dehydrated in ethanol and delipidated in chloroform before embedding in paraffin; sectioned with a microtome at 9-µm thickness, stained with Hematoxylin and Eosin) to reveal the geometry of the different skin layers.

### Light-sheet microscopy

We used a light-sheet microscope (Ultramicroscope from Miltenyi Biotec) to perform 3D imaging of cell densities in full rhinaria of embryos of all three species, as well as in noses from juvenile dogs, juvenile ferrets and adult ferrets. Whole-mount samples (maximum 2 cm³) were prepared using a clearing protocol for volumetric imaging, consisting of the following successive steps: dehydration (methanol), bleaching (hydrogen peroxide), permeabilization (saponin + triton) and blocking (bovin gelatin), Phosphorylated Histone H3 (pH3) and TOPRO3 staining (to mark proliferating cells and nuclear double-stand DNA, respectively), delipidation in dichloromethane, and clearing in dibenzyl ether. FastGreen staining was used to mark collagen fibers. Additional details on the protocol can be found in Timin et al.<sup>83</sup> Imaging was done at 4X magnification (for tissue 3D segmentation) or 12X (for segmentation of cells).

### Hydrogel analogues

Polyacrylamide solutions (10%) were prepared following Tanaka et al.,<sup>16</sup> with a ratio of cross-linker to polymer (bis-acrylamide to acrylamide) of 1:37.5. Gels were polymerized in glass containers (hemispherical molds or Petri dishes) covered with a glass plate to prevent oxygen from penetrating the surface and weakening the cross-linking process. After complete polymerization, the gels were removed from the molds and let to dry for several hours, thus shrinking to a minimal size. The hemispheres were super-glued face down on a rigid surface. Finally, the dried gels were immersed in water, causing the outer layer to swell gradually. This generates a bilayer system where the growth of the outer layer, constrained mechanically by the static core (and by the glued face or the borders of the Petri dish), induces the spontaneous formation of creasing patterns of polygonal units. The gels (Figure S1) were imaged using a VHX-700 digital microscope (Keyence, Itasca, IL, USA) and with a M205 fluorescence microscope (Leica, Heerbrugg, Switzerland), relying on the autofluorescence of polyacrylamide.

## QUANTIFICATION AND STATISTICAL ANALYSIS

The software Imaris (https://imaris.oxinst.com/imaris-viewer) was used to render plane sections or maximal intensity projections (MIPs) of the 3D images. The spot detection algorithm from that software was employed to count the number of cells per unit volume in different skin layers, e.g., dermis versus epidermis (either all cells marked by TOPRO3 or only dividing cells marked by pH3), as shown in Video S9.
### 3D segmentation of tissues

We segmented volumetric light-sheet images using a combination of the ImageJ plugin Labkit,⁸⁴ the online platform Biomedisa,⁸⁵ and manual tools from the segmentation editor of the software Amira (ThermoFisher Scientific, Massachusetts, USA). Labkit relies on a 'supervised random forest' method, which requires orders of magnitude less ground-truth data than a typical deep neural network. The training data set was created by drawing over sections of the 3D image, to label small zones for each segmentation class (here, dermis, epidermis, basal layer, blood vessels and cartilage). After training, the classifier is applied to all pixels to predict their labels, with an interactive curation step to correct mistakes until a satisfactory classification is obtained. The results were cleaned manually with the segmentation editor from Amira to remove any remaining inaccuracies. The Biomedisa online platform comprises a 3D U-Net, i.e., a convolution neural network (CNN) for pixel classification. We used this CNN to segment the most complex tissue layers (blood vessels and cartilage), which were not well reconstructed by the Labkit method. The training data set for the CNN was created by applying Labkit over a subset of the volumetric image (1/20 of all 2D XY slices), including a manual cleaning step with Amira. In the end, label fields of all tissue layers were combined into a single volume, and tetrahedral meshes were generated for each layer using the Computational Geometry Algorithms Library (CGAL)⁸⁶ with adaptive cell size, i.e., smaller tetrahedra in the thinner layers (epidermis, basal layer, blood vessels) and larger units in the bulk layers (dermis and cartilage). The result is shown in panel (ii) of Figure 7A.

### 3D surface reconstruction of large samples

The skin surfaces of adult cow noses were reconstructed as 3D triangular meshes using a structure-from-motion approach (also known as 'multiview photometric stereo scanning'). To do so, high-resolution pictures of the surfaces were taken with a Nikon D810 camera from a multitude of angles by a custom-built robotic system⁶³ and combined to recreate the 3D geometry using the COLMAP software.⁸⁷ Alternatively, for intermediate-sized samples (juvenile dogs), a structured-light method was used to reconstruct the surface, with a VR-6000 optical profilometer (Keyence, Itasca, IL, USA).

### Differential growth numerical simulations

We performed growth simulations using tetrahedral volumetric meshes (embryonic dog nose and flat-slab multilayered models), starting with a smooth-surface configuration. The relative thicknesses of the flat-slab models were chosen to reproduce the relative thickness of the basal layer compared to the upper epidermis in the true rhinarium geometry. The tangential growth factors, Poisson's ratios and Young's moduli of individual layers were varied to test the effect of each parameter on the final surface pattern. At each simulated time step, elastic forces induced by the deformation of the growing tetrahedral elements are calculated based on the minimization of the strain energy function. The skin is modeled as a soft neo-Hookean material. The material model and energy minimization process are formulated in the finite element framework with tetrahedral elements.

We define the bulk of the material as the collection of points forming a specific configuration in 3D space. The motion of these points fully describes its deformation. We define the deformation gradient, $F$, as a second-order tensor which relates the initial (reference) configuration of any material point $X$ to its current (deformed) configuration $x$:

$$F = \frac{dx}{dX} \qquad (Equation 1)$$

To model tissue growth, the deformation gradient tensor $F$ is decomposed into elastic and growth components $F_e$ and $F_g$, respectively:

$$F = F_e F_g, \qquad (Equation 2)$$

Here, $F_g$ describes the effective growth of the tissue and can be derived from the biological experiments/observations. In our model, we assume that the growth deformation gradient is saturating as a function of time:

$$F_g(t) = I + \lambda_N(1 - \exp(-\beta t))N \otimes N + \lambda_T(1 - \exp(-\beta t))(I - N \otimes N) \qquad (Equation 3)$$

in which, $N$ is the normal vector in the reference configuration, $\beta > 0$ is a constant that determines growth rate, and $I$ is the identity matrix. Furthermore, $\lambda_N$ and $\lambda_T$, determine the final relative growth in the normal and tangential directions, respectively.

Neo-Hookean material models (a sub-group of hyper-elastic models) are appropriate for soft material (bio-)mechanics because they are robust even under large deformations.³²,³⁶,³⁸ In hyper-elastic models, the stress tensor is computed from the strain energy density function:

$$\psi = \frac{\mu}{2} (\mathrm{tr}(\mathbf{F}_e \mathbf{F}_e^T) J^{-\frac{2}{3}} - 3) + K(J - \mathrm{ln}J - 1), \qquad (Equation 4)$$

where $J = \det F_e$ and $\mu$ and $K$ are shear and bulk moduli, respectively, which are related to Young's modulus $E$ and Poisson's ratio $\nu$ as follows:

$$\mu = \frac{E}{2(1+\nu)} \quad \text{and} \quad K = \frac{E}{3(1-2\nu)}. \qquad (Equation 5)$$
In the framework of hyper-elastic materials, we compute the Cauchy stress $\sigma$ directly from the strain energy density function:

$$ \sigma = \frac{1}{J} \frac{\partial \Psi}{\partial \mathbf{F}_e} \mathbf{F}_e^T \qquad (\text{Equation 6}) $$

Substituting (Equation 4) into (Equation 6) results in the following stress function

$$ \sigma = \mu J^{-5/3} \mathrm{dev}(\mathbf{F}_e \mathbf{F}_e^T) + K \left(1 - \frac{1}{J}\right) I \qquad (\text{Equation 7}) $$

Stress is computed for all tetrahedra, each of which is associated to well-defined material properties (e.g., Young's modulus and Poisson's ratio) and growth deformation.

Internal forces are computed from this stress term and used for inferring mesh node displacements with Newton's second law of motion. A low threshold ratio (=10⁻⁵) of kinetic energy over total energy is used to check for convergence of the steady-state solution.

All formulations are implemented in an in-house application which exploits NVIDIA GPUs for higher efficiency. This application is written in CUDA and C++ programming languages. CUDA is used to develop intensive-computation kernels whereas C++ is used for the data management, geometry processing, input/output operations and the graphical user interface. Our application integrate the following open-source libraries: Dear ImGui (https://github.com/ocornut/imgui, MIT licence) for graphical user interface, CUDA® C++ Core Libraries (https://github.com/NVIDIA/cccl, Apache-2.0, FreeBSD, BSD-3-Clause licences) for parallel algorithms, Eigen (https://gitlab.com/lib eigen/eigen, MPL-2.0, BSD licences) for linear algebra, and libigl (https://github.com/libigl/libigl, GPL-3.0, MPL-2.0 licences) for geometry processing. Source code and examples of meshes are provided in the form of a CodeOcean capsule at https://codeocean.com/capsule/6939939/tree/v1; the capsule's DOI is https://doi.org/10.24433/CO.5622496.v1.

## Statistical analysis of polygonal networks

We computed with Meshlab⁸⁹ the maximum curvature of the 3D surface meshes obtained from juvenile and adult rhinaria. The vertices of the 3D meshes were then treated as point clouds and projected onto regular volumetric grids, and the maximum curvature values were binarized with Matlab (Version R2021a) using an ad hoc threshold to skeletonize the networks of creases in 3D. Various attributes of the polygonal domes were then measured, to create a statistical description of the pattern: (i) the angular orientation of each unit relative to the global orientation of the surface (calculated as the angle between the principal orientation vector of the whole rhinarial surface, obtained by principal component analysis of the point cloud formed by the mesh vertices of the whole surface, and the principal orientation vector of each individual epidermal unit, also determined by PCA of the point cloud formed by the vertices of that particular unit); (ii) the mean height of each unit, measured from a base form approximated by fitting a polynomial surface of degree 4 to the skeletonized folds network (the folds being thus considered as the lowest points on the skin), inside a local patch circle with an ad-hoc radius; (iii) the surface area of each dome (projected on its base form); (iv) the length of each dome across its two principal directions and its corresponding aspect ratio; (v) the number of branch nodes (corners of the polygons, connecting different branches in the folds network); and (vi) the number of end nodes (unconnected edge node) in each unit. This process generated a statistical distribution for each geometrical attribute and for each nose sample. In order to compare the distributions between the different rhinaria, we used the Jensen-Shannon divergence ($D_{JS}$). It measures the level of discrepancy between two statistical distributions, defined as:

$$ D_{JS}(P \parallel Q) = \frac{1}{2} D_{KL}(P \parallel M) + \frac{1}{2} D_{KL}(Q \parallel M), \qquad (\text{Equation 8}) $$

where $D_{KL}$ is the Kullback-Leibler divergence:

$$ D_{KL}(P \parallel Q) = \sum_{x \in X} P(x) \log_2 \frac{P(x)}{Q(x)} \qquad (\text{Equation 9}) $$

where P and Q are two discrete probability distributions defined over the domain X, and M is the mean of the two distributions. The values of $D_{JS}$ are symmetric and bounded between 0 (identical distributions) and 1 (maximally divergent distributions). $D_{JS}$ values were calculated pair by pair for the 19 nose samples (9 cow clones, 4 cow controls, 2 cows with Ehlers-Danlos, 2 dogs, 2 ferrets), for every geometrical attribute. Finally, those $D_{JS}$ values were averaged over all samples belonging to the same group (cow clones, cow controls, cows with E-D, dogs, ferrets) and over all geometrical attributes, yielding a single parameter value between 0 and 1 to define the level of statistical divergence between the rhinarial patterns of the different groups. This quantification is illustrated in Figure S6.