Key Signatures of Magnetofossils Elucidated by Mutant Magnetotactic Bacteria and Micromagnetic Calculations

Abstract Magnetotactic bacteria (MTB) produce single‐stranded or multi‐stranded chains of magnetic nanoparticles that contribute to the magnetization of sediments and rocks. Their magnetic fingerprint can be detected in ancient geological samples and serve as a unique biosignature of microbial life. However, some fossilized assemblages bear contradictory signatures pointing to magnetic components that have distinct origin(s). Here, using micromagnetic simulations and mutant MTB producing looped magnetosome chains, we demonstrate that the observed magnetofossil fingerprints are produced by a mixture of single‐stranded and multi‐stranded chains, and that diagenetically induced chain collapse, if occurring, must preserve the strong uniaxial anisotropy of native chains. This anisotropy is the key factor for distinguishing magnetofossils from other populations of natural magnetite particles, including those with similar individual crystal characteristics. Furthermore, the detailed properties of magnetofossil signatures depend on the proportion of equant and elongated magnetosomes, as well as on the relative abundances of single‐stranded and multi‐stranded chains. This work has important paleoclimatic, paleontological, and phylogenetic implications, as it provides reference data to differentiate distinct MTB lineages according to their chain and magnetosome morphologies, which will enable the tracking of the evolution of some of the most ancient biomineralizing organisms in a time‐resolved manner. It also enables a more accurate discrimination of different sources of magnetite particles, which is pivotal for gaining better environmental and relative paleointensity reconstructions from sedimentary records.


Supporting
: Parameter values used for chain simulation. Statistical distributions ( =  Fischer-Von Mises, =  Beta, =  Normal, =  Uniform) indicate that the corresponding parameter is a random realization of the given distribution, e.g.

Parameter
Single-stranded chains Native doublestranded chains

Supporting Movies
Movie S1: Simulated zero-field reversal of a single-stranded chain. The chain is made of 15 SD particles without shape and magnetocrystalline anisotropy.
Movie S2: Simulated zero-field transition between the high-and low-moment states of a ring. The ring is made of 20 SD particles without shape and magnetocrystalline anisotropy.
Movie S3: Simulated zero-field transition between the high-and low-moment states of a doublestranded chain. The strands contain 11 and 10 SD particles, respectively, without shape and magnetocrystalline anisotropy.
Movie S4: Simulated zero-field transition between the high-and low-moment states of a foldcollapsed chain. The two strands contain 11 and 10 SD particles, respectively, without shape and magnetocrystalline anisotropy.

Glossary of magnetic terms
Anhysteretic remanent magnetization (ARM): Remanent magnetization obtained by placing the specimen in a decaying alternating magnetic field superimposed to a small direct current (DC) field. ARM is strongly selective towards non-interacting SD particles and magnetic systems that behave equivalently, such as chains. This selectivity is similar, but not identical, to that of the central ridge in FORC diagrams. The anhysteretic susceptibility, a d c ARM H = / is defined as the ARM normalized by the DC field dc H (in A/m) used during acquisition.

Applied field (B):
In the FORC measurement protocol, it is the magnetic field applied during magnetic measurements.
ARM ratio (χ a /IRM): The ARM ratio a IRM / , where IRM is the isothermal remanent magnetization acquired in a DC field equal to the maximum amplitude of the AC field used to obtain the ARM, is a magnetic grain size indicator. Non-interacting SD particles and intact magnetosome chains are characterized by a mm A >1 / .
Bias field (B u ): Vertical coordinate of the FORC diagram. It has a direct physical interpretation only in the Preisach-Néel model of interacting SD particles, where it corresponds to the internal interaction field acting on individual particles.

Biogenic soft (BS):
Typically used to describe the magnetic signature of conventional (i.e., not giant) magnetofossils made of equidimensional crystals. The coercivity distribution is characterized by a median field of ∼45 mT and a mean dispersion parameter of ∼0.18. The magnetic anisotropy required to explain the range of coercivities covered by this component (∼10-100 mT) originates largely from the chain structure, since individual crystals are almost isotropic.
Biogenic hard (BH): Typically used to describe the magnetic signature of conventional (i.e., not giant) magnetofossils made of elongated crystals. The coercivity distribution is characterized by a median field of ∼70 mT and a mean dispersion parameter of ∼0.1. The magnetic anisotropy required to explain the range of coercivities covered by this component (∼30-140 mT) originates from the chain structure and from magnetosome elongation.

Blocking volume:
The minimum volume that a SD magnetic particle must have to maintain a stable magnetization at a given temperature.
Boltzmann factor: The ratio between the energy barrier E ∆ that needs to be overcome for a transition between two given magnetic states at the absolute temperature T, and the energy Central ridge: Name given to a distinct feature of FORC diagrams, which consists of a narrow ridge extending along, or very close to the horizontal axis defined by u B = 0 . It is a distinctive feature of isolated magnetic systems possessing only few magnetic states. Notable examples are noninteracting SD and vortex particles, and systems of particles that mimic these domain states, such as magnetosome chains. Because of its quasi-unidimensional nature, the central ridge is separable from other contributions to the FORC diagram. In sediments, the central ridge is often uniquely associated with secondary SD magnetic particles and magnetofossils, since primary production of non-interacting SD magnetite requires very fast cooling and magnetic particles with larger sizes are rarely limited to the narrow range of vortex particles.

Coercivity (B c ):
In the context of magnetic hysteresis, it is defined as the absolute value of the applied field for which the magnetization is zero during the measurements of a major hysteresis loop; that is, for the ascending branch M − , and c M B + − = 0 ( ) for the descending branch M + . It is a common measure of the "magnetic hardness" of a specimen. In case of magnetic particle assemblages, c B depends strongly on the domain state and on magnetic interactions. In SD particles, c B is a measure of their magnetic anisotropy energy, with larger values in case of strong particle elongation or strong magnetocrystalline anisotropy. In MD particles, c B depends essentially on how strongly pinned are domain walls by crystal defects. c B is also used to denote the horizontal coordinate of FORC diagrams, where it is interpreted as the coercivity of rectangular hysterons associated with the FORC function.

