# Predictive model for cytoneme guidance in Hedgehog signaling based on Ihog- Glypicans interaction

Adrian Aguirre-Tamaral, Manuel Cambon, David Poyato, Juan Soler, Isabel Guerrero

► To cite this version:

Adrian Aguirre-Tamaral, Manuel Cambon, David Poyato, Juan Soler, Isabel Guerrero. Predictive model for cytoneme guidance in Hedgehog signaling based on Ihog- Glypicans interaction. Nature Communications, 2022, 13 (1), pp.5647. 10.1038/s41467-022-33262-4. hal-04134082

HAL Id: hal-04134082

https://hal.science/hal-04134082v1

Submitted on 3 Sep 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Predictive model for cytoneme guidance in Hedgehog signaling based on Ihog- Glypicans interaction

Received: 11 November 2021

Accepted: 9 September 2022

Published online: 26 September 2022

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

Check for updates

Adrián Aguirre-Tamaral ⒰<fnref n="1" />,², Manuel Cambón ⒰², David Poyato ⒰²,³,
Juan Soler ⒰²✉ & Isabel Guerrero ⒰<fnref n="1" />✉

During embryonic development, cell-cell communication is crucial to coor-
dinate cell behavior, especially in the generation of differentiation patterns via
morphogen gradients. Morphogens are signaling molecules secreted by a
source of cells that elicit concentration-dependent responses in target cells.
For several morphogens, cell-cell contact via filopodia-like-structures (cyto-
nemes) has been proposed as a mechanism for their gradient formation.
Despite of the advances on cytoneme signaling, little is known about how
cytonemes navigate through the extracellular matrix and how they orient to
find their target. For the Hedgehog (Hh) signaling pathway in *Drosophila*, Hh
co-receptor and adhesion protein Interference hedgehog ( Ihog) and the gly-
picans Dally and Dally-like-protein (Dlp) interact affecting the cytoneme
behavior. Here, we describe that differences in the cytoneme stabilization and
orientation depend on the relative levels of Ihog and glypicans, suggesting a
mechanism for cytoneme guidance. Furthermore, we have developed a
mathematical model to study and corroborate this cytoneme guiding
mechanism.

Cell-cell communication is crucial during organism development, and it is mainly mediated through specific signaling molecules called morphogens. These signals are distributed in a graded form within a morphogenetic field, activating different target genes in a concentration-dependent manner<fnref n="1" />. Classical models assume that the graded distribution of signals occurs through simple diffusion², although for many morphogens their biochemical properties inherently impede their free diffusion through the extracellular matrix (ECM)³. For those cases, experimental evidence suggests a mechanism for morphogen transport based on its distribution by specialized signaling filopodia called cytonemes⁴⁻⁷.

Cytonemes play a critical role in the process of cell communica-
tion in several biological systems, both during development and in
tissue homeostasis in adults⁸,⁹. Thus, understanding the underlying

mechanisms of cytoneme behavior turns essential. Here, we focus on
the study of Hedgehog (Hh) signaling cytonemes in the *Drosophila*
wing imaginal disc. This epithelium is divided into two
compartments:¹⁰ the posterior (P) compartment cells produce the Hh
signal that is received by the anterior (A) compartment cells. In this
system, cytonemes protrude from both producing and receiving cells,
driving both short- and long-distance cell signaling¹¹,¹². The interaction
of cytonemes from both compartments at specific contact sites, where
all the components of the reception complex interact, promotes Hh
reception¹¹.

Previous studies have identified several proteins involved simul-
taneously in Hh signaling and cytoneme dynamics¹¹. A key protein is
the Hh co-receptor Interference hedgehog (Ihog), an adhesion protein
containing four immunoglobulin (IG) domains, two fibronectins type

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

<footnote n="1">¹Tissue and Organ Homeostasis, Centro de Biología Molecular "Severo Ochoa" (CSIC-UAM), Nicolás Cabrera 1, Universidad Autónoma de Madrid, Cantoblanco, E-28049 Madrid, Spain. ²Departamento de Matemática Aplicada and Research Unit "Modeling Nature" (MNat), Facultad de Ciencias, Universidad de Granada, E-18071 Granada, Spain. ³Institut Camille Jordan (ICJ), UMR 5208 CNRS & Université Claude Bernard Lyon 1, F-69100 Villeurbanne, France.</footnote>
III (FNIII) domains, a single-pass transmembrane domain and an
intracellular domain<sup>13,14</sup>. The overexpression of Ihog is able to stabilize
(freeze) cytoneme dynamics. This stabilization has been described and
quantified in<sup>11,12,15</sup>, where cytonemes seem to be static.

Other proteins required for long-range Hh gradient formation and
also present in cytonemes are the ECM components Heparan sulphate
proteoglycans (HSPGs)<sup>16,17</sup>. In *Drosophila*, two glypicans belonging to
the HSPGs family have been identified: Division abnormally delayed
(Dally) and Dally-like protein (Dlp). Dally and Dlp also intervene in
cytoneme dynamics and their interaction with Ihog facilitates cyto-
neme mediated cell contact<sup>11</sup>. The interaction of Ihog with Dally and
Dlp occurs through the first fibronectin Fn1 domain of Ihog<sup>15,18</sup>. Thus,
the Fn1 domain is involved in both Ihog-Hh and Ihog-glypicans het-
erophilic interactions, Ihog-Ihog homophilic interaction for cell-cell
adhesion, and in the regulation of cytoneme dynamics. The double
function of Ihog and glypicans in controlling both signaling reception
and cytoneme behavior suggests that coordination between Hh
reception, cytoneme behavior and their spatial orientation is essential
for efficient signaling.

The dynamics of Hh cytonemes over time have already been
studied<sup>11,15,19</sup> together with the correlation of dynamic cytonemes with
the morphogen signaling gradient<sup>1,2,0</sup>. Here, we complement previous
studies exploring how cytonemes are oriented and guided for the
correct Hh signaling. In this work we observe that the interplay
between different levels of Ihog and glypicans can orient and guide
cytonemes in the Hh producing and receiving cells. To better under-
stand this interplay, we develop a theoretical model that combines
experimental data with mathematical tools and physical principles.
Our model is able to predict cytoneme trajectories under different
levels of Ihog and glypicans; these predictions have been experimen-
tally validated in *Drosophila* wing imaginal discs. We propose that the
glypicans-Ihog interaction can provide the spatial information
required for cytoneme guidance during Hh reception.

Results

Ihog overexpression reveals cytoneme spatial behavior

Both the thinnness and highly dynamic behavior of cytonemes hamper
their observation in fixed tissues<sup>12,21</sup>, so it is important to mention that
the non-visualization of cytonemes in fixed tissues in not equivalent to
the absence of cytonemes. Previous in vivo experiments in abdominal
histoblast nest, where the visualization of temporal dynamic cyto-
nemes is feasible in physiological conditions<sup>12</sup>, have shown that the
stabilization of cytonemes by overexpression of Ihog, in the same cells
(overexpression in *cis*), is around 70% of the total cytoneme popula-
tion (Supplementary Fig. 1a<sup>11</sup>) without changing their length and
orientation (Supplementary Fig. 1b<sup>11</sup>). By stabilization we mean that
cytonemes remain static ("frozen") for up to 9 hours<sup>12</sup>, while their
typical wild type lifetime is about 10 min<sup>11</sup>. We observed the same
stabilization effect in ex vivo experiment in the wing imaginal disc
(Supplementary Movie 1). This change in the dynamics of cytonemes,
making them static, has been previously referred to as cytoneme sta-
bilization by Ihog overexpression<sup>11,12,15</sup>, probably forming bundles of
two or more cytonemes<sup>18</sup>. This condition has been used as a genetic
tool to visualize and study cytonemes in previous works<sup>11,12,18,22</sup> and we
have used it here to study their spatial behavior.

In the wing imaginal disc, Ihog overexpression shows that cyto-
nemes from both Hh producing and receiving cells protrude perpen-
dicularly to the A/P compartment border (Fig. 1a<sup>12</sup>). This alignment is
also present in wild type<sup>12</sup> (Supplementary Fig. 2), probably serving to
facilitate the encounter of cytonemes protruding from the P com-
partment cells, carrying the Hh signal, with cytonemes protruding
from the A compartment cells receiving it<sup>11</sup>.

To experimentally dissect the orientation of cytonemes, we first
confronted the P compartment cell population, overexpressing Ihog,
to A compartment cell populations expressing different levels of Ihog.

We will refer as *trans* interaction to the effect on the behavior of
cytonemes protruding from a cell population when it is confronted to
a genetically different cell population. Cytonemes protruding from P
compartment cells overexpressing Ihog are stabilized when con-
fronted to A compartment cells in which Ihog levels are wild-type
(Fig. 1a<sup>12</sup>) or downregulated (Fig. 1b<sup>11</sup>), while no cytonemes are
observed when both confronted A and P cell populations equally
overexpress Ihog (Fig. 1c). We quantified and compared the stabiliza-
tion effect and orientation of cytonemes under the above conditions
(Fig. 1d-f). To ascertain that the apparent lack of cytonemes in the
confronted condition is due to non-stabilization instead of absence, we
performed in vivo experiments in abdominal histoblast nets, and
shown in Supplementary Movie 2, dynamic cytonemes are visible.

