Measurement of Sxa profiles

We extracted the shape of 31 Sxa from SEM images (see Methods). Since the Sxa are axisymmetric, we describe a Sxa’s shape using its “profile”, which is a set of points , i = 1 … 250 shown in Fig. 1(b,c). We measured the length, L m , and maximum cross-sectional radius, R m , of each Sxa from its profile (see Supplementary Section Details of Sxa profile measurements). By plotting the dimensionless profiles, , for i = 1 … 250, we see that the general nature of the taper appears uniform across different Sxa. To make a more quantitative comparison of the Sxas’ tapers we compute the aspect ratio, α m = L m /2R m , for each Sxa. The values of α m are plotted in Fig. 1(d). The mean and standard deviation of α m are 53.6 and 8.7, respectively. The small scatter of α m further supports our viewpoint that the tapered shape is uniform across different Sxa.

Mechanical testing of Sxa

The Sxa are primarily composed of silica35, which is a well-characterized ceramic material that behaves in a linear elastic fashion and fails through brittle fracture. Cursory inspection of the surfaces of fractured Sxa (see Fig. 2(b)) suggests that they are essentially homogeneous silica rods. However, spicules also possess a proteinaceous scaffold within their silica36,37. In some related species this protein forms distinct layers, which may affect the deformation and failure behavior of the spicules13,35,38,39. While T. aurantia’s Sxa do not contain separate layers of protein and silica, the influence of any underlying protein scaffold on their elastic behavior is unknown. Furthermore, the composition of the silica itself varies spatially within the Sxa40. To ascertain the effect of any potential elastic inhomogeneities on a Sxa’s deformation behavior, we performed three-point bending tests on 30 Sxa using a custom-built flexure device. Briefly, the Sxa were suspended over a trench with vertical, parallel walls and were indented by a cantilever that also acts as a force sensor. The details of our flexure device and three-point bending tests will be published elsewhere.

The magnitude of the transverse force, F, and deflection of the Sxa’s axis in the y-direction at midspan, w 0 (see Fig. 2(c)), were recorded until the Sxa failed. The w 0 –F data for the Sxa are shown in Fig. 2(a). The failure of every Sxa we tested was defined by a single fracture event. The fracture events are marked with red points in Fig. 2(a). The w 0 –F response of every Sxa was linear until failure. This observation indicates that the Sxa’s mechanical behavior is linear elastic until failure.

We compared a Sxa’s deformed shape during a bending test to that predicted by Euler-Bernoulli theory for an elastically homogeneous, tapered beam41. Euler-Bernoulli (EB) theory is a highly successful structural mechanics theory used for modeling the deformation of slender, linear elastic structures that primarily deform through bending. For details of how the shape predicted by EB theory was computed, see Supplementary Section Deflection of a tapered beam in three-point bending predicted by Euler-Bernoulli theory. The displacement of a Sxa’s axis in the y-direction was measured from images taken during the bending test (see Fig. 2(c)). A representative comparison of these measured displacements with those predicted by EB theory is shown in Fig. 2(d). The measurements and the theoretical predictions match very well for 27 of the 30 Sxa. This supports that the Sxa’s behavior is linear elastic and shows that a Sxa is elastically homogeneous along its length. Furthermore, it shows that a Sxa’s deformation can be described by an EB theory for an elastically homogeneous, tapered, axially symmetric beam.

Computational mechanics calculations

Being embedded within the sponge, a Sxa likely experiences a complex distribution of tractions along its length. However, using computational mechanics calculations we found that due to the Sxas’ arrangement within the sponge and the large mismatch between the compliance of the Sxa and the spongin, the tractions are localized at the ends of the Sxa (see Supplementary Section Computational mechanics model of a Sxa in its RoC). Thus, the most appropriate structural mechanics model based on EB theory would be a simply supported column, which is described by equations (1,2,3).

The Sxa are not uniformly scattered throughout the sponge’s body, rather they are grouped in bundles that extend radially from the sponge’s center to its outer surface (see Fig. 3(a))26. The Sxa are aligned along the bundles’ lengths and are staggered with respect to each other (see Fig. 3(b)). The bundles are 220–490 μm thick26 and a bundle’s cross-section contains approximately 50 Sxa27. From the average bundle thickness, number of Sxa per bundle, and Sxa diameter, we estimate the distance between the axes of neighboring Sxa to be ≈45 μm (see Supplementary Section Estimation of the distance between adjacent Sxa in a bundle). Thus, the Sxa within a bundle are separated from each other by a small amount, ≈8 μm, of spongin.