Coercivity of remanence (B cr ):
It is defined as the DC field that must be applied in the opposed direction to cancel the saturation remanent magnetization rs M , obtaining a zero remanent magnetization. In FORC measurements, it coincides with the reversal field r B of the curve that yields It is a measure of the "magnetic hardness" of the remanent magnetization of a specimen, having thus a similar meaning as c B . The so-called coercivity ratio cr c B B / is always ≥1 if FORC curves are limited by the major hysteresis loop, as usually observed in geologic materials. The coercivity ratio of particle assemblages depends on the domain state of particles, with values close to 1 for SD assemblages and values 1  for MD particles.
Coercivity or switching field distribution: The strict definition refers to the statistical distribution sw f B ( ) of switching fields in systems possessing large numbers of magnetic states (e.g., those of individual particles. Because switching fields are rarely directly measurable, a coercivity or switching field distribution is empirically identified with any field-dependent function f B ( ), whose integral coincides with a certain type of total magnetization (for instance the saturation remanent magnetization), or total magnetization of a subset of particles in a specimen. Vice versa, coercivity distributions can be defined as the first derivative of a so-called magnetization curve M B ( ) obtained with a specific measurement protocol. In these cases, the argument of the function, B , is a generic field that is not strictly identifiable with the coercive field of hysteresis, nor with a switching field. Coercivity distributions commonly used in the literature include those associated with IRM acquisition, DC demagnetization, AF demagnetization of IRM or ARM, and the central ridge.
Coercivity distribution (AF demagnetization, f af ): The derivative af / M B − of the AF demagnetization curve of a remanent magnetization (ARM or IRM), obtained in successive demagnetization steps with maximum AF amplitude af B . The integral of af f yields the initial remanent magnetization. Therefore, af f can be identified with the coercivity distribution of magnetic minerals carrying the specific type of remanent magnetization being demagnetized.
Coercivity distribution (IRM acquisition, f irm ): The derivative / M B of the IRM acquisition curve M B ( ) with respect to the IRM acquisition field r B . The integral of irm f yields the saturation remanent magnetization rs M . Therefore, irm f can be identified with the coercivity distribution of magnetic minerals carrying a saturation remanent magnetization. It is similar, but not identical, to the coercivity distribution obtained from DC demagnetization, the difference being that IRM acquisition starts from a fully demagnetized state.
Coercivity distribution (central ridge, f cr ): It is obtained by integrating the central ridge contribution cr c u ( , ) B B to the FORC function over u B . The integral of cr f over c B yields the total contribution of irreversible magnetization processes to the central ridge. This is a fraction of all irreversible contributions to the hysteresis loop, so that cr irr f f for any positive field, with cr irr f f = for non-interacting uniaxial SD particles.

Direct current demagnetization (DCD):
The procedure used to null the saturation remanence rs M of a specimen, consisting in the application and removal of a field cr B , called coercivity of remanence, in the opposed direction to rs M . Because cr B is not known a priori, increasingly large 'reversal' fields r B are applied and removed, and the resulting remanent magnetization r r ( ) M B is measured after each step.

Energy barrier (ΔE):
The minimum increase in the free energy E of a magnetic system required for a transition between two magnetic states while external conditions (e.g., temperature, applied field) are kept constant. If ΔE is sufficiently small, typically ~20 times the energy B k T of thermal fluctuations (see Boltzmann factor), the transition occurs spontaneously within a certain time.
Ferrimagnetism: Describes a class of magnetic materials that includes ferromagnetism (parallel atomic spins), and strict ferrimagnetism (two antiparallel spin lattices with unequal magnetic moments). The common characteristics of these materials is that they possess a spontaneous magnetization, defined as the magnetization of any microscopic (sub-domain) volume of the material. The saturation magnetization of ferrimagnetic particles is equal to the spontaneous magnetization of the material of which they are made, up to fields that are sufficiently large to overcome the given ferrimagnetic order.
First-order reversal curves (FORCs): Series of partial hysteresis curves directly originating (order 1) from the major hysteresis (order 0). In the most commonly used measurement protocol, a FORC curve is measured as follows: (1) a positive saturation field s B is applied, which resets the specimen to a positive saturation state, (2) the field is decreased until a predefined, so-called reversal field )/ represents the contribution of a fictive rectangular hysteresis loop (called hysteron) with coercive field c B and subjected to a fixed bias field u B . In the case of ensembles of interacting SD particle assemblages, c B and u B are identifiable with the coercivity and the interaction field of each particle. The central ridge is the limit case with no interactions, in which case is condensed along u B = 0 . In other systems, reflects reversible and irreversible magnetization processes responsible for the differences between consecutive FORCs, being no longer associable with the coercivity and interaction field distribution of discrete particles.
Fluctuation field (B f ): A fictive field defined as the difference between the fields required to trigger a given transition between magnetic states without and with the effect of thermal fluctuations (random excursions of the magnetization caused by thermal agitation). These fluctuations are then equivalent to a fluctuating external field with amplitude f B . Magnetic particles become effectively superparamagnetic when f B is as large as the field required for the transition without fluctuations.
Flux closure (FC): Magnetic configuration of a single particle, a group of particles, or any other ferrimagnetic system, which produces a closed internal magnetic flux i B H M = + . FC states tend to minimize or nullify the net magnetic moment of the system. The vortex state (VO) of single particles is an example of flux closure configuration. Magnetic domains in soft materials also tend to form flux closure configurations in sufficiently small external fields.  Isothermal remanent magnetization (IRM): Remanent magnetization obtained by applying and subsequently removing a direct current (DC) field to a previously demagnetized specimen, while the temperature is kept constant. Only magnetic states that are irreversibly changed by the applied field contribute to the IRM, while the other states maintain the randomly oriented magnetic moments of the original demagnetized state. The saturation remanent magnetization rs M is a particular IRM obtained with a saturating field, that his, a field that is sufficiently strong to irreversibly change all magnetic states.

Magnetic hysteresis:
The phenomenon by which the magnetization of a specimen cycled between two fields B 1 and B 2 depends on the direction of the field sweep (from B 1 to B 2 or vice-versa). Magnetic hysteresis arises from irreversible changes of magnetic states in the specimen. The major hysteresis loop is obtained with Magnetic reversibility and irreversibility: Ferrimagnetic materials can respond in a reversible or an irreversible manner to changes of the external conditions (e.g., the applied field). A small change of the external conditions produces an irreversible magnetization change if the magnetization does not revert to the original configuration after restoring the original conditions. Total magnetization changes are usually a mixture of reversible and irreversible processes. Magnetic hysteresis is a wellknown example of irreversibility: sweeping the applied field between two extreme values produces magnetization changes along two distinct paths, called the ascending and descending branches of the hysteresis loop. Magnetization changes are always reversible in field amplitudes larger than the saturation field. FORC diagrams are two-dimensional representations of mixed reversible and irreversible magnetization changes. Purely reversible systems, such as above the saturation field or in SP particle assemblages, do not contribute to the FORC function.

Magnetic shape anisotropy:
The tendency of the spontaneous magnetization of a ferrimagnetic particle to align with the long axis (e.g., in an iron needle).

Magnetic state (or configuration):
Describes a specific magnetization configuration inside a ferrimagnetic particle, a group of ferrimagnetic particles, or any other ferrimagnetic system, which is stable against small field changes: application and removal of a small additional external field brings the system back to its original configuration. A transition to another magnetic state occurs when the magnetization change induced by external conditions (e.g., the applied field) becomes irreversible: in this case, restoring the original conditions does not resume the original magnetic state. Examples of magnetic state transitions are those between the two single-domain configurations (parallel and antiparallel) of a uniaxial particle, and between single-domain and vortex configurations. Magnetic state transitions are usually, but not always, characterized by a sudden change of the net magnetic moment. The appearance and disappearance of a magnetic state are often called nucleation and annihilation, respectively.
Magnetic viscosity or thermal relaxation: Describes the delayed response of a magnetic system to any change of the external conditions (typically the applied field). Magnetic viscosity is mainly observed in small SD particles close to the limit where they become superparamagnetic, and, sometimes, in MD particles. The manifestation of viscosity depends on the timescale: the same particles can be stable during laboratory measurements, and viscous over geologic times.

Magnetic vortex (VO):
Describes a magnetization whose direction is described by a vortex vector field. In the ideal vortex state, atomic spins are parallel to all free surfaces of the enclosing volume, so that no stray field is produced, and the net magnetic moment is carried only by spins close to the vortex axis (the so-called vortex tube). Vortex states occur in material with weak magnetocrystalline anisotropy, such as magnetite, over a size range between the upper SD stability limit and the lower limit for the nucleation of domain walls (MD).

Magnetization (M):
In a magnetic material, it is defined as the total magnetic moment of atomic or molecular spins in an infinitesimal volume element, divided by volume. In the context of magnetic measurements, it is defined as the net magnetic moment of the specimen, normalized by the specimen volume, mass, or area.

Magnetocrystalline anisotropy:
The tendency of the spontaneous magnetization of a ferrimagnetic material to align with preferred crystalline axes, called easy axes of magnetization.

Magnetostatic interaction:
The influence of the stray field produced by one magnetic particle on another magnetic particle. Interacting magnetic particles are subjected to an internal, so-called interaction field produced by the particles themselves, in addition to the external field applied, for instance, during FORC measurements. Weak magnetostatic interactions, typically occurring between particles separated by several diameters, introduce a small bias to the fields required to promote transitions between magnetic states. In case of random particle assemblages, the internal biasing field is a locally defined random variate that produces a random vertical offset of single particle contributions to the FORC diagram. Therefore, the FORC diagram of weakly interacting particles is a vertically blurred version of the FORC diagram of the same isolated particles without interactions. Strong magnetostatic interactions between closely spaced particles, such as those occurring between magnetosomes in a chain, can produce a strong magnetic coupling that replaces individual magnetic states with a collective state. For instance, isolated magnetosomes possess individual SD behaviors with different switching fields (controlled mainly by crystal elongation), while magnetosomes in a linear chain display a collective SD behavior with a single switching field.