To quantitatively explore the spatial effect of Ihog on cytoneme
behavior, we induced clones of different sizes and positions using a
genetic tool that allows the random generation of two different types
of confronted clones. In these clones, fluorescent labeling, as well as
changes in Ihog levels, can be independently manipulated using the
LexA and Gal4 drivers (Fig. 2a). These experiments showed that
proximity is important for cytoneme stabilization (Fig. 2b). Cytoneme
length is affected by the proximity of another clone if it is close enough
(<15 microns Fig. 2c). This effect becomes statistically stronger the
shorter the distance (Fig. 2c); actually, cytonemes in close proximity to
a clone can only be stabilized if they change their orientation to cir-
cumvent the confronted clone (Fig. 2d).

We have also studied the dependence of cytoneme behavior on
Ihog levels by confronting cells overexpressing Ihog in the whole P
compartment with groups of cells expressing different levels of Ihog
(ectopic cell clones) within the A compartment. If an innocuous label
such as the actin reporter LifeActGFP<sup>23</sup> is used, cytonemes stabilized by
Ihog overexpression in the P compartment are seen (Supplementary
Fig. 3a), while they are not seen when the confronted clones express
high levels of Ihog (Supplementary Fig. 3b yellow arrows). Further-
more, when Ihog overexpressing cells are confronted with clones
ectopically expressing Ihog RNAi, thereby decreasing the levels of Ihog
in *trans*, cytonemes are visualized oriented towards the A receiving
cells (Supplementary Fig. 3c red and yellow arrows). These results
indicate that differences in Ihog levels between cytoneme membranes
in *trans* play a decisive role in cytoneme dynamics.

Based on these results, it could be proposed that the wild-type
distribution of Ihog in the wing disc, strongly downregulated in the Hh
receiving cells<sup>24</sup> (Supplementary Movie 3), could suffice to initiate
cytoneme orientation from P to A cells; that is, from regions with high
levels of Ihog in the P compartment towards regions with lower levels
in the Hh receiving area of the A compartment (see Supplementary
Movie 3). This hypothesis, solely based on the wild-type levels of Ihog,
cannot be valid for those cytonemes that orient from Hh receiving cells
towards Hh producing P cells, with higher levels of Ihog. Reasonably,
other factors might also be involved in achieving cytoneme orienta-
tion. Although Hh could have been a good candidate, at the moment
we lack convincing evidences of a possible role of Hh in cytoneme
generation or dynamics<sup>11</sup>.

Glypicans influence in *trans* the behavior of cytonemes expres-
sing ectopic Ihog

Cytoneme stabilization also depends on the heterophilic interaction
between Ihog and glypicans<sup>11,15,18</sup>; such stabilization in wing disc cells
requires the presence in *trans* (neighboring cells) of the HSPGs<sup>12</sup> and,
more specifically, of glypicans<sup>11</sup>. To test whether glypicans might also
be required within the same cells (interaction in *cis*), we induced *tout
vellue* (*ttv*) and *brother of ttv* (*botv*) double mutant clones while
expressing also Ihog in the P compartment (Fig. 3a). *ttv* and *botv* are
genes coding for enzymes that synthesize the Heparan Sulfate chains
of the glypicans Dally and Dlp, so the effect of their absence on Hh
signalling is equivalent to that of the lack of glypicans<sup>16,17</sup>. Cytonemes
<figure><img src="image_5.png" /><figcaption>**Fig. 1 | Cytoneme stabilization by Ihog overexpression.** a P compartment cells overexpressing (↑) Ihog in *cis* (grey) are able to stabilize cytonemes over wild type (WT) A compartment cells (in *trans*). b P compartment cells overexpressing (↑) Ihog (grey) in *cis* are able to stabilize cytonemes over the A compartment cells (in *trans*) where Ihog is downregulated (↓). c P compartment cells overexpressing (↑) Ihog (grey) do not stabilize over the A compartment cells where Ihog is also overexpressed (↑, green) in *trans*. The schemes at the right depict the cytoneme behavior in the different genetic conditions. d Summary table depicting how cytoneme stabilization depends on Ihog levels in *trans*. e Boxplot and statistical</figcaption></figure>