External forces acting on the sponge are transmitted by the spongin to the Sxa as tractions on their surfaces. To determine the distribution of these tractions, we performed a stress analysis on a continuum mechanics model of an individual Sxa embedded in a cylindrical section of spongin. We refer to this cylinder as a Sxa’s region of confinement (RoC) (see Fig. 3(c,d)). The diameter of the RoC is equal to the distance between neighboring Sxa in a bundle.

We model the spongin in the RoC as an isotropic, linear elastic solid with Young’s modulus and Poisson’s ratio of 600 KPa and 0, respectively. These values correspond to measurements of the mechanical properties of spongin in a related species28. Furthermore, measurements of the Young’s modulus of spicules from a related species42 indicate that the silica is between four and five orders of magnitude stiffer than the spongin. We will present Young’s modulus measurements of the Sxa from our own work in a future paper. Motivated by this large difference in stiffnesses, we model a Sxa as a rigid inclusion whose surface is bonded to the spongin in its RoC. We assume that external forces act normal to the sponge’s surface and result in axial compressive stresses in the Sxa bundles. Therefore, we apply compressive tractions to the ends of the RoC (see Fig. 3(d)). Since the spongin in a RoC is also connected to the spongin in the RoCs of neighboring Sxa, we constrain points on the lateral surface of the RoC from moving in the radial direction. Further details about this model can be found in Supplementary Section Computational mechanics model of a Sxa in its RoC.

We computed the stress field in the spongin using finite element procedures (see Fig. 3(e))43. We found that for a wide range of traction distributions applied to the ends of the RoC, the axial force per unit length acting on the Sxa is always localized on the Sxa’s ends (see Fig. 3(e) and Supplementary Section Computational mechanics model of a Sxa in its RoC). This localized force distribution contrasts with that predicted for an ellipsoidal inclusion embedded in a linear elastic solid subjected to far-field compressive stress. Specifically, a celebrated elasticity solution by Eshelby44 predicts that the axial force per unit length will vary in a piecewise affine fashion along an ellipsoidal inclusion. It is not necessary, however, for this result to hold true for non-ellipsoidal inclusions. Thus, our numerical results do not contradict Eshelby’s solution. In fact, they are consistent with results from computational models of short fiber reinforced composites45,46, full-field elasticity solutions for rigid line inclusions47, and photoelasticity experiments on line-like inclusions48. Based on the insight gained from our computational mechanics calculations, we modeled the effect of the spongin by replacing the tractions applied to the ends of the RoC with opposing point forces, ±P M at the Sxa’s ends (see Fig. 3(f)).

The structural mechanics model for the Sxa

Initially a Sxa behaves like a column with two free ends, which is unstable when subjected to the the axial forces ±P M . Even if these forces are aligned with the Sxa’s axis, small perturbations in the configuration will inevitably cause the Sxa to rotate about one of its transverse axes. However, after rotating by only a small amount (≈1.3°), the proximity of neighboring Sxa in the bundle will prevent further rotation (see Fig. 3(f)). Due to the spongin’s large compliance, it is unlikely that this small rotation will substantially change the stress state in the spongin and consequently the traction distribution on the Sxa’s surface. However, there will be non-negligible reaction forces, ±P N , at the points where a Sxa is restrained by its neighbors. The net force at a Sxa’s end, P, which includes contributions from P M and P N , must act in the direction of the Sxa’s axis (see Fig. 3(f)). This is a consequence of static equilibrium and can be deduced using a free body diagram.

Thus, a Sxa can be modeled using the EB theory in which the column’s ends are subjected to compressive, axial forces and cannot move in the direction perpendicular to the column’s axis. We refer to this model as a simply supported column (see Fig. 3(g)). In this model, the transverse deflection, w, is governed by the differential equation

for all z ∈ (0, L), and boundary conditions