Multi-domain (MD):
Domain state used to describe magnetic crystals containing several discrete, homogeneously magnetized domains separated by domain walls. In the absence of external fields, magnetic domains tend to be magnetized in a manner that cancels all stray fields, compatibly with size, shape, and magnetocrystalline anisotropy. In so-called magnetically soft materials, such as magnetite, domain walls can easily move when a small external field is applied. Therefore, the magnetization of magnetically soft MD crystals is typically much less stable than that of SD or VO particles and depends essentially on pinning by crystal defects. The lower MD size limit for magnetite is a few microns, depending on shape and other factors.

Remanent magnetization (M r ):
The magnetization of a specimen in a null external field. All ferrimagnetic materials except SP particles can hold a non-zero remanent magnetization that depends on their previous hystory.

Reversal field (B r ):
The field at which a FORC measurement starts. Application of r B to a positively saturated specimen brings the specimen in a mixed magnetic state that reverts to positive saturation while a FORC is measured.

Saturation field (B s ):
The field amplitude above which the major hysteresis loop becomes completely closed. Application of s B converts all magnetic states into a SD-like configuration. It is therefore applied before each FORC measurement to reset the specimen to the same initial magnetic state.

Saturation magnetization (M s ):
The magnetization acquired in a magnetic field that is sufficiently large to align all atomic or molecular spin moments. In a ferrimagnetic material, s M approaches asymptotically the spontaneous magnetization of the material. It is the only type of specimen magnetization that is simply proportional to the amount of magnetic material that is contained in it and is therefore a good indicator of magnetic concentration.

Saturation remanent magnetization (M rs ):
It is defined as the specimen magnetization in a null external field, which remains after removing a saturation field. It coincides with the zero-field magnetization of the upper and lower branches of the major hysteresis loop. It is usually also the largest remanent magnetization of a specimen in a null field. The ratio rs s / M M depends on the domain state of particles and the magnetic interactions existing between them. Assemblages of non-interacting, randomly oriented, uniaxial SD particles are characterized by Single domain (SD): Domain state describes crystals that are sufficiently small to contain only one magnetic domain. SD particles possess the maximum possible magnetic moment, equal to the product of the spontaneous magnetization and the particle volume. Every particle becomes SD in a saturating external field, regardless of size. Therefore, the domain state of a particle is usually referred to its remanent magnetization.

Stable single domain (SD):
Describes SD particles whose magnetization is stable over the experimental time range (seconds to hours for magnetic measurements, millions of years in paleomagnetism). The definition is extended to all systems with same magnetic behavior as SD particles, such as magnetosome chains.

Superparamagnetism (SP):
Describes SD particles whose magnetic anisotropy energy is small enough for thermal fluctuations to overcome the energy barrier required to rotate the magnetization between stable orientations (two for a uniaxial particle). Accordingly, superparamagnetic particles do not possess a stable magnetization and do not exhibit an open hysteresis loop. The transition from stable SD to SP depends on temperature and measurement time. The lower stability limit for SD magnetite particles is ∼20 nm at room temperature.

Switching field (B sw ):
The field that is required to nucleate a transition between magnetic states. It is often referred to uniaxial SD particles; in this case, it is the field required to switch between two opposed SD states (parallel and antiparallel to the easy magnetic axis).
Uniaxial Non-interacting single domain (UNISD): UNISD particles represent the simplest model of SD magnetic particles, which has been used to explain the central ridge signature of FORC diagram, as well as other magnetic properties of magnetofossil-rich sediments.

Supporting information: AMB-1 mutant
Deletion plasmid construction and generation of the AMB-1 mutant ΔmamJΔlimJ ΔMIS The MIS deletion construct was created by Gibson assembly. The upstream (AB recombination fragment, 904 bp) and downstream (CD recombination fragment, 842 bp) regions of the islet were amplified by PCR from the AMB-1 genomic DNA with the following primer sets: jwP1A (5'-gaattcctgcagcccgggggatccACTAGTgttcccctccatcacctatac-3') and jwP1B (5'-CCCATCCACTAAATTTAAATAagcggcggtatagcccatg-3') for the AB recombination fragment, jwP1C (5'-TATTTAAATTTAGTGGATGGGctattcgaaccgcct gctc-3') and jwP1D (5'-caccgcggtggcggccgctctagaACTAGTgatagcgagaaccgtcatac-3') for the CD recombination fragment. The Gibson assembly was used to clone AB and CD fragments into the SpeI site (underlined in the above primers) of pAK0 (Komeili et al., 2006), which is a suicide plasmid carrying a kanamycin resistance cassette and a counterselectable sacB gene, to generate pAK1121. A spacer (highlighted bold in the above primers) at the end of AB and the beginning of CD fragments was used to connect the AB and CD fragments during Gibson assembly. ΔmamJΔlimJΔMIS was generated by a two-step recombination method (Murat et al., 2010). The deletion plasmid pAK1121 was conjugated into ΔmamJΔlimJ AMB-1 strain (Draper et al., 2011) by using WM3064 as the donor strain. For the first recombination step, the colonies were selected on agar plates containing kanamycin (medium composition given by the ATCC under AMB-1 reference), to obtain the recombinants that have the deletion plasmid integrated into the genomic DNA. For the second recombination step, cells were grown in 10 ml of AMB-1 growth medium without kanamycin, and then counterselected on agar plates containing 2% sucrose, to obtain recombinants that lost the integrated plasmid. These sucrose-resistant colonies were then screened for the absence of MIS by culture PCR with the following primers: jwP2F (5'-catcaccatgaccctgaccg-3') and jwP2R (5'-gacgttttgaaggggctggac-3').

Supporting information: micromagnetic calculations
This section explains in detail how synthetic magnetosome chains are constructed, and the methods used for simulating thermally activated transitions between magnetic states and whole FORC measurements. Mathematical symbols used for simulation parameters and other parameters are listed in Table S2 and Table S3 at the end, respectively. Parameter values used for the simulations are listed in Table S1.

Magnetosome properties
Micromagnetic simulation of magnetofossils are controlled by two groups of parameters describing the properties of individual crystals on the one hand, and their arrangement in intact or collapsed chains on the other. Relevant crystal properties include size, shape, and, to a minor extent, crystallographic orientation with respect to shape.