study was performed using a Wilcoxon rank sum test per pairs over the number of P cytonemes per disc in different experimental conditions. f Orientation of P cytonemes in different conditions shows that P cytonemes are perpendicular to A/P compartment border. Data was quantified from *n* = 23 discs for section E and from 513 cytonemes of *n* = 16 discs for section F. In box plots, a box indicates the median (in red and 25 and 75 percentiles, whiskers indicate range of data and crosses indicate outliers. Raw data and p-values (ns; *p* > 0.05; * *p* ≤ 0.05; ** *p* ≤ 0.01; *** *p* ≤ 0.001; **** *p* ≤ 0.0001) are provided as a Source data file. Scale bars: 15 µm.

overexpressing Ihog emanating from HSPG-deficient cells (clone 2 in Fig. 3b, yellow arrowhead) are stabilized when navigating through wild-type territory, in contrast to cytonemes emanating from wild-type cells that are not stabilized when traversing glypican-deficient cells (clone 1 in Fig. 3b, red arrowhead)¹². These data demonstrate that glypicans have an effect on Ihog-induced stabilization of cytonemes in *trans* but not in *cis*, most probably by enhancing the amount of Ihog in those cytonemes that go across glypican-deficient clones.

Accordingly, cytonemes are stabilized, probably in bundles, and oriented when a cell population with high levels of Ihog is in the proximity of another cell population (*trans*) with high levels of either Dally (Fig. 3c) or Dlp (Fig. 3d). However, cytonemes are not stabilized in

wild type regions if cytonemes come from a region co-expressing Ihog and Dally or Ihog and Dlp in the same cells (*cis*) (Fig. 3c', d'). These results suggest that cytoneme stabilization may be due to molecular competition dependent on the availability of Ihog and glypicans.

In this context, we observed that stabilization is recovered when a co-expressing region (Ihog + Dally or Dlp) is confronted to a region where either Dally or Dlp are overexpressed (Supplementary Fig. 4, red arrows). Interestingly, in those conditions we can observe how cytonemes are oriented to reach regions with high levels of Ihog and glypicans (Supplementary Fig. 4, yellows arrows). Therefore, there appears to be an orientation effect mediated by the Ihog-glypicans interaction in *trans* providing directionality to cytonemes.
<figure><img src="image_6.png" /></figure>

**Fig. 2** | Cytoneme behavior depends on the distance to the Ihog source. a Ihog overexpression clones marked with different fluorophores. Blue arrows point to cytonemes stabilized towards areas without clones; red arrows point to touching clones showing no cytoneme stabilization when confronted; green arrows point to stabilized cytonemes crossing small clones (<2 cell diameters); <18% of the touching clones presented cytoneme stabilization (yellow arrows) for large clones (>2 cell diameters). b Cytoneme length (y-axis) versus distance between clones (x-axis) in µm. The red dotted line indicates the distances where the data cloud has a different behavior (<15 µm); the green dotted line represents the linear relation between the length of stabilized cytoneme and the distance between clones at short distances. c Statistical studies was performed using a Wilcoxon rank sum test

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

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

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

per pairs over different regions (boxplots of 5 µm) showed a significant change of cytoneme behavior if the clones were closer than 15 µm. d We studied in detail cytonemes over the green dotted line (i.e., stabilized cytonemes which are longer than the distance with the closest clone) and we observed that those cytonemes able to stabilize show a change in their orientation (α) to avoid the confronted clones (α = 0 means no change in orientation and α = 90 a change towards opposite direction). Data were quantified along 235 cytonemes from different clones (n = 21). In box plots, a box indicates the median (in red) and 25 and 75 percentiles, whiskers indicate range of data and crosses indicate outliers. Raw data and *p*-values (ns; *p* > 0.05; *p* ≤ 0.05; ***p* ≤ 0.01; **p* ≤ 0.001; ****p* ≤ 0.0001) are provided as a Source data file. Scale bars: 15 µm.

## The Ihog Fn1 domain is essential for the *trans* interactions with Ihog and glypicans

Previous works have demonstrated that the cell adhesion of both the Ihog-Ihog homophilic and the Ihog-glypicans heterophilic interactions occur through the first Ihog fibronectin domain (Fn1)¹⁵,¹⁸,²², the same domain involved in heparin- dependent Hh binding²². Therefore, to analyze the Ihog interactions leading to cytoneme orientation, we generated an Ihog form (LexAop.ihogΔFn1-RFP) deficient for the Fn1 domain, fused to RFP and under the control of the Lex inducible promoter. Supporting our previous results, we detected that the stability of cytoneme bundles observed in the *trans* Ihog-glypicans interaction (Fig. 4a¹¹) is lost when the IhogΔFn1 is expressed (Fig. 4b). Furthermore, in contrast to the absence of Ihog cytoneme stabilization, observed after confronting two cell populations both over-expressing wild type Ihog (Fig. 1c), cytonemes can now be stabilized when confronting a cell population overexpressing wild type Ihog to a

cell population overexpressing the Ihog Fn1 deficient form (Fig. 4c). These changes in cytoneme behavior are statistically significant (Fig. 4d, e). Altogether, we conclude that the Ihog-glypicans heterophilic interaction might compete with Ihog-Ihog homophilic interaction for cytoneme stabilization and that both interactions require an intact Ihog Fn1 domain.

## Modeling cytoneme guidance: a mathematical framework

To further understand the mechanisms that control cytoneme orientation, we used a mathematical approach aimed at designing an in silico model able to predict cytoneme "guidance" based on the protein levels of Ihog, Dally and Dlp.

The experimental results presented above showed that the effects of those protein levels on the cytoneme spatial behavior could give rise to regions where cytoneme stabilization is favored compared to regions where it is not. These favored/disfavored regions are in some
<figure><img src="image_10.png" /><figcaption>**Fig. 3** | Cytoneme behavior depends on Ihog and glypicans levels. a Scheme depicting induced *ttw*<sup>fl</sup> *botw*<sup>fl</sup> double mutant clones in a wing disc while expressing also Ihog in the P compartment (green). b Wing disc containing induced *ttw*<sup>fl</sup> *botw*<sup>fl</sup> double mutant clones marked by the absence of βgal (b') while expressing Ihog in the P compartment (b''). Cytonemes stabilization can be observed (white arrow) in A compartment cells except in glypicans loss of function clones induced in the A compartment cells, which do not express Ihog (red arrows). The yellow arrow indicates the presence of stabilized cytonemes when the *ttw*<sup>fl</sup> *botw*<sup>fl</sup> double mutant clone is originated in the P compartment, which overexpresses Ihog, indicating that glypicans are required in *trans* and not in *cis* for cytoneme stabilization. c Ihog and</figcaption></figure>

Dally ectopic clones. Strong cytoneme stabilization is observed when Dally clones are confronted in *trans* interaction with P compartment cells overexpressing Ihog. c' Cytonemes are not visualized if Dally clones are in P compartment cells where Ihog is overexpressed (*cis* interaction). d Ihog and Dlp ectopic clones. Strong cytoneme stabilization is observed when Dlp clones are confronted in *trans* interaction with P compartment cells overexpressing Ihog. d' Cytonemes are not visualized when Dlp clones are in the P compartment cells in which Ihog is overexpressed (*cis* interaction). Experimental sample size: *n*<sub>Dally</sub> = 17 and *n*<sub>Dlp</sub> = 16 discs. Scale bars: 15 µm.

sense reminiscent of the classical attractive/repulsive field in physics. With this idea in mind, we used a mathematical framework inspired by the physical theory of electric potentials where the “charges” in physics become “guiding factors” in our biological model. Thus, in our model, the optimal trajectories described by cytonemes across the tissue would be determined by the field lines associated to the “orientation” vector field, in which the levels of Ihog and glypicans are the guiding factors generating the potential of the orientation field.

In what follows, we will briefly describe the building up of the orientation field; further details can be found in Methods and in Supplementary Notes. First, we outline the “charge” distribution associated to the guiding proteins, which serves as the source of the orientation potential that guides cytonemes. These sources will have different behaviors depending on the various levels of Ihog and glypicans described for the wing disc pattern. The spatial dependence of protein distributions in the wing tissue emphasizes the importance of working with the complete 2D information. Accordingly, we define the “charge” distribution $\rho(r, t)$ at the position $r \in \mathbb{R}^2$ of the basal plane of the

imaginal wing disc epithelium as the total concentration of Ihog and glypicans Dlp and Dally in tissue cell and cytoneme membranes, that is,

$$ \rho(r,t) = \beta_{\text{Ihog}} \left( [\text{Ihog}](r) + \sum_{\text{cyt}} [\text{Ihog}]_{\text{cyt}}(r,t) \delta_{\text{cyt}}(r,t) \right) \\ + \beta_{\text{Dlp}} \left( [\text{Dlp}](r) + \sum_{\text{cyt}} [\text{Dlp}]_{\text{cyt}}(r,t) \delta_{\text{cyt}}(r,t) \right) \quad (1) \\ + \beta_{\text{Dally}} \left( [\text{Dally}](r) + \sum_{\text{cyt}} [\text{Dally}]_{\text{cyt}}(r,t) \delta_{\text{cyt}}(r,t) \right). $$

where $[X](r)$ is the density of protein X at the position r of the tissue cell membranes, $[X]_{\text{cyt}}(r,t)$ is the density of protein X at the position r of the cytoneme membranes, that can change spatially and temporally, and $\delta_{\text{cyt}}(r,t)$ is the Dirac function that localizes the density along each cytoneme. “$\sum_{\text{cyt}}$” refers to the sum over all cytonemes protruding from A or P compartment cells. Note that, in particular, $\rho(r, t)$ contains all the collective information of the guiding factors along all cytonemes. $\beta_{\text{Ihog}}$, $\beta_{\text{Dlp}}$, $\beta_{\text{Dally}}$ are proportionality constants that govern the correct guiding factor attractive/repulsive relation among the implicated
<figure><img src="image_11.png" /><figcaption>**Fig. 4** | Ihog and glypicans interaction for cytoneme stabilization in trans is mediated by the Fn1 domain of Ihog. a Cytonemes protruding from both A and P compartments are strongly stabilized when cells overexpressing glypican Dally (green) are confronted to cells overexpressing wild-type Ihog (red). b Cytoneme bundle stabilization is lost when cells overexpressing Dally (green) are confronted to cells overexpressing the Ihog ΔFn1 mutant form (red). The schemes at the right depict the cytoneme behavior in the different genetic conditions. c Cytonemes protruding from cells overexpressing wild-type Ihog (green) are stabilized when confronted to cells overexpressing the Ihog ΔFn1 mutant form (red). d Boxplot and</figcaption></figure>

the statistical study were performed using a Wilcoxon rank sum test per pairs over the number of cytoneme bundles per disc in the different experimental conditions.
e) Boxplot and statistical study of the number of P cytonemes per disc in the different experimental conditions shown in Fig. 1c and panel c. Data were quantified from *n* = 19 discs for section D and from *n* = 25 discs for section E. In box plots, a box indicates the median (in red) and 25 and 75 percentiles, whiskers indicate range of data and crosses indicate outliers. Raw data and *p*-values (ns; *p* > 0.05; *p* ≤ 0.05; ***p* ≤ 0.01; ** *p* ≤ 0.001; *** *p* ≤ 0.0001) are provided as a Source data file. Scale bars: 15 µm.
<figure><img src="image_12.png" /><figcaption>**Fig. 5** | Experimental validation of the guidance model. **a** P compartment cells overexpressing Ihog confronted to clones also overexpressing Ihog. White arrows indicate those cytonemes to be simulated. **b** Guidance field estimation (green small arrows) of the experimental conditions eliminating selected cytonemes. **c** In silico predicted trajectories for the simulated cytonemes. **a'** P compartment cells</figcaption></figure>

overexpressing Ihog confronted with clones overexpressing Dally. White arrows show the cytonemes to be simulated. **b'** Guidance field estimation (green small arrows) of the experimental conditions eliminating selected cytonemes. **c'** In silico predicted trajectories for the simulated cytonemes. Scale bars: 15µm.

proteins (see the table of parameters in section 4 of the Supplementary Note).

The potential $\phi(r, t)$, that contains the collective information for the orientation of cytonemes, is generated by the protein distribution $\rho(r, t)$ as follows:

$$ \phi(r, t) = (K^*\rho)(r, t) = \int K(|r - s|)\rho(t, s) ds. \quad (2) $$

Here, the potential kernel function $K$ defines the type of interactions between cytonemes according to the guiding factors and must be experimentally determined. In our case, since cytonemes only sense protein levels in the nearest cell membranes (Fig. 2b, c), the potential kernel function $K$ must be of short-range. This means that the guidance information is spatially limited to the neighboring cells (Fig. 2c). Even though this short-range phenomenon does not take place in the classical electric field, it is in agreement with most of the cases described in biological systems²⁵.

The choice of the sign in parameters $\beta_{\text{Ihog}} > 0$, $\beta_{\text{Dlp}} = \beta_{\text{Dally}} < 0$ in the charge distribution (1), and orientation field (3) has been carefully set in agreement with the above experimental results, that is, the observed effects of the Ihog-glypican interactions and their consequences on the stabilization and orientation of cytonemes. We refer to the discussion in section 1 of the Supplementary Note for further details.

Finally, from the potential $\phi(r,t)$ we can obtain the orientation vector field $O(r,t)$ that guides cytonemes:

$$ O(r, t) = -\nabla\phi(r, t) = -\int \nabla K(|r - s|)\rho(t, s) ds. \quad (3) $$

The dynamics of cytonemes is then deduced from a principle of least action. In other words, the equation of motion of each cytoneme

arises as a minimizer of the Lagrangian associated to its kinetic and the above (2) potential energy under the natural constraints of the biological system. This leads to a coupled system of non-linear and non-local hyperbolic equations for each cytoneme, which allows modeling the guidance process (see equation 12 in the Supplementary Note).

### In silico study of cytoneme orientation under overexpression conditions of Ihog and glypicans

In order to validate our model, we predicted in silico the trajectories of cytonemes protruding from experimental clones overexpressing Ihog (Fig. 5a) or glypicans (Fig. 5a') induced in the A compartment when confronted to P compartment cells expressing Ihog.

For simulations, we selected some of the experimental cytonemes (arrows in Fig. 5a-a'), removed them and predicted their orientation according to our model and within the simulated field (Fig. 5b-b'). Note that the orientation field corresponds to the frozen configuration after cytoneme removal and that the new simulated cytonemes grow and evolve in response to the global Ihog concentration.

The model simulations were able to predict the correct cytoneme orientation under conditions of overexpression of Ihog and glypicans (Fig. 5c-c'). It is important to note that the trajectories of simulated cytonemes are nearly identical to the experimental ones, even the curvature and orientation tendency were similar in both the in silico simulations (Fig. 5c-c') and in the real cases (Fig. 5a-a'). The temporal evolution of the cytonemes can be observed in more detail in Supplementary Movies 4 and 5.

### Model prediction of guidance for wild type cytonemes

Since our guidance model is able to predict cytoneme orientation in Ihog, Dally and Dlp overexpression conditions, we explored its prediction capacity for cytoneme orientation in wild type wing discs.
<figure><img src="image_13.png" /><figcaption>**Fig. 6 |** Wild type Ihog, Dally and Dlp spatial distributions at the basal plane of wing imaginal discs. **a** Wild type Ihog basal distribution in a Bac.IhogGFP wing disc. **b** Wild type Dlp basal distribution obtained by immunohistochemistry. **c** Wild type Dally basal distribution in a Dally.trap.YFP wing disc. **d** Characterization of Ihog (green), Dlp (blue) and Dally (red) profiles at the basal plane of imaginal wing discs. Dotted lines represent the average profile and shaded area and their experimental</figcaption></figure>

variability, quantified by the standard deviation over different samples ($n_{\text{Ihog}} = 17$, $n_{\text{Dlp}} = 17$, $n_{\text{Dally}} = 19$). **e** Surface representation of wild-type Ihog, Dlp and Dally basal patterns in the marked elliptical area of the wing pouch. These quantified 2D patterns (matrices) are the numerical inputs for the computational wild-type simulations. Scale bars: 30 µm.

To this end, the wild type spatial distribution of each protein was analyzed by fluorescence in Bac.IhogGFP wing discs (Fig. 6a), using the antibody for Dlp (Fig. 6b) and in Dally.trap.YFP (Fig. 6c) wing discs. Representative samples of their protein distributions have already been characterized in 1D profiles²⁴. Here we have made a statistical study of the profiles in several samples and examined the spatial correlation between Ihog, Dally and Dlp (Fig. 6d, 1D profiles). We observed that Ihog has a strong downregulation in the Hh reception region, Dlp increases its levels on the same region, and Dally is downregulated at both sides of the A/P compartment border.

A quantitative description in terms of 2D spatial distribution is needed for more accurate model simulations. To that end, we quantified the 2D spatial distribution of the basal protein pattern from confocal fluorescence images. In particular, we developed a script that allows visualization of the 2D protein patterns of the wing imaginal disc on 3D surfaces (Figs. 6e, 2D profiles) and generation of numerical data (matrices in Fig. 6e) for the levels of each protein, which have been used in (1) to model cytoneme orientation in wild type conditions. The 2D quantification showed that the wild-type modulation of the protein levels is not the same in all regions (Figs. 6e, 2D profiles), something that could not be deduced from 1D profiles.

We then tested if, using our guidance model, the wild-type spatial protein distributions of these proteins could be sufficient to predict cytoneme orientation within the Hh signaling region (Fig. 7). Our simulations show an increase in the parallel orientation of the field lines in the A compartment cells close to the A/P border (Fig. 7b, small red arrows) that change direction and become perpendicular to the A/P compartment border. Furthermore, as depicted in Fig. 7d, the model is able to predict the real observed cytoneme orientation in wild type cells (Fig. 7b).

To further study cytoneme behavior at the A/P compartment border region, we selected six cells in the Hh signaling area, three from the A and three from the P compartment. Cytonemes have a dynamic behavior of extension and retraction¹¹, although the signals triggering these dynamics have not yet been identified. For the simulations, we have considered two different plausible behaviors of retraction: cytonemes retracting after having contacted the tip (middle bottom cytonemes in Supplementary Movie 6) and cytonemes that retract after contacting other cytoneme while growing (middle top cytonemes in Supplementary Movie 6). The resulting predictions (Fig. 7d and Supplementary Movie 6) show that cytonemes protruding from the Hh-producing cells of the P compartment (green) are oriented towards the Hh receiving cells of the A compartment (red) to deliver the morphogen. At the same time, cytonemes protruding from receiving cells are oriented towards the producing cells to collect the morphogen.

## Discussion

Hh graded distribution across the receiving *Drosophila* epithelia is mediated by cytonemes¹². The transport of Hh is composed of two parallel and complementary processes. On the one hand, the transport along the cytonemes of vesicles containing Hh is associated with other proteins²⁶,²⁷; on the other, the elongation, orientation and contact of these cytonemes for reception¹¹,¹². In this work, we have proposed a model focused on the orientation of cytonemes.

It has been described that receptors, coreceptors and ligands are present in close proximity for the contact-dependent Hh reception involving both source and recipient cell cytonemes, in a way similar to a synaptic process¹¹,²⁷,²⁸. The purpose of contact would be the release of Hh at specific cytoneme sites. The molecular and cellular events that induce the interactions between cytonemes are critical determinants
<figure><img src="image_14.png" /><figcaption>**Fig. 7 |** *In silico simulation of cytoneme orientation in wild-type wing discs.*</figcaption></figure>

a Scheme depicting a wing imaginal disc. b Guidance field prediction (small arrows) in a wing disc using experimental wild-type distributions of Ihog, Dally and Dlp (top). Wild type orientation of P compartment cytonemes, marked with CD4-tomato (bottom). c The model predicts cytoneme orientation in wild-type conditions. d In silico predictions for cytoneme temporal growth on wild-type discs.

of signaling. Therefore, a cytoneme orientation and guidance mechanism must be needed for directed and efficient signaling. In this context, we have described that, in the Hh reception area, cytonemes are oriented perpendicular to the A/P compartment border. This preferential orientation, clearly shown by overexpressing Ihog, would give the shortest and most efficient contact-dependent reception. Using a multidisciplinary approach, we have developed a biomechanical guidance model to explain cytoneme orientation during Hh signaling.

The adhesion molecule Ihog plays a role in cytoneme dynamics since its overexpression stabilizes cytoneme temporal dynamics in different *Drosophila* tissues<sup>11,15,19</sup>. Here, we have analyzed if Ihog could also be involved in the control of cytoneme guidance. We first studied the effect of Ihog after confronting populations of cells expressing different levels of Ihog. The results have shown that cytoneme stabilization depends not only on the differential levels of Ihog in opposing cell populations but also on the distance between populations expressing different Ihog levels. We have observed that cytonemes are stabilized by Ihog overexpression even when Ihog is not present in the confronted cell population, indicating that Ihog-Ihog homophilic interaction in different cell populations is not sufficient for cytoneme orientation.

Since Dally and Dlp are also extracellular components of the Hh pathway, we investigated their possible involvement in the orientation of cytonemes. These proteins interact with Ihog for Hh signaling<sup>21,24,29</sup> cytoneme stabilization<sup>11,15</sup> and cell adhesion<sup>18,22</sup>. We studied cytoneme guidance using cytoneme stabilization as a tool and confronting clones with different levels of glypicans and Ihog. Cytonemes were stabilized and oriented whenever a cell population with high levels of Ihog was in the proximity of another cell population with high levels of either Dally or Dlp. However, cytonemes protruding from clones co-expressing Ihog and Dally or Ihog and Dlp were not stabilized. These results suggest a different cytoneme behavior depending on whether Ihog and glypicans interaction occurs in *trans* (in different cells) or in *cis* (in the same cells). Accordingly, glypican loss of function clones in the wing imaginal disc also indicates that *cis* interactions of Dally and Dlp with Ihog are not needed for Ihog-induced cytoneme stabilization, while *trans* interactions are necessary. Therefore, the specific protein levels of Ihog and glypicans in *cis* and *trans* can regulate both cytoneme stability and orientation. We then propose that the differences in levels of Ihog and glypicans could be responsible for cytoneme

orientation at the A/P compartment border and that these differences might be due to molecular competition of Ihog-Ihog homophilic and Ihog-glypicans heterophilic interactions. The “freezing” effect on cytonemes is probably due to Ihog sequestering glypicans in *trans*. In addition, the observed lack of cytoneme stabilization in confronted cell populations overexpressing Ihog might indicate that an excess of Ihog sequestering the glypicans available in *cis* could block the stabilization of the opposing cytonemes in *trans*.

It has recently been described that Ihog-Ihog interaction for cell-cell adhesion occurs via the fibronectin FNIII domains<sup>18</sup>, which are also responsible for the Ihog-glypican interaction: after expression of an Ihog mutant form lacking the Fn1 domain (IhogΔFn1), the ability of Ihog to form homodimers and to interact with glypicans disappears<sup>15,18</sup>. We have observed here that: 1) There is no cytoneme stabilization after expressing IhogΔFn1, in contrast to the cytoneme stabilization seen in cells overexpressing wild type Ihog. 2) Cytoneme stabilization, absent when Ihog-overexpressing cell populations are confronted, is recovered when cytonemes from cells overexpressing wild type Ihog are confronted to cells expressing IhogΔFn1. 3) Due to the lack of interaction, no strong stabilization is observed when cells expressing IhogΔFn1 are confronted to cells overexpressing glypicans.