where P, E, L and I are the magnitude of P, the column’s Young’s modulus, length and second moment of area, respectively. Based on the results of Sections Measurement of Sxa profiles and Mechanical testing of Sxa we take E to be constant and I(z) = πr(z)4/4, where r(z) is the radius of the Sxa’s cross-section—i.e. its profile.

Comparison with the Clausen profile

The buckling strength of a simply supported column is the smallest P for which there exists a solution to equations (1,2,3) other than w = 0 for all z ∈ [0, L]. For an elastically homogeneous column, the buckling strength can be modulated by varying I, or in this case r, along the column’s length31. Our hypothesis would gain support if the profile of the simply supported column with the greatest buckling strength resembled the measured profiles of the Sxa.

The profile that maximizes a simply supported column’s buckling strength for a given length, L, and volume, V, was first sought by Lagrange in the late 1700s49. The correct solution, however, was discovered in 185150, and an accessible proof that it is in fact optimal was given in 196234. This optimal profile, which we refer to as the Clausen profile, is given by

where ρ = r/L and ζ = z/L are the dimensionless radial and axial coordinates, respectively, and θ is a parameter that lies between 0 and π33,34. The parameter α = (3πL3/16V)1/2 is a measure of the column’s aspect ratio. We refer to a column whose taper is described by the Clausen profile as a Clausen column (see Fig. 4(A)).

To test our hypothesis, we compared the Clausen profile to the Sxa profiles. We did this by fitting equations (4) and (5) to each Sxa profile in the least-squares sense by varying the parameter α (see Supplementary Section Fitting profiles to the Sxas’ shape). The best fit Clausen profile for a representative Sxa is shown in Fig. 4(B). We also fit three other prototypical profiles; a semiellipse, an isosceles triangle and a constant to the Sxa profiles (see Supplementary Section Fitting profiles to the Sxas’ shape and Fig. 4(A)). We use the sum of squared residuals for a fitted profile, mSSR, to indicate how well that profile describes a Sxa’s shape. The medians, means and standard deviations of each profile’s mSSR are shown in Table 1, from which we see that the Clausen profile has the lowest mean and median mSSR. Furthermore, a two-sided Wilcoxon signed rank test indicates that the median mSSR for the Clausen profile differs from that of the semiellipse profile at the 1% significance level (p = 0.0002). Thus, using the median mSSR as a metric, we conclude that the Clausen profile describes the Sxas’ tapers the best out of the different profiles that we considered. Further investigation using another comparison criterion also supports this conclusion (see Supplementary Section Additional profile comparison using the Akaike information criterion for details).

Table 1 mSSR of the candidate profiles (N = 31). Full size table

Direct estimates of the Sxas’ buckling strengths

The fact that the Clausen profile describes the measured Sxa profiles the best among the prototypical tapered profiles that we considered gives strength to our hypothesis. However, it is still possible that there may exist some other profile, which corresponds to an alternate hypothesis, that describes the Sxa’s taper even better than the Clausen profile. If such a profile exists, would our hypothesis remain viable?

To answer this question, we numerically estimated the Sxas’ buckling strengths, P s , using the measured profiles and our structural mechanics model. Briefly, we computed a Sxa’s second moment of area I(z) = πr(z)4/4 from its profile, r(z), and used the Rayleigh-Ritz method51 to find an approximate value for the smallest P for which there exists a solution to equations (1,2,3) other than w = 0. We computed P s for each of the 31 Sxa whose profiles we measured in Section Measurement of Sxa profiles and compared it to the buckling strength P c = πEV2/(4L4) of the equivalent cylinder—i.e., the cylinder with the same length, volume, and elastic properties (see Fig. 5). Taken as a group, we found that the median buckling strength of the Sxa is 13.4% greater than that of their equivalent cylinders. Furthermore, some Sxa achieve values of (P s − P c )/P c as large as 0.3 which is close to the enhancement of 0.33 provided by the Clausen column33 (see Fig. 5).

So, even if there existed a profile that better resembled the Sxas’ tapers, the fact still remains that the Sxas’ tapers substantially enhance their buckling strengths. Therefore, even if there existed a better matching profile based on an alternate hypothesis, the support for our hypothesis would still remain strong. Such a scenario would only mean that the Sxa serve more than one function.