Sizes and shapes
Magnetosomes occur in a variety of complex, strain-specific crystal shapes (Mann et al., 1984(Mann et al., , 1987Sparks et al., 1990;Meldrum et al., 1993a,b;Lefèvre et al., 2011;Clemett et al., 2002). Because magnetic properties of single-domain (SD) magnetite particles are mainly controlled by their elongation, magnetofossils can be subdivided into three main categories for magnetic modeling purposes: (1) equant or isometric crystals, such as cuboctahedra and octahedra, which do not possess a systematic elongation, (2) prismatic crystals, which are moderately elongated, and (3) tooth-or bullet-shaped crystals, which are highly elongated in their mature stage. These three categories form three distinct clusters in a so-called Butler-Banerjee diagram of short-to-long axis ratio vs. length (Fig. S1a). The size and shape distribution within each cluster is determined by different bacterial strains and natural variations within each strain, related to the growth stage or environmental conditions. The apparent elongation of cuboctahedral and octahedral crystals is due in part to projection effects of the TEM micrographs from which shape and size data are retrieved, and in part to random deviation from a perfectly isometric shape.
Magnetofossil sizes are mostly comprised within the limits required for intact magnetosome chain to possess a stable SD-like remanent magnetization, with the magnetic moment of each crystal being parallel to the chain axes and not switchable by thermal activations at room temperature (Muxworthy and Williams, 2006;Newell, 2009). The stability range for isolated crystals is much smaller (Butler and Banerjee, 1975), with a non-negligible proportion of isometric crystals being either too small or too large to carry a stable single-domain remanence without the stabilizing effect of magnetostatic interactions within a chain. Unlike the crystal length, widths within each shape category are almost independent of elongation. Therefore, individual crystal shapes are effectively described by independent width and axis ratio distributions (Fig. S1b-d). These distributions are well approximated by a rescaled versions of the Beta function where is the Gamma function and , 1 2 two shape parameters.

Shape axes
The shape of simulated magnetosomes is controlled by three orthogonal axes , , a a a  Magnetosomes in real chains feature small, random variations of their axis orientation and length, relative to the ideal case. In order to reproduce such variations, the orientation of the a 3 -axes is generated with a Fisher-Von Mises distribution where is the angle between a 3 and the chain axis, and s a so-called concentration parameter governing the dispersion of around the chain axis. Large values of s correspond to small deviations of a 3 from the chain axis, with the standard deviation of being proportional to s −1 for s 1  . The other two axes are oriented randomly.
The axes lengths of the n-th magnetosome are given by where n x , n z are random Gaussian variates with standard deviation s 1  , limited to ±0.8, and t n < 0 1 ( ) is a size tapering function, which depends on the position of the n-th magnetosome with respect to the closest chain extremity. Cylindrical symmetry is imposed through , , i i a a = 1 2 , since the axes are used merely to simulate a uniaxial magnetic anisotropy due to the crystal elongation. The magnetosome shapes are thus controlled by the following set of chainspecific parameters: the crystal width s 0 and elongation e 0 of mature crystals, the dispersion s of long axes orientations with respect to the chain axis, the standard deviation s of random size variations, and the tapering function t . Shape axes are constructed by applying the following rules to each crystal: (1) A set of orthogonal axes , , x y z ( ) is generated, with z being parallel to the chain axis.
(2) The axes are rotated about z by a random angle with uniform distribution between 0 and 2 .
(3) The axes are rotated about a random axis perpendicular to z , by an angle generated by a random realization of the Fisher-Von Mises distribution in eq. (2).

Crystallographic orientations
Crystal orientation is identified with the 100 [ ] magnetite axes. In the case of ideal cuboctahedral or prismatic magnetosomes, it is assumed that one of the 111 [ ] axes is parallel to the chain axis, since this is the typical orientation observed for these shapes (Simpson et al., 2005;Pósfai et al., 2013;Körnig et al., 2014). The other 111 [ ] axes appear to be randomly oriented (Simpson et al., 2005;Orue et al., 2018). The crystal orientation of simulated magnetosomes is a randomized version of the ideal case, with the 111 [ ] axis closest to the chain axis being governed by a Fisher-Von Mises distribution (eq. 2) with concentration parameter c . The crystal axes are generated with the following steps: ( (2) The axes are rotated about z by a random angle with uniform distribution between 0 and 2 .
(3) In case of elongated magnetosomes, the axes are rotated about a random axis perpendicular to z in the exact same manner as for the shape axes (i.e., using the same rotation axis and rotation angle). This rotation generates the whole crystal misalignment with respect to the chain axis.
(4) The axes are rotated about a random axis perpendicular to z , by an angle generated by a random realization of a Fisher-Von Mises distribution with chain-specific concentration parameter c . This rotation generates an additional misalignment of the crystalline axes with respect to the shape axes.

Magnetosome shapes
Magnetosome shape determines the uniaxial magnetic anisotropy used in micromagnetic simulations. This anisotropy is expressed by the so-called demagnetizing tensor D . The calculation of D is performed through an intermediate step where the real magnetosome shape is approximated by an idealized shape. The idealized shape of equidimensional (cuboctahedral) magnetosomes is a rotation ellipsoid with axes lengths a a = 1 2 and a 3 . Elongated prismatic magnetosomes are represented by chamfered circular cylinders with diameter a a = 1 2 and length a 3 . The circular cylinder approximates a hexagonal prism with 110 { } side faces and 111 { } top/bottom faces, as typically seen in MV-1 magnetosomes (Weyland et al., 2006).

Additional 100 { } and 111
{ } faces eliminating right angles between side and top/bottom faces (Clemett et al., 2002;Pósfai et al., 2013) are simulated by 45° chamfers of width s w 0 , such that the diameter of top/bottom faces is reduced by s w 0 2 (Fig. S9). The chamfer width parameter w is subjected to the same randomization as the shape axes length, so that actual chamfer widths are given by , where x is a normally distributed random variate with standard deviation s 1  , limited to ±0.8, and w 0 is the chain-specific mean relative chamfer width.

Construction of single-stranded chains
The construction of a single-stranded chain with N magnetosomes begins with the definition of chain-specific parameters that control the characteristics of individual crystals and their arrangement in space (Table S2). Chain-specific parameters for individual crystals are: the width s 0 , the elongation e 0 , and the relative chamfer width w 0 (prismatic crystals only) of ideal mature magnetosomes, the concentration parameters s and c of the Fisher-Von Mises distributions used to model shape and crystalline axes orientations with respect to the chain axis, and the standard deviation s of relative variations of axis lengths and chamfer widths. Other parameters control the chain geometry: the chain length, through the number N of magnetosomes and the gaps between them, the off axis magnetosome offsets, the crystal size tapering towards the chain extremities, and the chain bending.
The size tapering of the n-th magnetosome t t t t min , , n N n t n s s n n  The tapering function is thus controlled by the relative size t s of the end magnetosomes and the number t n N 2 / of end magnetosomes affected by tapering.
The arrangement of magnetosomes in a straight chain is controlled by the mean gap g 0 between neighbor crystals, and by the longitudinal and transversal randomization of their positions. Longitudinal gaps are defined as the mean distance between facing cylinder surfaces used to model prismatic crystals (Fig. S9c), or the shortest distance between the ellipsoidal surfaces used to model equant crystals (Fig. S9d). The longitudinal randomization of magnetosome positions causes individual gaps to differ from the expected value g 0 , and is realized by defining these gaps as g g g = + 0 1 ( ) , where g is a normally distributed random variate with standard deviation z 1  , limited to ±1 . Furthermore, g must be sufficiently large to guarantee a predefined minimum distance min g 0 for the facing surfaces. In intact native chains, the minimum distance is imposed by the thickness of the magnetosome membrane, which has been reported to vary between ~3 and ~7 nm (Gorby et al., 1988;Werckmann et al., 2017;Mickoleit et al., 2018), yielding min g values comprised between ~6 and ~14 nm, depending on magnetosome size. It is not clear whether the original gap is preserved after cell dissolution, until what remains of the original chains is fixed in the sediment matrix, for instance by adhesion on clay platelets (Galindo-González et al., 2009).
Transversal randomization of the magnetosome positions with respect to the original axial alignment along the z-direction of a cartesian coordinate system is described by distance vectors , ) from the chain axis, where x and y are normally distributed random variates with standard deviation x .
The straight shape of ideal chains can be altered in several ways, which include bending, kinking, and various forms of collapse that might occur after cell dissolution. Chain bending is the elastic response to small external loads and is commonly observed in intact MTB cells (Shcherbakov et al., 1997;Orue et al., 2018). The maximum elastic bending angle, defined as the angle between segments connecting consecutive crystal centers, is of the order of g s 0 0 / 0.1 (Shcherbakov et al., 1997). Above this limit, kinking occurs, until the chain extremities become close enough to form collapsed structures with magnetic flux closure, such as folded chains (e.g., Fig. 2b in Lins and Farina, 2004), rings, and handles (Kiani et al., 2018). Chain bending is therefore the only deviation of non-collapsed single chains from the ideal straight shape. The typical length scale over which a single magnetosome chain remains straight under the influence of thermal activation is >30 times larger than magnetosome size (Kiani et al., 2015), which means that the maximum value of the bending angle b between consecutive magnetosomes is of the order of few degrees. Chain bending is realized by rotating each magnetosome by the bending angle b with respect to the previous one, about a rotation axis that intersects the original chain axis while being perpendicular to it and tangential to the crystal surface (Fig. S10b-c). As a result, the mean gap between crystals does not change, as long as the bending angle is sufficiently small to maintain the minimum distance min g between crystals. Otherwise, the longitudinal gap is increased until this condition is met.
In practice, chain construction is based on the following steps: (1) Build N magnetosomes as described in Section 1.2, with predefined crystal orientations with respect to a straight chain axis parallel to z.
(2) Define N 1 − longitudinal gaps using the parameters g 0 , min g , z , and taking the magnetosome orientations into account.
(3) Construct a straight chain by placing the center of the first magnetosome at the origin, and subsequently adding the remaining magnetosomes with their centers on the z-axis, and distances along z defined by the N 1 − gaps calculated in (2). At this stage, the chain axis coincides with the z-axis of the coordinate system and all magnetosome centers lie exactly on the chain axis. Gaps are random realizations of a normal distribution with mean g 0 and standard deviation z , with upper limit g 0 2 and lower limit set by the minimum distance min g .
(4) Transversal magnetosome positions are randomized by adding , , x y 0 ( ) to the center coordinates of each magnetosome, where x and y are normally distributed random variates with standard deviation x .
(6) Chain bending is applied. For this purpose, the straight chain is first rotated randomly about the z-axis. Then, all n 2 ≥ magnetosomes are rotated by b about the y-axis. Next, all n 3 ≥ magnetosomes are rotated using the same procedure, and so on, until the last magnetosome has been rotated. The bended chain is rotated back about the z-axis by the same azimuthal angle used to prepare chain bending, so that the bending axis orientation is random within the xy-plane.

Construction of double-stranded chains
Several types of MTB produce strands of two or more closely arranged magnetosome chains made of equidimensional, prismatic, or tooth-shaped magnetosomes (Spring et al., 1994;Pósfai et al., 2006;Isambert et al., 2007;Faivre and Schüler, 2008;Zhang et al., 2017). In these double-or multi-stranded chains, the arrangement of magnetosomes is staggered, such that crystals of one chain face the gap between consecutive crystals in the other chain. This arrangement is energetically most favorable (Hanzlik et al., 2002;Ruder et al., 2012) to the formation of chain bundles with a consistent magnetic polarity (Simpson et al., 2005). Double-stranded chains can also form after cell death if a single chain is bent beyond the elastic limit: in this case the magnetostatic attraction of the two extremities causes the chain to fold around a kink point (Pósfai et al., 2006). If this type of collapse occurs in a weak field, complete folding produces two parallel chain fragments with opposed magnetic moments that form a closed magnetic loop. Strands with opposed magnetic moments tend to produce a side-by-side arrangement of the crystals that maximizes lateral attractive forces (Varón et al., 2013), as seen on attracted segments of extracted chains (Zhu et al., 2016). Side-by-side arrangement and size tapering only at one extremity (e.g., Fig. 2b in Kiani et al., 2015) are distinctive features of fold-collapsed chains, which contrasts with those of native double-stranded chains.
The construction of a double-stranded chain of N N N = + 1 2 magnetosomes, with two parallel strands of N 1 and N 2 magnetosomes each, begins with the definition of the magnetosome properties of the two strands and the corresponding longitudinal gaps l g as in section 1.2. In addition, transversal gaps t g are used to set the distance between the two subchains. Transversal gaps are defined in the same manner as longitudinal gaps, using the parameters g 0 , min g , and z . The overall shape of double-stranded chains is governed by two degrees of freedom: bending, as with single chains, and twisting. Twisting consists in the rotation of the two sub-chain axes around a common axis of the double chain. Details of the double chain construction depend on its origin (intact or collapsed), as specified in the following.

Intact double-stranded chains
Ideal native double-stranded chains consist of two staggered and identically tapered subchains of similar length, whose magnetosomes face a common axis at a constant distance nearly equal to that of the longitudinal gaps (Fig. S11a). The longitudinal position of the second strand relative to the first one is defined by the integer lag number l , which describes the lag between the end magnetosomes at one extremity of the double chain in units of magnetosome length. Accordingly, the difference between the z-coordinates of these end magnetosomes is given by z l s g ∆ = + + 0 0 0 1 2 ( / )( ), where the half-integer factor l +1 2 / accounts for the staggered position of magnetosomes in strand 2 facing the gaps between magnetosomes in strand 1. Double-stranded chains can be twisted (Fig. S11b) and bended (Fig. S11c) by rotating and rearranging individual magnetosomes so to maintain, on average, similar distances between crystals. Deviations from the ideal double-stranded chain structure include random variations of magnetosome shape and position, length of the two strands, and longitudinal lag of one strand with respect to the other.
Twisting converts the two straight and parallel strand axes (Fig. S11a) into a double helix whose pitch t s g + 0 0 0 2 ( ) / is controlled by the pitch angle t , defined as the angle by which the distance vector between the strand axes rotates over the length occupied by a single magnetosome. The helix pitch is much longer than the size of individual magnetosomes, so that a half-turn is typically completed over ~10 magnetosomes or more (Fig. S11b). Twisting requires to rotate the magnetosome centers about the common chain axis, and the magnetosomes about the distance vectors connecting their centers with the chain axis. As with singlestranded chains, magnetosome rotations are performed so to maintain the predefined gaps, compatibly with the given minimum gap min g . Bending of double-stranded chains is more complicated than the single-stranded chain case, because the bending axis might be in-plane with the two sub-chains: in this case, the gap between consecutive magnetosomes lying further away from the bending axis needs to be increased to maintain the original staggered arrangement (Fig. S11c). Figure S11: Construction of ideal double-stranded chains by twisting and bending of an initial straight configuration. In these examples, prismatic magnetosomes with elongation e 0 are modeled by chamfered cylinders. (a) Straight double-stranded chain of staggered magnetosomes with size tapering at both ends. This example with N = 1 9 , N = 2 8 , and l = 0 is similar to the double chain shown in Fig. 3 of (Simpson et al., 2005). (b) Same as (a), after twisting the double chain plane by a half-turn ( t =°20 ). (c) Same as (b), after bending the twisted chain with a bending angle b =°6 , which is close to the elasticity limit calculated in (Shcherbakov et al., 1997). Note the larger gaps on the external bending side (right). (d-f) Same as (a-c) for a double-stranded chain obtained by folding a single-stranded, tapered chain of 16 magnetosomes. This example with N N = = 1 2 8 , N = 2 8 is similar to the double-stranded chain shown in Fig. 2 of Lins and Farina (2004).
The construction of staggered double-stranded chains is based on the following steps: (1) Build N 1 magnetosomes for strand 1 and N 2 magnetosomes for strand 2 as described in Section S1.2, with predefined crystal orientations with respect to a straight chain axis parallel to z and sizes determined by a common tapering function.
(2) Construct N 1 -and N 2 -vectors of the form , , x y 0 ( ) defining the transversal randomization of the magnetosome positions in the two strands, with x and y being normally distributed random variates with standard deviation x .
(3) Define N − 1 1 and N − 2 1 longitudinal gaps using the parameters g 0 , min g , z , and taking the magnetosome orientations into account. Define N 1 and N 2 transversal halfgaps with parameters g 0 , min g , z , which will be used to set the minimum distance between the two strands.
(4) Construct the first strand using the longitudinal gaps defined in (2). Unlike singlestranded chains, the right crystal edges are aligned with the chain axis, instead of the magnetosome centers through where k x is the x-coordinate of the center of the k-th magnetosome in strand 1, k x the corresponding transversal randomization, k s the crystal width, and t k g 2 / the k-th transversal half-gap. Magnetosome center coordinates along y are defined by k k y y = as for single-stranded chains.
(5) Construct the second strand using the longitudinal gaps defined in (2). The left crystal edges are aligned with the chain axis through where k x is the x-coordinate of the center of the k-th magnetosome in strand 2, k x the corresponding transversal randomization, k s the crystal width, and t k g 2 / the k-th transversal half-gap.
Magnetosome center coordinates along y are defined by k k y y = as for single-stranded chains. The z-coordinate of the first magnetosome of strand 2 is identified with that of the midpoint of the gap between the (l +1)-th and the (l +2)-th magnetosome of strand 1, where the integer l 0 is the lag between the two strands.
(6) Because of random variations in magnetosome dimension and gap sizes, magnetosomes are not perfectly staggered. The best possible staggering is obtained by moving strand 2 along z until the sum of the squared differences , between the longitudinal position of the magnetosomes in strand 2 and the corresponding gap midpoint in strand 1 is minimized. Finally, the overall separation between the two strands is adjusted by adding or subtracting the same offset along x, until transversal gaps between corresponding magnetosomes are as small as possible, but larger than the minimum gaps prescribed by the transversal half-gaps t k g 2 / defined in (5).
(7) Chain twisting is applied. The starting point is a double-stranded chain with a straight axis along z and magnetosome centers lying in the xz-plane, up to small random displacements along y. Twisting consists in rotating the center coordinates of all magnetosomes about the z-axis by an angle t , (8) Chain bending is applied. For this purpose, the twisted double-stranded chain is first rotated randomly about the z-axis. Then, all n 2 magnetosomes in strand 1 and all n l + 1 magnetosomes in strand 2 are rotated by b about the y-axis. The rotation axis is fixed with respect to the second magnetosome of strand 1 or the (l +1)-th magnetosome of strand 2, depending on which of the two is closer to the bending axis. Next, all n 3 ≥ magnetosomes in strand 1 and all n l 2 ≥ + magnetosomes in strand 2 are rotated using the same procedure, and so on, until the last magnetosome has been rotated. Finally, the bended chain is rotated back about the z-axis by the same angle used to prepare chain bending, so that the bending axis orientation is random within the xy-plane.

Fold-collapsed double-stranded chains
The ideal fold-collapsed chain consists of two parallel strands obtained by complete folding of a tapered single-stranded chain (Fig. S11d). Magnetosomes in the two strands are facing each other, whereby this condition is exactly met at the folded side. Tapering occurs only at the other side of the chain, which correspond to the two extremities of the original singlestranded chain. In the ideal case, folding occurs in the middle, so that the number of magnetosomes in the two strands differs at most by 1. The remaining properties of fold-collapsed chains, including twisting, bending, and random deviations from the ideal shape, are defined in the same manner as for native double-stranded chains (section 3.1).
The construction of fold-collapsed chains is based on the following steps: (1-4) Same as the corresponding steps for staggered double-stranded chains.
(5) Construct the second strand using the longitudinal gaps defined in (2). The left crystal edges are aligned with the chain axis through where k x is the x-coordinate of the center of the k-th magnetosome in strand 2, k x the corresponding transversal randomization, k s the crystal width, and t k g 2 / the k-th transversal half-gap.
Magnetosome center coordinates along y are defined by k k y y = as for single-stranded chains. The z-coordinate of the first magnetosome in strand 2 is the same as that of the first magnetosome in strand 1.
(6-7) Chain twisting and bending is applied in the same manner as for staggered doublestranded chains.

Energy barrier calculations
The role of chain geometry for the stabilization of high-moment states is investigated by considering idealized chain geometries made of identical, equidimensional magnetosomes with the spontaneous magnetization of magnetite ( s kA/m μ = 480 ) and no magnetocrystalline anisotropy. In this case, the total free magnetic energy F in a null field is entirely controlled by point dipole interactions between crystals (Hanzlik et al., 2002). For a chain of N crystals where k m is the magnetic moment vector of the k-th magnetosome, ki r the center-to-center distance between magnetosomes k and i, and ki r the unit vector parallel to the distance vector. In the context of energy barrier calculations, it is convenient to normalize eq. (5) by the energy B k T of thermal fluctuations at the absolute temperature T . In this case, and for spherical magnetosomes of diameter D separated by a gap gD , eq. (5) becomes In case of magnetite at room temperature ( The thermally activated transition between two stable magnetic states, defined by local minima of eq. (6), is represented by a continuous path connecting the two local minima in the N 2 -dimensional parameter space ( , , , , , ) formed by the polar angles k and azimuthal angles k of the N magnetic moment vectors. This path is a function of the fractional path length s , with ( ) s = 0 and ( ) s =1 yielding the coordinates of the two local minima. The energy required for a transition along is given by max ( ( )) ( ( )) s s − 0 . Because the transition probability is strongly dependent on the energy, transitions will occur along paths that minimize the energy increase. These paths are characterized by an energy barrier ∆ that minimizes the energy increase. Finding ∆ requires is equivalent to the search of a global minimum for max ( ( )) ( ( )) s s − 0 with respect to all possible paths with same starting and ending points in a N 2 parameter space. This problem is solved by minimizing the geometric action (64) is an infinitesimal change along , and ∇ is the gradient of the normalized free energy given in eq. (6). In practice, eq. (8) is solved numerically starting from an initial path = x 0 0 that is close enough to the final solution and updating this path iteratively with an updating scheme that decreases S . Here, we use the updating scheme proposed by Fabian and Shcherbakov (2018), where ( ) i s t is the tangential vector of i x at the position s along the path, and k an empirical step size that is updated at every iteration on the basis of the observed change of ( ) i S x . In practice, eq. (9) moves laterally (i.e., perpendicular to x  ) along the direction that decreases the total energy, until i t becomes parallel to the energy gradient. This ensures that the final path contains one or more saddle points that define the height of the energy barrier.
The choice of the initial path is crucial for the convergence of eq. (9) to a global minimum of the energy increase along . Short linear chains of 2-4 spheres reverse by a fanning mechanism (Jacobs and Bean, 1955) given by ( ) k k s = −1 and k = 0 , where even and odd magnetic moments rotate clockwise and counterclockwise, respectively, by 180°. In this case, the energy barrier is given by the difference between the interaction energy of dipoles arranged side-byside and head-to-tail, respectively. This reversal mode yields, to a first approximation, if only nearest-neighbor interactions are considered. The energy barrier associated with symmetric fanning grows linearly with chain length and is replaced by more favorable reversal modes for chains containing more than four particles. Newell (2009) found analytically that long single-stranded chains reverse through the nucleation of a reversed domain at one extremity. The reversed domain grows with little additional energy, pushing the boundary between normal and reversed magnetic moments towards the other chain extremity, until the original domain is denucleated. In this case, the energy barrier is mainly controlled by the domain nucleation process and is only marginally affected by the number of magnetosomes when N >5 . With these considerations in mind, the initial reversing path for a singlestranded chain is constructed from a two-domain configuration obtained by numerical minimization of eq. (6) with starting parameters k = 0 for k n < , k = for k n > , / , n = 2 and k = 0 , where n denotes the middle magnetosome (Fig. S7a). If necessary (e.g., for ( )/ n N +1 2), the two-domain configuration is stabilized by fixing / n = 2 . The energy associated with this configuration is larger than that of stable states where all magnetic moments are characterized by k = 0 or k = . Relaxation to one of the two stable states is thus obtained by relaxing the two-domain configuration after introducing a little modification of the dipole interaction strength that favor the propagation of the transition zone between the two domains to one of the two extremities. For instance, relaxation to k = 0 can be obtained by multiplying all terms with k n < in eq. (6) by a fixed factor <1. The initial transition path 0 is then obtained by joining the reversed path for the relaxation to k = 0 with the path for the relaxation to k = , realizing a continuous transition from an initial state with k = 0 to a final state with k = through the previously defined two-domain configuration. Because the relaxation of a two-domain configuration is readily triggered by values of very close to unit, 0 represents already a good approximation of the true transition path, and convergence of eq. (9) to a global minimum can be expected.
The energy barrier of magnetosome rings and double-stranded chains is calculated in a similar manner. The high-moment state of a large ring of spherical particles, obtained by decreasing a saturating magnetic field to zero, is a two-domain state that divides the ring into two halves with tangential magnetic moments arranged clockwise and counterclockwise, respectively, separated by two magnetic moments pointing to the same direction along the ring diameter (Fig. S7b). In case of large rings ( N 1  ) the domain boundaries can be moved along the circumference with little additional energy. Therefore, the transition to a zeromoment state with clockwise or counterclockwise tangential moments is obtained by reducing the extension of one of the two domains until it is fully denucleated. Double-stranded chains and their magnetic configurations are topologically equivalent to those of a ring of particles "squeezed" by transforming it into an ellipse with unit eccentricity, whose sides coincide with the two strands. Accordingly, the high-and low-moment states of a double-stranded chain are obtained when the two strands contain axial magnetic moment with same and opposite polarity, respectively (Fig. S7c). Magnetic moments at the extremities deviate from the chain axis, due to the stray field produced by the dipolar interactions between the two strands. Migration of one domain boundary in a ring structure away from the symmetric position of the high-moment state is equivalent to the nucleation of a domain boundary in one of the two strands of a double-stranded chain, and subsequent displacement along the chain axis. Therefore, the transition between low-and high-moment states of a double-stranded chains is calculated by relaxation of an intermediate configuration where one of the two stains is divided into two domains as it was the case for a single-stranded chain.
In summary, a similar procedure is used to calculate the zero-field transition between two stable states for all four chain structures considered in this article: single-stranded chains, rings, double-stranded chains with staggered magnetosomes, and fold-collapsed chains. This procedure consists of the following steps: (1) For all non-ring configurations: construct a transition configuration in which one strand of the chain is divided into two domains by setting / = 2 for the middle magnetosome and relaxing the remaining magnetic moments.
(2) Calculate the initial transition path by making one of the two domains grow at the cost of the other. This is accomplished by decreasing the strength of dipole interactions within the shrinking domain. Depending on the choice of the shrinking domain, denucleation of this domain yields the high-and the low-moment state.
(3) Extend the transition path calculated in (2) by relaxing the path ends to the final highand the low-moment configurations after resetting all dipole interaction strengths to their original strength according to eq. (6).

Micromagnetic FORC simulations
Micromagnetic calculations are needed to obtain the component of the total chain magnetic moment along the direction of the applied field while the applied field changes according to a predefined FORC measurement protocol. These calculations simulate VSM measurements, where the specimen magnetic moment parallel to the applied field direction is measured. Full micromagnetic simulations require each magnetosome to be divided into volume elements whose size does not exceed the exchange length et al., 2013), where ex A and s μ are the exchange constant and the spontaneous magnetization of the magnetic mineral ( ex nm l 10 for magnetite). The minimum size for the nucleation of vortex states is ex l 6 , which corresponds to ~60 nm in equant magnetite particles (Enkin and Williams, 1994;Newell and Merrill, 2000). The critical size of magnetosomes is 20-50% larger (Fig. S1), due to the stabilizing effect of magnetostatic interactions (Muxworthy and Williams, 2006). Accordingly, vortex states are not nucleated by FORC measurements if the magnetosome width does not exceed 70 nm. In this case, the magnetization of individual crystals is nearly homogeneous, with some deviations resembling a flower bundle appearing between 50 and 70 nm (Berndt et al., 2020).
Regular magnetofossil sizes fall within the stable single-domain limits (Fig. S1), so that their magnetization can be effectively considered homogeneous. In this case, the magnetization of a whole magnetosome is fully described by the bulk magnetic moment vector s V μ = m u, where V is the crystal volume and u a unit vector representing the direction of m. This simplification greatly accelerates micromagnetic simulations, enabling the calculation of highresolution FORC diagrams for thousands of magnetosome chains.

Energy minimization
The magnetic state of a single or double chain containing N magnetosomes is completely specified by the vectors , , , N = 1 2  ( ) and , , , N = 1 2  ( ) of polar and azimuthal angles of the magnetic moment orientations , , , N u u u 1 2  , with cos sin ,sin sin , Starting from a well-defined initial state, such as positive saturation at the beginning of each FORC, simulated measurements in the external field H dictated by the FORC protocol are obtained by successive minimizations of the total free magnetic energy of the chain. Here a distinction is made between the magnetic field H (SI unit: A/m) and the magnetic induction B (SI unit: T) improperly identified with the field applied during magnetic measurements through the fixed proportionality μ = B H 0 between the two quantities in air. The free magnetic energy is conveniently normalized by The field direction is specified by cos sin ,sin sin ,cos = h ( ) , where is the polar angle between the tangent of the chain axis at its midpoint and the field direction, and the azimuthal angle of the field direction with respect to the chain plane.
The normalized free energy of shape anisotropy, or demagnetizing energy, is given by where D is the magnetometric (i.e., volume-averaged) demagnetizing tensor (Fukuma and Dunlop, 2006) of a homogeneously magnetized body with given shape. The demagnetizing tensor is defined so, that D M is the volume-averaged internal field caused by the homogeneous magnetization . M Analytical solutions exist for ellipsoids and cylinders. In case of a triaxial ellipsoid with principal axes a a a 3 2 1 , the internal field is homogeneous, and the demagnetizing tensor has the principal components (Osborn, 1945) [ with D D D + + = In the particular case of prolate rotation ellipsoids, which are used to model equidimensional magnetosomes, the limit of eq. (15) for a a − 2 1 0 yields and The demagnetizing tensor of a circular cylinder has two identical principal components perpendicular to the cylinder axis, and a third principal component (Joseph, 1966;Chen and Brug, 1991) parallel to the cylinder axis, with , K F = 2 ( / ) and , E E = 2 ( / ) being the complete elliptic integrals of the first and second kind, respectively, and the length-to-diameter ratio ‡ . The demagnetizing tensor of chamfered cylinders cannot be resolved analytically; however, the effect of chamfering on D 3 does not exceed 10% for the usual geometries required to approximate prismatic magnetosomes. Therefore, excellent numerical approximations can be obtained using the following empirical correction factor: where , D w 3 ( ) is the longitudinal demagnetizing factor of a cylinder with length-to-diameter ratio and normalized chamfering width w , and , D 3 0 ( ) longitudinal demagnetizing factor of the same cylinder without chamfering, as given by eq. (19) (Fig. S12). The maximum error of this approximation is <1% for the typical ranges of λ and w that characterize prismatic magnetosomes.
The interaction energy between two magnetized particles 1 and 2 is equivalent to the sum of the potential energies of the interaction matrix , F 1 2 .
In case of equidimensional magnetosomes, each particle is represented by a rotation ellipsoid with diameter i s and elongation i , whose stray field coincides with that of a centered point dipole with same dipole moment vector. In this case, eq. (22) Using this description of the cylinder surfaces, the interaction matrix is given by where positive and negative signs apply to the top and bottom faces of , where , S 1 2 are the surfaces specified in eq. (26-27). The matrix coefficients in eq. (29) need to be calculated only once for a complete FORC simulation of a given chain configuration, since they depend only on geometric factors. Interaction coefficients are very sensitive to the geometry, relative position and orientation of nearest neighbor magnetosomes, while converging to the dipolar approximation of eq. (25) over longer ranges.

FORC simulations
FORC simulations of individual magnetosome chains are made of few curves (Fig. S13) from which the whole set of measurements corresponding to a complete FORC protocol is reconstructed. This simplification is possible because isolated chains possess only few magnetic states. Each simulated FORC starts in a positive saturating field sat H that is sufficiently large for all chains to possess only one stable magnetic state. Next, the field is decreased in small regular steps H corresponding to the resolution of the FORC simulation. At each step, a new state is calculated using the , , , N = 1 2  ( ) and , , , N = 1 2  ( ) values of the previous state as starting parameters for the energy minimization. This operation is repeated until the minimum reversal field min r H specified by the FORC protocol is reached. For a complete FORC simulation, the amplitude of min r H needs to exceed the field r,max H above which the major hysteresis loop closes (Fig. S13). The magnetic state calculated for each applied field step defines the component s sin cos ,sin sin ,cos of the magnetic moment parallel to the positive applied field direction, which simulates the contribution of a single chain to the bulk measurement.
Because calculations start at a positive saturating field, the pairs ,  (Lang and Schüler, 2006;Chen et al., 2007), while the application of large magnetic fields to aqueous suspensions, such as during magnetic extraction, promotes the formation of particle strings resembling magnetosome chains (Zhang et al., 2005;Barrett et al., 2011;Wang et al., 2013). On the other hand, the complete collapse of magnetosome chains in natural environments might be prevented by spontaneous adhesion of magnetite to larger sediment particles (Tombácz et al., 2001;Galindo-González et al., 2009).
Simulations in this article are limited to six types of hypothetical magnetofossil configurations that might be encountered in natural sediment. The first four types include intact single and double chains of equant and prismatic magnetosomes with size and shape distributions taken from 3441 magnetite crystals identified as magnetofossils (Fig. S1). The last two types include double chains of equant and prismatic magnetosomes resulting from the collapse of single chains when kinked beyond their elastic limit. This form of collapse is expected to occur in nature, since it is triggered by minor lateral forces that can be expected to occur in bioturbated sediment once the cell material has been dissolved. Other forms of collapse induced by specific treatments used for TEM observation, such as clumping (Kobayashi et al., 2006) and shrinkage ( Fig. 9 in Shcherbakov et al., 1997) are not considered, since they result from the application of strong mechanical and magnetic forces to aqueous suspensions, which are not present in nature.
Chain folding, however, is not necessarily the only type of chain alteration that might occur during early sediment diagenesis. The complete collapse of isolated chains might produce isolated small magnetosome clusters similar to those produced by the ΔmamJ mutant of M. gryphiswaldense; however the small remanence ratio of these clusters ( rs s M M / 0.23, Katzmann et al., 2013), which is typical for dense clusters of single-domain magnetite (Muxworthy et al., 2003), is incompatible with typical magnetofossil values (Ludwig et al., 2013). Therefore, significant contributions from this form of total collapse can be excluded. Other forms of alterations of the original chain structure include fragmentation, that is, the subdivision into shorter chain fragments. Certain magnetotactic bacteria strains, such as MV-1 (Kalirai et al., 2013) produce chains containing large gaps, at least under fast growing conditions typical of cultures. Subsection of those chains might be easily ripped apart after cell dissolution, through adhesion to different sediment particles. Fragmentation is thus expected to decrease the mean length of magnetofossil chains. The abundance of short chain fragments in magnetofossil extracts might reflect this fragmentation process, but it might also be an artifact of the sample preparation procedure for TEM observations. For modeling purposes, fragmented chains differ from intact chains only by the number of magnetosomes. The effect of chain length on bulk magnetic properties is large for chains fragments containing 2 or 3 magnetosomes and becomes negligible with 6 or more crystals (Chang et al., 2019). A minimum of 4 magnetosomes has been assumed for single chains or chain fragments, based on the typical minimum length of natural sections of chains containing large gaps (Meldrum et al., 1993b;Bazylinski et al., 1995;Kalirai et al, 2013;Blondeau et al., 2018;Monteil et al., 2018;Zhu et al., 2016).
The maximum length of simulated single chains is limited to 9 magnetosomes since longer chains are magnetically undistinguishable. A uniform chain length distribution between these two limits is assumed (Table S1).
The other parameters are chosen according to the criteria discussed in Sections 1-3, so that the generated chains resemble TEM observation reported in the literature, with respect to the average shape and random deviations from this shape (Table S1). Random examples of double-stranded chains generated with these parameters are shown in Fig. S2-5. For each of the abovementioned six chain configuration, ~10 5 examples like those in Fig. S2-5 have been generated, and for each of these examples, a single set of FORC measurements has been simulated assuming a pure magnetite composition, random orientation with respect to the applied fields, a maximum applied field of 0.35 T, and 1 mT field steps.
The full FORC simulation of a single-stranded chain takes ~4 min on average in case of equant magnetosomes, and ~11 min in case of prismatic magnetosomes, on a single processor core running at 2.7 GHz. The longer time required for prismatic magnetosomes is due to the more complex handling of magnetostatic interactions when the point dipole approximation cannot be used. Double-stranded chains take a significant longer time: 4.5 and 7.6 hours on average for equant and prismatic magnetosomes respectively, under the same conditions, due to the larger number of particles and nearest neighbor magnetostatic interactions to account for. The total CPU time of ~3×10 6 hours·core·GHz required to complete all simulations has been shared among 6 desktop and 20 Raspberry Pi 4 computers.
Results obtained for individual chains have been merged into a single set of simulated FORC measurements with unit saturation magnetization ( s M =1) for each of the six configurations (Fig. S6). Simulated measurements are apparently continuous; however, they contain small magnetization jumps from a finite number of individual chain contributions of the type shown in Fig. S14. Statistical fluctuations in the contribution of these jumps to the bulk magnetization become visible in the finite difference estimate ( ) M B = 2 / of the FORC function (Fig. S14). Therefore, FORC diagrams have been calculated with the same variable smoothing procedure VARIFORC (Egli, 2013) used to process noisy measurements. The smoothing factor . S B B = + 5 0 07 / , where B is the FORC diagram coordinate and B =1 mT measurement resolution, with a limitation to S = 2 along the central ridge, and to S = 9 at places with maximum curve slope suppress the statistical noise of simulated measurements while preserving high-resolution features such as the central ridge.
The full FORC simulation of a single chain takes ~4 min on average in case of equant magnetosomes, and ~11 min in case of prismatic magnetosomes, on a single processor core running at 2.7 GHz. The longer time required for prismatic magnetosomes is due to the more complex handling of magnetostatic interactions when the point dipole approximation cannot be used. Applied magnetic field vector, and parallel unit vector.
Coercive field or horizontal FORC coordinate, and its maximum value Reversal field, and corresponding minimum and maximum values.
Vertical FORC coordinate, and its minimum and maximum values Field step used in the FORC simulation First and second cubic magnetocrystalline anisotropy constant.
Lag of sub-chain 2 with respect to sub-chain 1 Number of tapered magnetosomes at each chain end.
Number of magnetosomes in a single chain, or in the two chains composing a double chain.
Full magnetosome width (size of mature magnetosomes), and width of end magnetosomes in tapered chains.
Chain tapering function, defining the size reduction of magnetosomes as a function of their position n inside the chain.

Mean chamfer width.
Chain bending angle.
Polar angle of magnetic moment, and vector of polar angles.
Angle between midpoint chain axis and applied field.
Concentration parameter for the Fisher-Von Mises distribution governing crystal axes misalignment with respect to shape axes.
Concentration parameter for the Fisher-Von Mises distribution governing shape axes misalignment with respect to the chain axis.
Standard deviation of the normal distribution governing the transversal randomization of magnetosome position.
Standard deviation of the normal distribution governing the longitudinal (gap) randomization of magnetosome position.
Standard deviation of the normal distribution governing the shape axes length randomization.
Azimuthal angle of magnetic moment, and vector of azimuthal angles Azimuthal angle between field direction and chain plane.

Parameter Explanation
, , a a a Magnetosome width.
Chain tapering function, defining the size reduction of magnetosomes as a function of their position n inside the chain.
Unit vector specifying the magnetic moment direction.
Individual magnetosome volume, and mean volume of mature magnetosomes

Exchange length
Magnetosome elongation