The Fn1 domain has been described to be essential for both Ihog-Ihog and Ihog-Hh interactions<sup>18</sup>. Since the purpose of cytoneme contact recognition at reception is the release of Hh at specific cytoneme contact sites, the competition between Ihog-Ihog homophilic interaction and Ihog-Hh heterophilic interaction has been proposed to regulate the reception process at contact sites<sup>22</sup>. Interestingly, the expression of an Ihog form with point mutations in the Fn1 domain that abrogate the interaction with Hh but not with glypicans stabilizes and orients cytonemes as the overexpression of wild-type Ihog does<sup>15,18</sup>. Therefore, although both Ihog-Hh and Ihog-Ptc interactions are essential for Hh transfer during reception, Hh by itself does not seem to be necessary for the generation or the stabilization of cytonemes<sup>11</sup>. However, the expression patterns of Ihog, Dlp and Dally are regulated in Hh signaling regions, probably due to a modulation by the Hh pathway itself<sup>24,30</sup>. Consequently, a regulatory loop might be established in which the Hh responses could indirectly influence the orientation of cytonemes. This autoregulatory loop is consistent with the one described for the orientation of cytonemes in the FGF signaling pathway in *Drosophila*<sup>31–33</sup>.

Our experimental results permit the generation of a framework for cytoneme guidance dependent on different levels of Ihog, Dally and Dlp. To this end, using a multidisciplinary approach, we designed a mathematical model based on the potential theory in physics, which resembles the observed favored/disfavored regions for cytoneme stabilization. We validated our model comparing the in silico predictions with the experimental results obtained from clones expressing different levels of Ihog and glypicans. The model simulations were able to reproduce with high accuracy the experimental cytoneme trajectories.

To close the circle, after validating our model in gain of function conditions, we studied the orientation of cytonemes observed in the wild-type wing imaginal disc. Our model predictions showed that the wild type protein distributions of Ihog, Dally and Dlp can generate a field of directions that can be enough to orient both producing and receiving cytonemes. In other words, for proper cytoneme guidance, the interactions of glypicans and Ihog might allow the exchange of information about cellular relative locations. Therefore, we propose a new role of Ihog and glypicans in cell-cell signaling, influencing the orientation of cytonemes for an efficient contact during the reception process. Inefficient cytomeme contacts could explain the lack of Hh responses in $ttv^{-/−}$ *btv^{-/−}* double mutant clones located in the reception region<sup>16,17</sup>.

Our model is general enough to be applied to predict and study cytoneme guidance in other signaling processes. For instance, we have observed that cytonemes are also oriented perpendicular to the D/V
compartment border, where Wg and Notch, but not Hh, are expressed³⁴. This could suggest that Wg and Notch signaling can be dependent on cytonemes in the wing disc, as it has been shown in other *Drosophila* tissues³⁵⁻³⁷ and in vertebrates³⁸⁻⁴⁰. Since glypicans have a specific expression pattern at the D/V border⁴¹, they might help the driving of cytonemes perpendicular to this border.

Likewise, our model can also be applicable to vertebrate systems. Other receptor proteins with structural domains similar to Ihog, such as DCC and Robo, have been described to be important in axon guidance⁴², what emphasizes the similarities between cytoneme-mediated cell signaling and neuronal communication¹¹,²⁷,²⁸. Moreover, the role in neuronal guiding of some extracellular matrix proteins, including glypicans, has been described⁴³,⁴⁴. Similar orientation of cytonemes after ectopic expression of CDO and BOC (vertebrate homologs of Ihog and Boi) has been reported for Shh signaling in the chicken limb bud⁴⁵, suggesting that the control mechanisms to orient cytonemes might be evolutionarily conserved.

## Methods

### Experimental material and methods

**Drosophila stock maintenance.** *Drosophila melanogaster* stocks were maintained according to protocols described in Ashburner manual⁴⁶. In particular, flies were reared under standard lab conditions at 25 °C in vials and fed with yeast for their life circle. For transient expression of transgenic constructs, the fly crosses were maintained at 18 °C and then incubated at restrictive temperature (29 °C) to inactivate the Gal80ts repressor 24 hours before dissection. The description of mutations, insertions and transgenes is available at Fly Base (http://flybase.org).

**Binary expression systems.** The transient expression (overexpression or downregulation) of proteins was performed under the control of the binary systems: Gal4/UAS⁴⁷ LexA/LexAop⁴⁸. Although each binary system works in a similar manner, each system only recognizes its specific target sequences. This specificity allows the combination of both systems in the same sample.

For those cases in which there are not specific drivers to express the transgenes in a specific tissue region, or in cases in which the mutant phenotype is very invasive, it is possible to study the mutant condition limited to certain cell populations (clones). This is achievable through the use of the FLP/FRT system⁴⁹, which uses the enzyme flippase (FLP) to induce recombination between FRTs (Flippase Recombination Target) sequences. The mitotic recombination was induced through a heat shock (HS) at 37°C at a specific developmental stage and the *Drosophila* tissue cells respond by generating at random mosaic recombinant clones with a genotype different than the rest of the tissue. In combination with the previous binary system, we generated random clones in which the transgenes are activated. In this work, we have used the CoinFLP system⁵⁰, which is a useful Flip-out recombination system that introduces another FRT (FRT3) to allow the induction of two types of clones, each one with a different driver. In our case, we have used CoinFLIP-LexA::GAD.Gal4 to induce clones expressing either the Gal4 or the LexA drivers.

The following drivers were used to induce ectopic expression using the Gal4/UAS and LexA/LexO systems: tubGal80ts, hs-Flp122 (Bloomington Drosophila Stock Center, BDSC), hs-FLP⁴⁹, hh.Gal4⁵¹, F.Act > CD2 > Gal4⁵², Ptc Gal4⁵³; ap.Gal4⁵⁴, hh.LexA¹¹, generated by Ernesto Sánchez-Herrero, CBMSO), CoinFLIP-LexA::GAD.Gal4 (BDSC 59270) and dpp.LexA⁵⁵.

### Overexpression stocks

**The UAS-transgene strains.** UAS.ihog-YFP²¹, UAS.ihog-CFP²⁹, UAS.CD4-tdTomato⁵⁶, UAS.ihog-RNAi (VDRC 102602), UAS.boi-RNAi (VDRC 108265), UAS.LifeActGFP (BDSC 35544), UAS.LifeActRFP (BDSC 58362), UAS.dlp-GFP⁵⁷, UAS.dally-GFP⁵⁸, UAS.GMA-GFP (Actin-binding

domain of moesin tagged with GFP)⁵⁹ and UAS.CD63-Cherry (Provided by Clive Wilson, Oxford University).

**The LexAop-transgene strains.** LexAop.ihog-RFP¹¹. The LexAop.ihogΔFn1-RFP construct used IhogΔFn1-RFP from the pTWR vector¹⁵, which was introduced into the pLOTattB plasmid⁴⁸ carrying the lexA operator (LexAop). Transgenic strains were recovered using standard protocols.

**Other stocks.** Dally.trap.YFP (DGRC 115511), EnhancerPtcRed (Kyoto stock center, DGRC 109138), Bac.IhogGFP: Rp49 promoter³⁰, ttv ⁵²⁴ ¹⁷, and botv⁵¹⁰ ¹⁷. Since the Ihog homolog, Brother of Ihog (Boi) can substitute Ihog in some functions⁶⁰, Ihog downregulation experiments were also done by expressing together Ihog (ref: 102602) and Boi (ref: 108265) RNAis (VDRC; http://stockcenter.vdrc.at).

**Transient expression of transgenes.** The binary system can be used under a temporal control of the GAL80ts that acts as a Gal4 repressor blocking the interaction between Gal4 and the UAS sequence. The GAL80ts have also been described to repress the LexA/LexAop system when it contains the Gal4 sequences repressible by GAL80. For transient expression of transgenic constructs, the crosses were maintained at 18 °C and then incubated at the restrictive temperature (29 °C) to inactivate the Gal80ts repressor 24–48 h before dissection. After dissection, third instar imaginal discs were fixed in 4% paraformaldehyde (PF) in PBS for 20 min at room temperature (RT) and mounted in Vectashield medium (Vector labs).

Mutant clones were generated, as described, by a flippase induced by heat shock (HS) incubating larvae at 37 °C. The HS was given at 48–72 h after egg laying (AEL) for 7–10 min in the case of the CoinFLIP-LexA::GAD.Gal4 system and for 10–15 min in the case of the F.Act > CD2 > Gal4 system. *ttv*⁻/⁻ *botv*⁻/⁻ double mutant clones (*ttv* ⁵²⁴ *botv*⁵¹⁰) were generated using FRT42D cassette under the control of hsFlp¹²² flippase and individuals were incubated at 37 °C for 45 min.

**Immunohistochemistry of wing imaginal discs.** Immunostaining was performed according to standard protocols⁶¹. Imaginal discs from third instar larvae were fixed in 4% paraformaldehyde (PFA) in PBS for 20 min at room temperature (RT) and permeabilized with 0,1% Triton in PBS (PBT) before incubation with 1% BSA (bovine serum albumin) in PBT for blocking the tissue to decrease the unspecific background (1 h at RT) and follow by primary antibody incubations (overnight at 4 °C). Incubation with fluorescent secondary antibodies was performed for 1 h at RT and then washed and mounted in mounting media (Vectashield).

The mouse monoclonal α-Dlp primary antibody (from¹³) and the rat monoclonal α-Ci primary antibody (from⁶²) were used at 1:30 and 1:20 dilutions, respectively, for regular immunostaining and/or extracellular immunostaining. The protocol for the extracellular labelling was done following⁶³; imaginal discs from third instar larvae were dissected on ice, transferred immediately to ice-cold M3 medium containing α-GFP (rabbit anti GFP polyclonal antibody from Chromotek, ref: PABG1, at 1:1000 dilution) and α-Dlp (mouse monoclonal antibody from¹³ at 1:30 dilution) antibodies and incubated at 4 °C for 1 hour. The incubation with the primary antibody was done under these ex vivo conditions, without detergents prior fixation, preventing antibody penetration inside cells. Imaginal discs were then washed in ice-cold PBS, fixed in PBS 4% PF at 4 °C, washed in PBT and incubated with the secondary fluorescent antibody as above. The secondary antibodies used in this work were anti-mouse Pacific Blue from ThermoFischer (ref: P-31582) at 1:400 dilution and anti-rabbit Alexa 647 from ThermoFischer (ref: A-31573) at 1:400 dilution.

**Data acquisition for fixed, in vivo and ex vivo experiments.** Laser scanning confocal microscopes (LSM700 and LSM800 Zeiss) were
used for confocal fluorescence imaging of the fixed imaginal discs. The in vivo abdominal histoblasts imaging was performed in a chamber to seat and orient the pupae to look under the microscopy as described in⁶⁴. The dorsal abdominal segment A2 was filmed using 40x magnification; Z-stacks of around 10 µm of thickness with a step size of 0.8 µm were taken, every 2 min, using a LSM800 confocal microscope. For the ex vivo experiments we induced 40 h of Ihog overexpression in P posterior cells and we dissected third instar larvae as described in⁶⁴. The ex vivo imaging was performed in standard microscope slides, where the discs were mounted in M3 insect medium and covered with a 15 × 15 mm coverslip. The samples were filmed using 40x magnification; Z-stacks of around 8 µm of thickness with a step size of 0.62 µm were taken every minute using a LSM800 confocal microscope.

**Data analysis and statistical study.** The data of abdominal histoblast nests are published in¹¹. The quantification of the number of cytosemes and their orientation angles have been obtained manually using Fiji. The data representation, in boxplots and rose polar diagrams, was performed using Matlab2015a functions. For the statistical study of the data in the figures, we first performed a Shapiro-Wilk test of normality, since the S-W test showed that our data distributions were not parametric, to study the statistical significance between datasets, we used a Wilcoxon rank sum test per pairs in Matlab2015a (See Supplementary excel file for the quantified data and the calculated p-values). Finally, some cytosemes were marked in ex vivo and in vivo experiments to better show the dynamics and the stabilization of cytosemes, those tracks were performed using the Fiji plugging: Manual Tracking.

**Quantification of cytoneme orientation.** The experimental data of cytoneme length and orientation due to Ihog overexpression in LexA and Gal4 clones were quantified manually using FJJI along 235 cytosemes from different clones (n = 21). In each measurement, we quantified the length (λ) and the angle (β) of a cytoneme and the distance (d) and angle (γ) of the tip of that cytoneme with the closest clone. The study of the effect in the length was performed analyzing λ vs d, and the study of the orientation was calculated using α = γ - β. Please notice that the increase in α means that cytosemes are avoiding the closest confronted clone. Limit cases, α = 0 and α = 90, can be respectively interpreted as no change in the direction and opposite direction.

**Quantification of Ihog, Dally and Dlp protein distribution profiles.** The profiles of each protein were obtained measuring, with Fiji profile tool, the fluorescence intensity over a ROI of 15–80 µm² in n = 17 Bac.IhogGFP wing discs, where the Dlp profile was also measured in the same ROI using α-Dlp primary antibody¹³ and Dally profiles were obtained from the same ROI size in n = 19 Dally.trap.YFP wing discs. (See supplementary excel file for the quantified raw data). The 1D average profiles, together with the experimental variability (measured with the standard deviation of the quantified profiles), were computed following the mathematical protocol described in²⁰.

## Theoretical material and methods

Since cytosemes mainly move over the basal plane of the tissue we can simplify our model to two spatial dimensions as considered above. In addition, since Ihog and glypican proteins are attached to cell membranes, we assumed that the guiding factors are spatially constrained to cytosemes and cells. In order to deduce the cytoneme kinetics (behavior), we have chosen the variational principle approach, which facilitates a robust mathematical formulation of the equations. Note that cytosemes are dynamic structures and the amount of proteins located on those cytoneme membranes changes spatially over time. This implies that the orientation field is dynamic and not only the protein concentration at cell membranes is responsible for the orientation of cytosemes, but also the motion of cytosemes themselves can

modify the generated field. This feedback loop, coupled with the biological complexity of the process, introduces a large non-linearity in the model, which forces us to solve the equations via numerical methods. Here, we briefly describe the theoretical and numerical protocol used for in silico simulations, a detailed description can be found in Supplementary Note.

Note that the dynamics of cytosemes cannot be reduced simply to each single protein moving freely under the force field $O(r, t) = -\nabla\phi(r, t)$. In fact, we must take into account the geometric information that filopodia structure imposes spatial constraints to those proteins that are attached to their membranes. We return here to the aforementioned idea suggesting that energy must be optimized during the movement of cytosemes according to the principle of least action. Let $L(t)$ be the variable length of the curve $\gamma(\xi, t)$ that represents the cytoneme. Then, the minimization problem

$$ \min_{\gamma \in \mathcal{G}} \int_0^T \int_0^L \left( \underbrace{\frac{1}{2} |\dot{\gamma}(\xi, t)|^2}_{\text{kinetic energy}} - \underbrace{\bar{\rho}\phi(\gamma(\xi, t), t)}_{\text{potential energy}} \right) d\xi dt, $$

$$ \mathcal{G} = \{\gamma \in \mathcal{C}^1 : |\dot{\gamma}(\xi, t)| = 1, + \text{boundary conditions}\}, $$

leads to the differential equation that governs the dynamics of the cytoneme, (see supplementary note for a more exhaustive description):

$$ \frac{\partial}{\partial t}(\dot{\gamma}) = \frac{\partial}{\partial \xi}(\lambda\dot{\gamma}) - \frac{1}{t}\dot{\gamma} - \bar{\rho}\nabla\phi(\gamma, t), \quad (4) $$

being τ the relaxation time due to linear friction of the cytosemes with the extracellular matrix. Here, ρ is the distribution of proteins over the specific cytoneme γ, whilst the orientation potential φ is collectively generated by all protein distribution in the system according to (1).

In the sequel, we present a discretized version of the dynamics of cytosemes that allowed us to computationally solve the previous continuous model (4). In order to do this, we defined the discretized cytoneme as an inextensible chain of N(t) connected nodes carrying a specific protein concentration. Then, we used generalized coordinates in order to eliminate inextensibility constraints. Hence, the position $\mathbf{r}_i(t)$ of the i-th node of a cytoneme $\gamma(\xi, t)$ at time t is determined from the position $\mathbf{r}_{i-1}(t)$ of the (i - 1)-th node by $\mathbf{r}_i(t) = \mathbf{r}_{i-1}(t) + l_i(\cos \theta_i(t), \sin \theta_i(t))$, for some angle $\theta_i(t)$ and some length $l_i$ between nodes i and i - 1. Then, we obtain:

$$ \mathbf{r}_i(t) = \mathbf{r}_0 + \sum_{j=1}^{j=i} l_j(t) (\cos \theta_j(t), \sin \theta_j(t)), \quad i=1, \dots, N(t), \quad (5) $$

where $\mathbf{r}_0$ denotes the base point where the cytoneme is attached to the cell membrane. Note that since cytosemes can elongate and retract, the length $l_{N(t)}$ of the last piece of the chain can grow (up to a maximal length) or decrease, and the number N(t) of nodes will change with time. Then, we can discretize the continuous differential Eq. (4) at each discrete parameter value $\xi_i = \sum_{j \le i} l_j$ representing the i-th bond $\mathbf{r}_i(t) \approx \mathbf{y}(t, \xi_i(t))$. Then, using the change of variables in (5), and the generalized coordinates $\mathbf{\Theta}(t) = (\theta_1(t), \dots, \theta_{N(t)}(t))$ and $\dot{\mathbf{\Theta}}(t) = (\dot{\theta}_1(t), \dots, \dot{\theta}_{N(t)}(t))$ the constrained dynamics of the whole discrete cytoneme in generalized coordinates are:

$$ \mathbf{M}(\mathbf{\Theta}(t), t)\ddot{\mathbf{\Theta}}(t)^T = -\frac{1}{\tau}\mathbf{M}(\mathbf{\Theta}(t), t)\dot{\mathbf{\Theta}}(t)^T + \mathbf{G}(\mathbf{\Theta}(t), t)\dot{\mathbf{\Theta}}(t), \quad (6) $$

Here M and G are a matrix and a vector, respectively, whose explicit expression is derived in the Supplementary Note. In particular,
**G** contains all the collective information of the orientation potential generated by the full network of cytonemes. An alternative but equivalent derivation of (6) could have been to implement a discrete version of the principle of least action. Specifically, for the discrete kinetic energy $E_{Kin}(\mathbf{\Theta}(t), \mathbf{\Theta}(t))$ and potential energy $E_{Pot}(\mathbf{\Theta}(t),t)$ in generalized coordinates, we can formulate the discrete minimization problem of the corresponding discrete action functional $\mathcal{L} = E_{Kin} - E_{Pot}$ and (6) arises as its corresponding Euler-Lagrange equation. Although in this paper, for simplicity, we have set the relaxation time $\tau$ as a constant, we could also introduce modulation factors that modify $\tau$. For example, it might include the required time to stabilize the cytonemes in the membrane transmitting the information to other cells. In that case, the appropriate modulation should depend, among other factors, on the concentration of Hedgehog that the cytoneme carries, the Ptc concentration of the cell membrane, the co-receptor Ihog or other proteins with adhesion properties.

The matrix Eq. (6) contains the dynamics of each cytoneme. In our numerical simulations, we solved it numerically by first inverting *M* with LU-decomposition and then computing the time evolution of the guidance with a Runge-Kutta 4th order method for all the cytoneme nodes at the same time (see Supplementary Movie 7 for an illustrative example of the numerical resolution).

**Experimental data and model inputs for computational simulations - Model simulations for clone experiments.** To study and validate our model prediction in mutant conditions, we first selected an experimental clone, obtained previously by confocal imaging. We then selected an area of interest in the clone, that contained the cytonemes that we wanted to simulate. In order to maintain the closest conditions to the experimental case, the selected cytonemes were removed from the image with FIJI tools with the following protocol:

* Select the cytoneme area manually with the Freehand ROI tool.

* Save the ROI.

* Save the cytoneme origin position.

* Measure the average and standard deviation of the signal in the surrounding background area around previous ROI.

* Multiply the cytoneme ROI by zero to remove the cytoneme concentration.

* Add the average background signal to the ROI.

* Introduce a Gaussian noise with the standard deviation background value.

The resulting image is the one obtained experimentally but devoid of cytonemes. Next, we introduced the protein concentration of the experimental clones in (1) using the pixel distribution, and we simulated the path that each removed cytoneme would follow according to our model, departing from positions of origin previously saved.

**- Model simulation for wild-type wing imaginal disc.** The experimental protein concentrations $p$ (1) were obtained from the pixel intensities in each channel in gray-scale units of the experiments described in the Results section. For wild-type simulations, the protein spatial distribution was measured in the elliptical area at the basal plane of the wing pouch; the elliptical shapes were scaled to the same lattice to be uploaded in the in silico wing disc. Dlp pattern was quantified using Dlp monoclonal antibody in the wing imaginal disc expressing either Ihog tagged with GFP (Bac.IhogGFP) or Dally tagged with YFP (Dally.trap.YFP). Extracellular staining was performed for Dally tag and Dlp protein to confirm that the same pattern is present in the extracellular matrix where the cytonemes sense the levels of glypicans and Ihog.

Since our tissue model works in 2D and Ihog and glypicans proteins are present along all tissue membranes, we obtained the quantitative description of 2D spatial distribution at the basal plane of the imaginal wing disc, plane where cytonemes are present. To that end,

the data were treated using a computational script that allows extracting the experimental data from confocal images and visualizing it in a 3D surface the 2D spatial distribution of the particular fluorescence protein pattern along the disc tissue. The script also generates numerical matrices that contain the quantifications of the levels of a specific protein in different regions. Those matrices are the experimental inputs that we used in our model to simulate in silico cytoneme orientation.

## Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

## Data availability

Source data generated in this study is provided with this paper in the Source Data file, and more supplementary information from the corresponding authors upon request. Source data are provided with this paper.

## Code availability

The cytoneme simulation code was performed using a custom-made code on C# language using Unity Engine. The source code is available at https://github.com/AdrianA-T/Guidance-Model. More details of this computational code can be obtained from the corresponding authors on request.

## References

1. Lawrence PA. Morphogens: How big is the big picture? *Nat Cell Biol.* **3**, https://doi.org/10.1038/35083096 (2001).

2. Kondo, S. & Miura, T. Reaction-diffusion model as a framework for understanding biological pattern formation. *Sci. (80-).* **329**, 1616-1620 (2010).

3. Akiyama, T. & Gibson, M. C. Morphogen transport: theoretical and experimental controversies. *Wiley Interdiscip. Rev. Dev. Biol.* **4**, 99-112 (2015).

4. Ramírez-Weber, F. A. & Kornberg, T. B. Cytonemes: cellular processes that project to the principal signaling center in Drosophila imaginal discs. *Cell* **97**, 599-607 (1999).

5. Gradilla, A. C. & Guerrero, I. Cytoneme-mediated cell-to-cell signaling during development. *Cell Tissue Res.* **352**, 59-66 (2013).

6. Kornberg, T. B. Cytonemes and the dispersion of morphogens. *Wiley Interdiscip. Rev. Dev. Biol.* **3**, 445-463 (2014).

7. González-Méndez L, Gradilla AC, Guerrero I. The cytoneme connection: Direct long-distance signal transfer during development. *Dev.* **146**, https://doi.org/10.1242/dev.174607 (2019).

8. Sagar, PrölsF., Wiegreffe, C. & Scaal, M. Communication between distant epithelial cells by filopodia-like protrusions during embryonic development. *Dev* **142**, 665-671 (2015).

9. Pröls, F., Sagar & Scaal, M. Signaling filopodia in vertebrate embryonic development. *Cell Mol. Life Sci.* **73**, 961-974 (2016).

10. Garcia-Bellido, A., Ripoll, P. & Morata, G. Developmental compartmentalisation of the wing disk of drosophila. *Nat. N. Biol.* **245**, 251-253 (1973).

11. González-Méndez L, Seijo-Barandiarán I, Guerrero I. Cytoneme-mediated cell-cell contacts for hedgehog reception. *Elife.* **6**, https://doi.org/10.7554/eLife.24045 (2017).

12. Bischoff, M. et al. Cytonemes are required for the establishment of a normal Hedgehog morphogen gradient in Drosophila epithelia. *Nat. Cell Biol.* **15**, 1269-1281 (2013).

13. Lum, L. et al. Identification of Hedgehog pathway components by RNAi in Drosophila cultured cells. *Sci.* (80-). **299**, 2039-2045 (2003).

14. Beachy, P. A., Hymowitz, S. G., Lazarus, R. A., Leahy, D. J. & Siebold, C. Interactions between Hedgehog proteins and their binding partners come into view. *Genes Dev.* **24**, 2001-2012 (2010).
15. Simon E, et al. Glypicans define unique roles for the Hedgehog co-receptors boi and ihog in cytoneme-mediated gradient formation. *Elife.* **10**, https://doi.org/10.7554/ELIFE.64581 (2021).

16. Lin, X. Functions of heparan sulfate proteoglycans in cell signaling during development. *Development* **131**, 6009–6021 (2004).

17. Takei, Y. et al. genes shape morphogen gradients through synthesis of heparan sulfate proteoglycans. *Development* **131**, 73–82 (2004).

18. Yang S, et al. Competitive coordination of the dual roles of the hedgehog co-receptor in homophilic adhesion and signal reception. *Elife.* **10**, https://doi.org/10.7554/ELIFE.65770 (2021).

19. Rojas-Ríos, P., Guerrero, I., González-Reyes, A. Cytoneme-mediated delivery of Hedgehog regulates the expression of bone morphogenetic proteins to maintain germline stem cells in Drosophila. *PLoS Biol.* **10**, https://doi.org/10.1371/journal.pbio.1001298 (2012).

20. Aguirre-Tamaral, A. & Guerrero, I. Improving the understanding of cytoneme-mediated morphogen gradients by in silico modeling. *Haugh JM, ed.* PLOS Comput Biol. **17**, e1009245 (2021).

21. Callejo, A. et al. Dispatched mediates Hedgehog basolateral release to form the long-range morphogenetic gradient in the Drosophila wing disk epithelium. *Proc. Natl Acad. Sci. USA* **108**, 12591–12598 (2011).

22. Wu, X. et al. The Drosophila Hedgehog receptor component Interference hedgehog (Ihog) mediates cell–Cell interactions through trans-homophilic binding. *J. Biol. Chem.* **294**, 12339–12348 (2019).

23. Riedl, J. et al. Lifeact mice for studying F-actin dynamics. *Nat. Methods* **7**, 168–169 (2010).

24. Yan, D. et al. The cell-surface proteins Dally-like and Ihog differentially regulate Hedgehog signaling strength and range during development. *Development* **137**, 2033–2044 (2010).

25. Co, A. D., van Vliet, S., Kiviet, D. J., Schlegel, S. & Ackermann, M. Short-range interactions govern the dynamics and functions of microbial communities. *Nat. Ecol. Evol.* **4**, 366–375 (2020).

26. Gradilla, A.-C. et al. Exosomes as Hedgehog carriers in cytoneme-mediated transport and secretion. *Nat. Commun.* **5**, 5649 (2014).

27. González-Méndez L, et al. Polarized sorting of Patched enables cytoneme-mediated Hedgehog reception in the Drosophila wing disc. *EMBO J.* **39**, https://doi.org/10.15252/embj.2019103629 (2020).

28. Huang, H., Liu, S. & Kornberg, T. B. Glutamate signaling at cytoneme synapses. *Sci.* (80-). **363**, 948–955 (2019).

29. Bilioni, A. et al. Balancing Hedgehog, a retention and release equilibrium given by Dally, Ihog, Boi and shifted/DmWif. *Dev. Biol.* **376**, 198–212 (2013).

30. Hsia, E. Y. C. et al. Hedgehog mediated degradation of Ihog adhesion proteins modulates cell segregation in Drosophila wing imaginal discs. *Nat. Commun.* **8**, 1–15 (2017).

31. Du, L., Sohr, A., Yan, G., & Roy, S. Feedback regulation of cytoneme-mediated transport shapes a tissue-specific FGF morphogen gradient. *Elife.* **7**, https://doi.org/10.7554/ELIFE.38137 (2018).

32. Patel, A. et al. Cytonemes coordinate asymmetric signaling and organization in the Drosophila muscle progenitor niche. *Nat. Commun.* **2022** **131**, 1–17 (2022).

33. Du L, Sohr A, Li Y, Roy S. GPI-anchored FGF directs cytoneme-mediated bidirectional contacts to regulate its tissue-specific dispersion. *Nat Commun.* **13**, https://doi.org/10.1038/S41467-022-30417-1 (2022).

34. Curtiss, J., Halder, G. & Mlodzik, M. Selector and signalling molecules cooperate in organ patterning. *Nat. Cell Biol.* **2002** **43**, 4, E48–E51 (2002).

35. De Joussineau, C. et al. Delta-promoted filopodia mediate long-range lateral inhibition in Drosophila. *Nature* **426**, 555–559 (2003).

36. Cohen, M., Georgiou, M., Stevenson, N. L., Miodownik, M. & Baum, B. Dynamic filopodia transmit intermittent delta-notch signaling to

drive pattern refinement during lateral inhibition. *Dev. Cell.* **19**, 78–89 (2010).

37. Huang H, Kornberg TB. Myoblast cytonemes mediate Wg signaling from the wing imaginal disc and Delta-Notch signaling to the air sac primordium. *Elife.* **4**, (MAY). https://doi.org/10.7554/eLife.06114 (2015).

38. Hamada, H. et al. Involvement of Delta/Notch signaling in zebrafish adult pigment stripe patterning. *Development* **141**, 1418–1418 (2014).

39. Eom DS, Bain EJ, Patterson LB, Grout ME, Parichy DM. Long-distance communication by specialized cellular projections during pigment pattern development and evolution. *Elife.* **4**, (DECEMBER 2015) https://doi.org/10.7554/eLife.12401 (2014).

40. Stanganello, E. et al. Filopodia-based Wnt transport during vertebrate tissue patterning. *Nat. Commun.* **6**, 5846 (2015).

41. Han, C., Yan, D., Belenkaya, T. Y. & Lin, X. Drosophila glypicans Dally and Dally-like shape the extracellular Wingless morphogen gradient in the wing disc. *Development* **132**, 667–679 (2005).

42. Dickson, B. J. Molecular mechanisms of axon guidance. *Sci.* (80-). **298**, 1959–1964 (2002).

43. Rawson, J. M. et al. The heparan sulfate proteoglycans Dally-like and Syndecan have distinct functions in axon guidance and visual-system assembly in Drosophila. *Curr. Biol.* **15**, 833–838 (2005).

44. Bülow, H. E. et al. Extracellular sugar modifications provide instructive and cell-specific information for axon-guidance choices. *Curr. Biol.* **18**, 1978–1985 (2008).

45. Sanders, T. A., Llagostera, E. & Barna, M. Specialized filopodia direct long-range transport of SHH during vertebrate tissue patterning. *Nature* **497**, 628–632 (2013).

46. Ashburner M. Drosophila: A laboratory Handbook. Published online 1989.

47. Brand, A. H. & Perrimon, N. Targeted gene expression as a means of altering cell fates and generating dominant phenotypes. *Development* **118**, 401–415, https://pubmed.ncbi.nlm.nih.gov/8223268/ (1993).

48. Lai, S. L. & Lee, T. Genetic mosaic with dual binary transcriptional systems in Drosophila. *Nat. Neurosci.* **9**, 703–709 (2006).

49. Golic, K. G. & Lindquist, S. The FLP recombinase of yeast catalyzes site-specific recombination in the drosophila genome. *Cell* **59**, 499–509 (1989).

50. Bosch, J. A., Tran, N. H. & Hariharan, I. K. Coinflp: a system for efficient mosaic screening and for visualizing clonal boundaries in drosophila. *Dev.* **142**, 597–606 (2015).

51. Tanimoto, H., Itoh, S., Ten Dijke, P. & Tabata, T. Hedgehog creates a gradient of DPP activity in Drosophila wing imaginal discs. *Mol. Cell.* **5**, 59–71 (2000).

52. Pignoni, F. & Zipursky, S. L. Induction of Drosophila eye development by decapentaplegic. *Development* **124**, 271–278 (1997).

53. Hinz, U., Giebel, B. & Campos-Ortega, J. A. The basic-helix-loop-helix domain of Drosophila lethal of scute protein is sufficient for proneural function and activates neurogenic genes. *Cell* **76**, 77–87 (1994).

54. Calleja, M., Moreno, E., Pelaz, S. & Morata, G. Visualization of gene expression in living adult drosophila. *Sci.* (80-). **274**, 252–255 (1996).

55. Yagi, R., Mayer, F. & Basler, K. Refined LexA transactivators and their use in combination with the Drosophila Gal4 system. *Proc. Natl Acad. Sci. USA* **107**, 16166–16171 (2010).

56. Han, C., Jan, L. Y. & Jan, Y. N. Enhancer-driven membrane markers for analysis of nonautonomous mechanisms reveal neuron-glia interactions in Drosophila. *Proc. Natl Acad. Sci. USA* **108**, 9673–9678 (2011).

57. Han, C. H., Belenkaya, T. Y., Wang, B. & Lin, X. Drosophila glypicans control the cell-to-cell movement of Hedgehog by a dynamin-independent process. *Development* **131**, 601–611 (2004).
58. Eugster, C., Panáková, D., Mahmoud, A. & Eaton, S. Lipoprotein-Heparan Sulfate Interactions in the Hh Pathway. *Dev. Cell.* **13**, 57-71 (2007).

59. Bloor, J. W. & Kiehart, D. P. zipper Nonmuscle Myosin-II Functions Downstream of PS2 Integrin in Drosophila Myogenesis and Is Necessary for Myofibril Formation. *Dev. Biol.* **239**, 215–228 (2001).

60. Camp, D., Currie, K., Labbé, A., Meyel, D. J. & van, Charron, F. Ihog and Boi are essential for Hedgehog signaling in Drosophila. *Neural Dev.* **5**, 28 (2010).

61. Capdevila J, Guerrero I. Targeted expression of the signaling molecule decapentaplegic induces pattern duplications and growth alterations in Drosophila wings. *EMBO J.* **13**, 4459, (accessed 2 August 2021); /pmc/articles/PMC395378/?report= abstract (1994).

62. Motzny, C. K. & Holmgren, R. The Drosophila cubitus interruptus protein and its role in the wingless and hedgehog signal transduction pathways. *Mech. Dev.* **52**, 137–150 (1995).

63. Torroja, C., Gorfinkiel, N. & Guerrero, I. Patched controls the Hedgehog gradient by endocytosis in a dynamin-dependent manner, but this internalization does not play a major role in signal transduction. *Development* **131**, 2395–2408 (2004).

64. Seijo-Barandiarán, I., Guerrero, I., & Bischoff, M. In vivo imaging of hedgehog transport in drosophila epithelia. In: *Methods in Molecular Biology*. Vol. 1322. Humana Press Inc. 9–18 https://doi.org/10.1007/978-1-4939-2772-2_2 (2015).

## Acknowledgements

We are grateful to and Pedro Ripoll and Ana-Citlali Gradilla for comments on the manuscript. We also thank to Laura González-Méndez, Carlos Jiménez-Jiménez, David Sánchez-Hernández and Eléanor Simon in I.G laboratory for their help with some experiments, the Confocal Facilities of the CBMSO and Bloomington and Vienna stock centers for fly stocks. This work was supported by grants BFU2017-83789-P, PID2020-114533GB-C21 and TENTACLES consortium RED2018-102411-T to I.G from the Spanish Ministry of Science, Innovation and Universities and by institutional grants from the Fundación Areces and Banco de Santander to the CBMSO. FPI fellowship from the Spanish Ministry of Science, Innovation and Universities supported A.A-T (BFU2017-83789-P). This work was also supported by grants RTI2018-098850-B-I00 to J.S from the MINECO-Feder (Spain), PY18-RT-2422 & A-FQM-311-UGR18 to M.C, D.P and J.S from the Junta de Andalucía (Spain), MECD (Spain) research grant FPI2015/074837 to M.C, and partially supported by the MECD (Spain) research grant FPU14/06304 and the European Research Council (Europe) Project ERC-COG-2019 WACONDY (grant agreement No 865711) to D.P.

## Author contributions

A.A.-T. performed *Drosophila* experiments, immunostaining and imaging of confocal images of Drosophila imaginal disc and in vivo abdominal histoblast nets, as well as data quantification and statistical analysis, and wrote the manuscript. M.C. and D.P. developed the mathematical model and wrote the supplementary note. M.C. developed the computational code and simulations. J.S. designed and supervised the mathematical framework and wrote the manuscript. I.G. designed and supervised the *Drosophila* experiments and wrote the manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-022-33262-4.

**Correspondence** and requests for materials should be addressed to Juan Soler or Isabel Guerrero.

**Peer review information** *Nature Communications* thanks the anonymous reviewers for their contribution to the peer review of this work. Peer reviewer reports are available.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022