diff --git a/CHANGELOG.md b/CHANGELOG.md
index e92ec4573..de0d42023 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,15 @@
 
 In development.
 
+**New features**:
+
+- A new mutation model, ``SLiMv6MutationModel``, is introduced to match changes in
+  SLiM 6.0: this mutation behaves as ``SLiMMutationModel``, but without the ``type``
+  and ``slim_generation`` arguments, and stores different information in metadata
+  (the list of mutation IDs, in binary). ``SLiMMutationModel`` will remain for
+  backwards compatibility, but code for use with SLiM versions 6.0 and above should
+  use ``SLiMv6MutationModel`` ({issue}`2534`, {pr}`2539`, {user}`petrelharp`)
+
 ## [1.4.2] - 2026-05-20
 
 Bugfix release.
diff --git a/docs/api.md b/docs/api.md
index 94b772ba0..e1ecd58ef 100644
--- a/docs/api.md
+++ b/docs/api.md
@@ -48,6 +48,7 @@ for discussion and examples of individual features.
   TPM
   EL2
   InfiniteAlleles
+  SLiMv6MutationModel
   SLiMMutationModel
 ```
 
@@ -217,6 +218,10 @@ for discussion and examples of individual features.
 .. autoclass:: msprime.InfiniteAlleles()
 ```
 
+```{eval-rst}
+.. autoclass:: msprime.SLiMv6MutationModel()
+```
+
 ```{eval-rst}
 .. autoclass:: msprime.SLiMMutationModel()
 ```
diff --git a/docs/mutations.md b/docs/mutations.md
index 33d4d6638..6775948e2 100644
--- a/docs/mutations.md
+++ b/docs/mutations.md
@@ -105,8 +105,11 @@ for more information.
 {class}`.InfiniteAlleles` (Integers)
 : A generic infinite-alleles mutation model
 
+{class}`.SLiMv6MutationModel` (Integers)
+: An infinite-alleles model producing SLiM-style mutations for SLiM v6+
+
 {class}`.SLiMMutationModel` (Integers)
-: An infinite-alleles model producing SLiM-style mutations
+: An infinite-alleles model producing SLiM-style mutations for older versions of SLiM
 
 {class}`.MatrixMutationModel` (General finite state model)
 : Superclass of mutation models with a finite set of states
@@ -485,7 +488,8 @@ Here are the available models; they are documented in more detail below.
 - {class}`.TPM`: Two-phase mutation model for microsatellite repeat copy number. DiRienzo et al. ('94)
 - {class}`.EL2`: Two-phase mutation model, equal rate, linear bias model for microsatellite repeat copy number. Garza et al. ('95)
 - {class}`.InfiniteAlleles`: A generic infinite-alleles mutation model
-- {class}`.SLiMMutationModel`: An infinite-alleles model of mutation producing SLiM-style mutations
+- {class}`.SLiMv6MutationModel`: An infinite-alleles model of mutation producing SLiM-style mutations for SLiM v6+
+- {class}`.SLiMMutationModel`: An infinite-alleles model of mutation producing SLiM-style mutations for older versions of SLiM
 
 (sec_mutations_matrix_mutations_models)=
 
@@ -855,30 +859,34 @@ to set the starting allele appropriately, and to make sure the results make sens
 
 A special class of infinite alleles model is provided for use with [SLiM](<https://messerlab.org/slim/>),
 to agree with the underlying mutation model in SLiM.
-As with the InfiniteAlleles model, it assigns each new mutation a unique integer,
-by keeping track of the `next_id` and incrementing it each time a new mutation appears.
+As with the {class}`.InfiniteAlleles` model, it assigns each new mutation
+a unique integer, by keeping track of the `next_id` and incrementing it each
+time a new mutation appears. For more information,
+see {ref}`the pyslim documentation<pyslim:sec_tutorial_adding_neutral_mutations>`.
 
-This differs from the {class}`.InfiniteAlleles` because mutations
+This differs from the {class}`.InfiniteAlleles` model because mutations
 in SLiM can "stack": new mutations can add to the existing state, rather than
 replacing the previous state. So, derived states are comma-separated lists of
 mutation IDs, and the ancestral state is always the empty string. For instance,
 if a new mutation with ID 5 occurs at a site, and then later another mutation
 appears with ID 64, the sequence of alleles moving along this line of descent
-would be `""`, then `"5"`, and finally `"5,64"`. Furthermore, the mutation
-model adds SLiM metadata to each mutation, which records, among other things,
-the SLiM mutation type of each mutation, and the selection coefficient (which
-is always 0.0, since adding mutations in this way only makes sense if they are
-neutral). For this reason, the model has one required parameter: the `type`
-of the mutation, a nonnegative integer. If, for instance, you specify
-`type=1`, then the mutations in SLiM will be of type `m1`. For more
-information, and for how to modify the metadata (e.g., changing the selection
-coefficients), see
-{ref}`the pyslim documentation<pyslim:sec_tutorial_adding_neutral_mutations>`.
-For instance,
+would be `""`, then `"5"`, and finally `"5,64"`. The same list of integers is
+recorded in metadata, as binary integers rather than ASCII text.
+This redundancy between derived state and metadata is useful
+because the model behaves as an infinite-alleles model out-of-the-box; but if
+ancestral and derived states are modified (for instance, to be VCF-compliant)
+then there is no loss of information.
+
+There are in fact two mutation models for SLiM:
+{class}`.SLiMv6MutationModel` and {class}`.SLiMMutationModel`,
+because the format that SLiM uses to store mutations changed with version 6.0.
+Both store the same thing in derived state; the difference is in metadata.
+So, you probably want to use {class}`.SLiMv6MutationModel`;
+see {ref}`sec_mutations_mutation_slim_mutations_versions` for more on the differences.
 
 ```{code-cell} python
 
-model = msprime.SLiMMutationModel(type=1)
+model = msprime.SLiMv6MutationModel()
 mts = msprime.sim_mutations(
     ts, rate=1, random_seed=1, model=model)
 t = mts.first()
@@ -888,14 +896,16 @@ SVG(t.draw_svg(mutation_labels=ml, node_labels={}, size=(400, 300)))
 ```
 
 These resulting alleles show how derived states are built.
+(We're looking at derived states for convenience, but could equivalently
+look at metadata; see {ref}`sec_mutations_mutation_slim_mutations_metadata`.)
 
 The behaviour of this mutation model when used to add mutations to a previously mutated
 tree sequence can be subtle. Let's look at a simple example.
-Here, we first lay down mutations of type 1, starting from ID 0:
+Here, we first lay down mutations starting from ID 0 (the default):
 
 ```{code-cell} python
 
-model_1 = msprime.SLiMMutationModel(type=1)
+model_1 = msprime.SLiMv6MutationModel()
 mts_1 = msprime.sim_mutations(ts, rate=0.5, random_seed=2, model=model_1)
 t = mts_1.first()
 ml = {m.id: m.derived_state for m in mts_1.mutations()}
@@ -903,15 +913,16 @@ SVG(t.draw_svg(mutation_labels=ml, node_labels={}, size=(400, 300)))
 
 ```
 
-Next, we lay down mutations of type 2.
-These we assign starting from ID 100,
-to make it easy to see which are which:
-in general just need to make sure that we start at an ID greater than any
+Now suppose we wanted to add more mutations.
+(This is essentially what happens when we use msprime to add mutations
+to a tree sequence produced by SLiM.)
+We assign the new mutations starting from ID 100 to make it easy to see which are which:
+in general we just need to make sure that we start at an ID greater than any
 previously assigned.
 
 ```{code-cell} python
 
-model_2 = msprime.SLiMMutationModel(type=2, next_id=100)
+model_2 = msprime.SLiMv6MutationModel(next_id=100)
 mts = msprime.sim_mutations(
     mts_1, rate=0.5, random_seed=3, model=model_2, keep=True)
 t = mts.first()
@@ -928,7 +939,110 @@ with IDs `100` and `102`, between these two mutations.
 These were added to mutation `0`, obtaining alleles `0,100` and `0,100,102`.
 But then, moving down the branch, we come upon the mutation with ID `3`.
 This was already present in the tree sequence, so its derived state is not modified:
-`0,3`. We can rationalise this, post-hoc, by saying that the type 1 mutation `3`
-has "erased" the type 2 mutations `100` and `102`.
-If you want a different arrangement,
-you can go back and edit the derived states (and metadata) as you like.
+`0,3`. We can rationalise this, post-hoc, by saying that the mutation with ID `3`
+from the earlier batch of mutations has "erased" the mutations with IDs `100` and `102`
+from the later batch. If you want a different arrangement,
+you will need to edit the derived states and metadata directly.
+
+(sec_mutations_mutation_slim_mutations_metadata)=
+
+#### Metadata
+
+The most common use of this model is to add mutations to a tree sequence produced
+by SLiM, or by {func}`pyslim.annotate`; if so, the metadata will already be decoded.
+(In general, we recommend using {func}`pyslim.annotate` after
+{func}`.sim_ancestry` but before {func}`.sim_mutations` if you're using msprime
+to create a tree sequence for SLiM.)
+
+For completeness, here is a self-contained example
+that shows how to decode the metadata directly
+(but we emphasize that you probably *don't* need this in your code;
+it's provided for pedagogical reasons only!).
+If we look at the metadata in the example we've produced so far, we simply see binary
+(the metadata is not decoded because a metadata schema is not present):
+
+```{code-cell} python
+mut = mts.mutation(10)
+mut.metadata
+```
+
+To see the actual values, we need a {class}`tskit.MetadataSchema`,
+which we can add as follows:
+
+```{code-cell} python
+metadata_schema = tskit.MetadataSchema(
+    {
+        "codec": "struct",
+        "type": "object",
+        "properties": {
+            "derived_states": {
+                "items": {
+                    "binaryFormat": "q",
+                    "type": "number",
+                },
+                "noLengthEncodingExhaustBuffer": True,
+                "type": "array",
+            }
+        },
+        "required": ["derived_states"],
+        "additionalProperties": False,
+    }
+)
+
+tables = mts.dump_tables()
+tables.mutations.metadata_schema = metadata_schema
+mts = tables.tree_sequence()
+```
+
+(This is a simplified version of the schema used by SLiM 6,
+so it matches the metadata produced by the {class}`.SLiMv6MutationModel`.)
+Now, the metadata will be properly decoded, into lists of integer IDs:
+
+```{code-cell} python
+mut = mts.mutation(10)
+mut.metadata
+```
+
+(sec_mutations_mutation_slim_mutations_stacking)=
+
+#### More on "mutation stacking"
+
+Both SLiM models simply create the derived state by appending to the inherited state,
+and so if you use {func}`.sim_mutations` with ``keep=True`` to add SLiM mutations
+to a tree sequence that already has non-SLiM mutations, you may end up with
+mutations with odd derived states. For instance, if a SLiM mutations occurs
+at the same site as, and inheriting from, a previous mutation
+with derived state ``"A"``, the SLiM mutation might have derived state ``"A,345"``.
+For a mutation without a parent, this inherited state is the ancestral state
+of the site. If the site was not already present in the tree sequence, this model
+will create a new site with the empty string as the ancestral state,
+but if the site is already present and has a nonempty ancestral state,
+then this ancestral state will appear at the beginning of any derived states
+produced by these models at that site.
+
+The same thing is true of metadata: metadata is simply appended (in binary)
+to the "inherited" metadata, which is from the parent mutation (if any) or from
+the site (otherwise). So, if in your tree sequence your sites have nonempty metadata,
+then any SLiM mutations that occur on those sites will have site metadata
+prepended to their metadata (and so the metadata schema above will not work).
+However, these SLiM mutation models are designed specifically to work with SLiM-format
+tree sequences, which don't have any site metadata;
+if you're trying to use them in another situation then you probably
+don't have the right tool for the job.
+
+(sec_mutations_mutation_slim_mutations_versions)=
+
+#### Previous versions
+
+The technical difference between 
+{class}`.SLiMv6MutationModel` and {class}`.SLiMMutationModel` is as follows.
+{class}`.SLiMv6MutationModel` stores the SLiM IDs in binary metadata
+as well as in the derived state,
+so that the actual metadata about mutations can be looked up from where
+SLiM expects it in top-level metadata.
+The older version, {class}`.SLiMMutationModel`,
+stores the metadata about SLiM mutations directly in each tskit mutation's metadata
+(and since individual SLiM mutations can appear in more than one tskit mutation,
+this leads to some redundancy; thus the change
+to how mutation metadata is handled in SLiMv6). For more information, see
+{ref}`the pyslim documentation<pyslim:sec_previous_versions>`.
diff --git a/lib/msprime.h b/lib/msprime.h
index 742d17564..3a47b34e6 100644
--- a/lib/msprime.h
+++ b/lib/msprime.h
@@ -439,6 +439,8 @@ typedef struct {
     double *transition_matrix;
 } mutation_matrix_t;
 
+// this one is used by both SLiM mutation models
+// (with mutation_type_id and slim_generation unused by v6)
 typedef struct {
     int32_t mutation_type_id; // following SLiM's MutationMetadataRec
     int64_t next_mutation_id; // following SLiM's slim_mutationid_t
@@ -584,6 +586,8 @@ int matrix_mutation_model_alloc(mutation_model_t *self, size_t num_alleles,
     double *transition_matrix);
 int slim_mutation_model_alloc(mutation_model_t *self, int32_t mutation_type_id,
     int64_t next_mutation_id, int32_t slim_generation, size_t block_size);
+int slim_v6_mutation_model_alloc(
+    mutation_model_t *self, int64_t next_mutation_id, size_t block_size);
 int infinite_alleles_mutation_model_alloc(
     mutation_model_t *self, uint64_t start_allele, tsk_flags_t options);
 int mutation_model_free(mutation_model_t *self);
diff --git a/lib/mutgen.c b/lib/mutgen.c
index 946e1dec8..820168972 100644
--- a/lib/mutgen.c
+++ b/lib/mutgen.c
@@ -255,7 +255,94 @@ mutation_matrix_free(mutation_model_t *self)
 }
 
 /***********************
- * SLiM mutation model */
+ * SLiM mutation models */
+
+/* Code shared between the SLiMMutationModel and SLiMv6MutationModel */
+
+static int
+slim_both_mutator_choose_root_state(
+    mutation_model_t *MSP_UNUSED(self), gsl_rng *MSP_UNUSED(rng), site_t *site)
+{
+    site->ancestral_state = NULL;
+    site->ancestral_state_length = 0;
+    return 0;
+}
+
+static int
+slim_both_mutator_transition(mutation_model_t *self, const char *parent_allele,
+    tsk_size_t parent_allele_length, const char *parent_metadata,
+    tsk_size_t parent_metadata_length, mutation_t *mutation,
+    const tsk_size_t new_metadata_length, char **metadata_buffer)
+{
+    /* This function performs shared work for both SLiM models:
+     * - copying the inherited state into the derived state,
+     *   and appending the SLiM ID to that derived state
+     * - copying over the previous metadata
+     *   but NOT appending the new metadata (since that differs between models);
+     *   the *metadata_buffer argument here exists so that the model-specific
+     *   function knows where to put the new bit of metadata
+     * - it is that client code's responsibility to make sure that no
+     *   more than parent_metadata_length + new_metadata_length
+     *   gets copied into *metadata_buffer
+     * - it is also the client code's responsibility to increment
+     *   params->next_mutation_id, since that is needed for the new metadata
+     */
+    int ret = 0;
+    slim_mutator_t *params = &self->params.slim_mutator;
+    char *buff = NULL;
+    int len;
+    /* The maximum number of digits for a signed 64 bit integer (including
+     * the leading "-") */
+    const size_t max_digits = 20;
+    /* We allow for a possible comma to separate the previous element
+     * in the list, as well as a NULL byte added by snprintf. We don't bother
+     * trying to alloc the exact number of bytes needed. */
+    const size_t alloc_size = parent_allele_length + max_digits + 2;
+    const char *sep = parent_allele_length == 0 ? "" : ",";
+    const tsk_size_t metadata_length = parent_metadata_length + new_metadata_length;
+
+    /* Append to derived_state */
+    buff = tsk_blkalloc_get(&params->allocator, alloc_size);
+    if (buff == NULL) {
+        ret = MSP_ERR_NO_MEMORY;
+        goto out;
+    }
+    len = snprintf(buff, alloc_size, "%.*s%s%" PRId64, (int) parent_allele_length,
+        parent_allele, sep, params->next_mutation_id);
+    if (len < 0) {
+        /* Technically this can happen. Returning TSK_ERR_IO should result
+         * in Python code checking errno. */
+        ret = TSK_ERR_IO;
+        goto out;
+    }
+    tsk_bug_assert(len < (int) alloc_size);
+    mutation->derived_state = buff;
+    mutation->derived_state_length = (tsk_size_t) len;
+
+    /* Append to metadata */
+    buff = tsk_blkalloc_get(&params->allocator, metadata_length);
+    if (buff == NULL) {
+        ret = MSP_ERR_NO_MEMORY;
+        goto out;
+    }
+    memcpy(buff, parent_metadata, parent_metadata_length);
+    // tell client code where to put the new metadata
+    *metadata_buffer = buff;
+    mutation->metadata_length = (tsk_size_t) (metadata_length);
+out:
+    return ret;
+}
+
+static int
+slim_both_mutator_free(mutation_model_t *self)
+{
+    slim_mutator_t params = self->params.slim_mutator;
+    tsk_blkalloc_free(&params.allocator);
+    return 0;
+}
+
+/*******************************
+ * (Older) SLiM mutation model */
 
 /* Typedefs from MutationMetadataRec in slim_sim.h:
  * line 125 in v3.4, git hash b2c2b634199f35e53c4e7e513bd26c91c6d99fd9
@@ -337,79 +424,97 @@ slim_mutator_check_validity(slim_mutator_t *self)
     return ret;
 }
 
-static int
-slim_mutator_choose_root_state(
-    mutation_model_t *MSP_UNUSED(self), gsl_rng *MSP_UNUSED(rng), site_t *site)
-{
-    site->ancestral_state = NULL;
-    site->ancestral_state_length = 0;
-    return 0;
-}
-
 static int
 slim_mutator_transition(mutation_model_t *self, gsl_rng *MSP_UNUSED(rng),
     const char *parent_allele, tsk_size_t parent_allele_length,
     const char *parent_metadata, tsk_size_t parent_metadata_length, mutation_t *mutation)
 {
     int ret = 0;
+    char *buff;
     slim_mutator_t *params = &self->params.slim_mutator;
-    char *buff = NULL;
-    int len;
-    /* The maximum number of digits for a signed 64 bit integer (including
-     * the leading "-") */
-    const size_t max_digits = 20;
-    /* We allow for a possible comma to separate the previous element
-     * in the list, as well as a NULL byte added by snprintf. We don't bother
-     * trying to alloc the exact number of bytes needed. */
-    const size_t alloc_size = parent_allele_length + max_digits + 2;
-    const char *sep = parent_allele_length == 0 ? "" : ",";
+    const tsk_size_t new_metadata_length = SLIM_MUTATION_METADATA_SIZE;
 
-    /* Append to derived_state */
-    buff = tsk_blkalloc_get(&params->allocator, alloc_size);
-    if (buff == NULL) {
-        ret = MSP_ERR_NO_MEMORY;
-        goto out;
-    }
-    len = snprintf(buff, alloc_size, "%.*s%s%" PRId64, (int) parent_allele_length,
-        parent_allele, sep, params->next_mutation_id);
-    if (len < 0) {
-        /* Technically this can happen. Returning TSK_ERR_IO should result
-         * in Python code checking errno. */
-        ret = TSK_ERR_IO;
+    ret = slim_both_mutator_transition(self, parent_allele, parent_allele_length,
+        parent_metadata, parent_metadata_length, mutation, new_metadata_length, &buff);
+    if (ret < 0) {
         goto out;
     }
-    tsk_bug_assert(len < (int) alloc_size);
+
+    // copy the final bit of metadata in (this is different than the v6 model)
+    copy_slim_mutation_metadata(params, buff + parent_metadata_length, mutation->time);
+    mutation->metadata = buff;
+
+    // finally, set up for the next one
     if (params->next_mutation_id == INT64_MAX) {
         ret = MSP_ERR_MUTATION_ID_OVERFLOW;
         goto out;
     }
     params->next_mutation_id++;
-    mutation->derived_state = buff;
-    mutation->derived_state_length = (tsk_size_t) len;
 
-    /* Append to metadata */
-    buff = tsk_blkalloc_get(
-        &params->allocator, parent_metadata_length + SLIM_MUTATION_METADATA_SIZE);
-    if (buff == NULL) {
-        ret = MSP_ERR_NO_MEMORY;
+out:
+    return ret;
+}
+
+/*************************
+ * SLiMv6 mutation model */
+
+/* To be fully SLiM-compatible, there also needs to be some information in
+ * top-level metadata.  But, it is not our job to add it in here, and that
+ * should be added after the fact with the method pyslim.add_mutation_metadata().
+ * */
+
+static void
+slim_v6_mutator_print_state(mutation_model_t *self, FILE *out)
+{
+    slim_mutator_t params = self->params.slim_mutator;
+    fprintf(out, "SLiMv6 mutation model :: next mutation ID = %d\n",
+        (int) params.next_mutation_id);
+}
+
+static int MSP_WARN_UNUSED
+slim_v6_mutator_check_validity(slim_mutator_t *self)
+{
+    int ret = 0;
+
+    if (self->next_mutation_id < 0) {
+        ret = MSP_ERR_BAD_SLIM_PARAMETERS;
         goto out;
     }
-    memcpy(buff, parent_metadata, parent_metadata_length);
-    copy_slim_mutation_metadata(params, buff + parent_metadata_length, mutation->time);
 
-    mutation->metadata = buff;
-    mutation->metadata_length
-        = (tsk_size_t) (parent_metadata_length + SLIM_MUTATION_METADATA_SIZE);
 out:
     return ret;
 }
 
 static int
-slim_mutator_free(mutation_model_t *self)
+slim_v6_mutator_transition(mutation_model_t *self, gsl_rng *MSP_UNUSED(rng),
+    const char *parent_allele, tsk_size_t parent_allele_length,
+    const char *parent_metadata, tsk_size_t parent_metadata_length, mutation_t *mutation)
 {
-    slim_mutator_t params = self->params.slim_mutator;
-    tsk_blkalloc_free(&params.allocator);
-    return 0;
+    int ret = 0;
+    char *buff;
+    slim_mutator_t *params = &self->params.slim_mutator;
+    const tsk_size_t new_metadata_length = sizeof(int64_t);
+
+    ret = slim_both_mutator_transition(self, parent_allele, parent_allele_length,
+        parent_metadata, parent_metadata_length, mutation, new_metadata_length, &buff);
+    if (ret < 0) {
+        goto out;
+    }
+
+    // copy the final bit of metadata in (this is different than the non-v6 model)
+    memcpy(
+        buff + parent_metadata_length, &params->next_mutation_id, new_metadata_length);
+    mutation->metadata = buff;
+
+    // finally, set up for the next one
+    if (params->next_mutation_id == INT64_MAX) {
+        ret = MSP_ERR_MUTATION_ID_OVERFLOW;
+        goto out;
+    }
+    params->next_mutation_id++;
+
+out:
+    return ret;
 }
 
 /***********************
@@ -558,10 +663,10 @@ slim_mutation_model_alloc(mutation_model_t *self, int32_t mutation_type_id,
 
     memset(self, 0, sizeof(*self));
 
-    self->choose_root_state = &slim_mutator_choose_root_state;
+    self->choose_root_state = &slim_both_mutator_choose_root_state;
     self->transition = &slim_mutator_transition;
     self->print_state = &slim_mutator_print_state;
-    self->free = &slim_mutator_free;
+    self->free = &slim_both_mutator_free;
     if (block_size == 0) {
         /* 8K is a good default, but we need to have the
          * block_size argument here because the size of the allocations
@@ -588,6 +693,43 @@ slim_mutation_model_alloc(mutation_model_t *self, int32_t mutation_type_id,
     return ret;
 }
 
+int MSP_WARN_UNUSED
+slim_v6_mutation_model_alloc(
+    mutation_model_t *self, int64_t next_mutation_id, size_t block_size)
+{
+    int ret = 0;
+    slim_mutator_t *params = &self->params.slim_mutator;
+
+    memset(self, 0, sizeof(*self));
+
+    self->choose_root_state = &slim_both_mutator_choose_root_state;
+    self->transition = &slim_v6_mutator_transition;
+    self->print_state = &slim_v6_mutator_print_state;
+    self->free = &slim_both_mutator_free;
+    if (block_size == 0) {
+        /* 8K is a good default, but we need to have the
+         * block_size argument here because the size of the allocations
+         * we need to make are in principle unbounded since SLiM
+         * copies the entire parent state for every mutation as it
+         * goes down along the tree.
+         */
+        block_size = 8192;
+    }
+    ret = tsk_blkalloc_init(&params->allocator, block_size);
+    if (ret != 0) {
+        ret = msp_set_tsk_error(ret);
+        goto out;
+    }
+    params->next_mutation_id = next_mutation_id;
+
+    ret = slim_v6_mutator_check_validity(params);
+    if (ret != 0) {
+        goto out;
+    }
+out:
+    return ret;
+}
+
 int MSP_WARN_UNUSED
 infinite_alleles_mutation_model_alloc(
     mutation_model_t *self, uint64_t start_allele, tsk_flags_t MSP_UNUSED(options))
diff --git a/lib/tests/test_mutations.c b/lib/tests/test_mutations.c
index ccc8d78f4..48c5e5abf 100644
--- a/lib/tests/test_mutations.c
+++ b/lib/tests/test_mutations.c
@@ -37,6 +37,9 @@ insert_single_tree(tsk_table_collection_t *tables, int alphabet)
         "0  1   6   4,5\n";
     */
     int ret;
+    const char *mut_metadata = "lmnop";
+    const char *site_metadata = "abc";
+
     tables->sequence_length = 1.0;
     ret = tsk_node_table_add_row(
         &tables->nodes, TSK_NODE_IS_SAMPLE, 0.0, 0, TSK_NULL, NULL, 0);
@@ -73,18 +76,18 @@ insert_single_tree(tsk_table_collection_t *tables, int alphabet)
     ret = tsk_population_table_add_row(&tables->populations, NULL, 0);
     CU_ASSERT_FATAL(ret == 0);
 
-    /* Add a site and a mutation */
+    /* Add a site and a mutation with metadata, above node 4 */
     if (alphabet == ALPHABET_BINARY) {
-        ret = tsk_site_table_add_row(&tables->sites, 0.1, "0", 1, NULL, 0);
+        ret = tsk_site_table_add_row(&tables->sites, 0.1, "0", 1, site_metadata, 3);
         CU_ASSERT_FATAL(ret >= 0);
         ret = tsk_mutation_table_add_row(
-            &tables->mutations, 0, 0, -1, 0.0, "1", 1, NULL, 0);
+            &tables->mutations, 0, 4, -1, 2.5, "1", 1, mut_metadata, 5);
         CU_ASSERT_FATAL(ret >= 0);
     } else if (alphabet == ALPHABET_NUCLEOTIDE) {
-        ret = tsk_site_table_add_row(&tables->sites, 0.1, "A", 1, NULL, 0);
+        ret = tsk_site_table_add_row(&tables->sites, 0.1, "A", 1, site_metadata, 3);
         CU_ASSERT_FATAL(ret >= 0);
         ret = tsk_mutation_table_add_row(
-            &tables->mutations, 0, 0, -1, 0.0, "C", 1, NULL, 0);
+            &tables->mutations, 0, 4, -1, 2.5, "C", 1, mut_metadata, 5);
         CU_ASSERT_FATAL(ret >= 0);
     }
 
@@ -864,13 +867,19 @@ parse_text_int64(char *ds, tsk_size_t n)
 static void
 verify_slim_mutation_ids(int64_t *mut_ids, size_t mut_ids_length, int64_t min_mut_id)
 {
-    int j;
+    int j = 0;
     CU_ASSERT_FATAL(mut_ids_length > 0);
     qsort(mut_ids, mut_ids_length, sizeof(*mut_ids), &cmp_int64);
-    CU_ASSERT_EQUAL_FATAL(mut_ids[0], min_mut_id);
+    // skip any -1s we might have put in
+    while (mut_ids[j] < 0) {
+        j++;
+    }
+    CU_ASSERT_EQUAL_FATAL(mut_ids[j], min_mut_id);
+    j++;
     if (mut_ids_length > 1) {
-        for (j = 1; j < mut_ids_length; j++) {
+        while (j < mut_ids_length) {
             CU_ASSERT_EQUAL_FATAL(mut_ids[j] - mut_ids[j - 1], 1);
+            j++;
         }
     }
 }
@@ -1056,6 +1065,265 @@ test_mutgen_slim_mutation_large_values(void)
     gsl_rng_free(rng);
 }
 
+static void
+verify_slim_v6_derived_states(tsk_site_table_t *sites, tsk_mutation_table_t *mutations,
+    int64_t *all_mut_ids, size_t all_mut_ids_length, int skip_mut)
+{
+    int j, k, s;
+    size_t len, parent_len;
+    char *ds;
+    int64_t mut_id;
+    char *parent_allele;
+    CU_ASSERT_FATAL(mutations->num_rows <= all_mut_ids_length);
+    for (j = 0; j < mutations->num_rows; j++) {
+        ds = mutations->derived_state + mutations->derived_state_offset[j];
+        len = (mutations->derived_state_offset[j + 1]
+               - mutations->derived_state_offset[j]);
+        if (j == skip_mut) {
+            parent_len = 0;
+            all_mut_ids[j] = -1;
+        } else {
+            k = mutations->parent[j];
+            if (k == TSK_NULL) {
+                s = mutations->site[j];
+                parent_len = (sites->ancestral_state_offset[s + 1]
+                              - sites->ancestral_state_offset[s]);
+                parent_allele
+                    = sites->ancestral_state + sites->ancestral_state_offset[s];
+            } else {
+                parent_len = (mutations->derived_state_offset[k + 1]
+                              - mutations->derived_state_offset[k]);
+                parent_allele
+                    = mutations->derived_state + mutations->derived_state_offset[k];
+            }
+            CU_ASSERT_FATAL(len > parent_len);
+            if (parent_len > 0) {
+                CU_ASSERT_EQUAL_FATAL(memcmp(ds, parent_allele, parent_len), 0);
+                CU_ASSERT_EQUAL_FATAL((ds + parent_len)[0], 44); // 44 is ',' in ascii
+                parent_len++;
+            }
+            CU_ASSERT_FATAL(len > parent_len);
+            mut_id = parse_text_int64(ds + parent_len, len - parent_len);
+            all_mut_ids[j] = mut_id;
+        }
+    }
+}
+
+static void
+verify_slim_v6_metadata(tsk_site_table_t *sites, tsk_mutation_table_t *mutations,
+    int64_t *all_mut_ids, size_t all_mut_ids_length, int skip_mut)
+{
+    int j, k, s;
+    size_t len, parent_len;
+    int64_t mut_id;
+    char *md, *pmd;
+    CU_ASSERT_FATAL(mutations->num_rows <= all_mut_ids_length);
+    for (j = 0; j < mutations->num_rows; j++) {
+        md = mutations->metadata + mutations->metadata_offset[j];
+        len = (mutations->metadata_offset[j + 1] - mutations->metadata_offset[j]);
+        if (j == skip_mut) {
+            parent_len = 0;
+            all_mut_ids[j] = -1;
+        } else {
+            k = mutations->parent[j];
+            if (k == TSK_NULL) {
+                s = mutations->site[j];
+                pmd = sites->metadata + sites->metadata_offset[s];
+                parent_len = (sites->metadata_offset[s + 1] - sites->metadata_offset[s]);
+            } else {
+                pmd = mutations->metadata + mutations->metadata_offset[k];
+                parent_len = (mutations->metadata_offset[k + 1]
+                              - mutations->metadata_offset[k]);
+            }
+            CU_ASSERT_FATAL(len > parent_len);
+            if (parent_len > 0) {
+                CU_ASSERT_EQUAL_FATAL(memcmp(md, pmd, parent_len), 0);
+            }
+            CU_ASSERT_FATAL(len == parent_len + sizeof(int64_t));
+            mut_id = (int64_t) md[parent_len];
+            all_mut_ids[j] = mut_id;
+        }
+    }
+}
+
+static void
+test_mutgen_slim_v6_mutations(void)
+{
+    int ret = 0;
+    int j;
+    mutgen_t mutgen;
+    gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);
+    tsk_table_collection_t tables;
+    mutation_model_t mut_model;
+    int64_t *all_mut_ids;
+    int64_t next_mutation_id = 23;
+
+    CU_ASSERT_FATAL(rng != NULL);
+    ret = slim_v6_mutation_model_alloc(&mut_model, next_mutation_id, 0);
+    CU_ASSERT_EQUAL_FATAL(ret, 0);
+
+    ret = tsk_table_collection_init(&tables, 0);
+    CU_ASSERT_EQUAL_FATAL(ret, 0);
+    insert_single_tree(&tables, -1);
+
+    ret = mutgen_alloc(&mutgen, rng, &tables, &mut_model, 0);
+    CU_ASSERT_EQUAL_FATAL(ret, 0);
+    ret = mutgen_set_rate(&mutgen, 1);
+    CU_ASSERT_EQUAL_FATAL(ret, 0);
+    ret = mutgen_generate(&mutgen, MSP_DISCRETE_SITES);
+    CU_ASSERT_EQUAL_FATAL(ret, 0);
+    CU_ASSERT_FATAL(tables.mutations.num_rows > 0);
+
+    // should have empty ancestral states
+    for (j = 0; j < tables.sites.num_rows; j++) {
+        CU_ASSERT_EQUAL_FATAL(tables.sites.ancestral_state_offset[j], 0);
+    }
+    // check that derived states append unique integers,
+    // counting up from next_mutation_id
+    // in both derived state (text) and metadata (binary)
+    all_mut_ids = malloc(tables.mutations.num_rows * sizeof(int64_t));
+    CU_ASSERT_FATAL(all_mut_ids != NULL);
+    // derived state
+    verify_slim_v6_derived_states(
+        &tables.sites, &tables.mutations, all_mut_ids, tables.mutations.num_rows, -1);
+    verify_slim_mutation_ids(all_mut_ids, tables.mutations.num_rows, next_mutation_id);
+    // metadata
+    verify_slim_v6_metadata(
+        &tables.sites, &tables.mutations, all_mut_ids, tables.mutations.num_rows, -1);
+    verify_slim_mutation_ids(all_mut_ids, tables.mutations.num_rows, next_mutation_id);
+
+    mutgen_print_state(&mutgen, _devnull);
+
+    mutgen_free(&mutgen);
+    free(all_mut_ids);
+    mutation_model_free(&mut_model);
+    tsk_table_collection_free(&tables);
+    gsl_rng_free(rng);
+}
+
+static void
+test_mutgen_slim_v6_mutations_keeps(void)
+{
+    int ret = 0;
+    int j;
+    mutgen_t mutgen;
+    gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);
+    tsk_table_collection_t tables, orig_tables;
+    mutation_model_t mut_model;
+    int64_t *all_mut_ids;
+    int64_t next_mutation_id = 23;
+
+    CU_ASSERT_FATAL(rng != NULL);
+    ret = slim_v6_mutation_model_alloc(&mut_model, next_mutation_id, 0);
+    CU_ASSERT_EQUAL_FATAL(ret, 0);
+
+    ret = tsk_table_collection_init(&tables, 0);
+    CU_ASSERT_EQUAL_FATAL(ret, 0);
+    insert_single_tree(&tables, ALPHABET_NUCLEOTIDE);
+    // this should have a site and a mutation and only one position
+    // so we know future mutations will hit it
+    CU_ASSERT_EQUAL_FATAL(tables.sites.num_rows, 1);
+    // change the position of that site to 0.0
+    tables.sites.position[0] = 0.0;
+    CU_ASSERT_FATAL(tables.mutations.num_rows > 0);
+    CU_ASSERT_EQUAL_FATAL(tables.sequence_length, 1.0);
+
+    ret = tsk_table_collection_copy(&tables, &orig_tables, 0);
+    CU_ASSERT_EQUAL_FATAL(ret, 0);
+    CU_ASSERT_EQUAL_FATAL(tables.sequence_length, 1.0);
+
+    ret = mutgen_alloc(&mutgen, rng, &tables, &mut_model, 0);
+    CU_ASSERT_EQUAL_FATAL(ret, 0);
+    ret = mutgen_set_rate(&mutgen, 1);
+    CU_ASSERT_EQUAL_FATAL(ret, 0);
+    ret = mutgen_generate(&mutgen, MSP_DISCRETE_SITES | MSP_KEEP_SITES);
+    CU_ASSERT_EQUAL_FATAL(ret, 0);
+    CU_ASSERT_FATAL(tables.mutations.num_rows > orig_tables.mutations.num_rows);
+    CU_ASSERT_EQUAL_FATAL(tables.sites.num_rows, 1);
+
+    // check site metadata has been preserved
+    ret = tsk_site_table_equals(&tables.sites, &orig_tables.sites, 0);
+    CU_ASSERT_FATAL(ret);
+
+    // check for derived states and metadata that each mutation
+    // is something appended onto its parent's
+    all_mut_ids = malloc(tables.mutations.num_rows * sizeof(int64_t));
+    CU_ASSERT_FATAL(all_mut_ids != NULL);
+    for (j = 0; j < tables.mutations.num_rows; j++) {
+        // find the mutation that's kept; it doesn't follow the rules
+        if (tables.mutations.time[j] == orig_tables.mutations.time[0])
+            break;
+    }
+    verify_slim_v6_derived_states(
+        &tables.sites, &tables.mutations, all_mut_ids, tables.mutations.num_rows, j);
+    verify_slim_mutation_ids(all_mut_ids, tables.mutations.num_rows, next_mutation_id);
+    verify_slim_v6_metadata(
+        &tables.sites, &tables.mutations, all_mut_ids, tables.mutations.num_rows, j);
+    verify_slim_mutation_ids(all_mut_ids, tables.mutations.num_rows, next_mutation_id);
+
+    free(all_mut_ids);
+    mutgen_free(&mutgen);
+    mutation_model_free(&mut_model);
+    tsk_table_collection_free(&tables);
+    tsk_table_collection_free(&orig_tables);
+    gsl_rng_free(rng);
+}
+
+static void
+test_mutgen_slim_v6_mutation_large_values(void)
+{
+    int ret = 0;
+    mutgen_t mutgen;
+    tsk_mutation_t mut;
+    gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);
+    tsk_table_collection_t tables;
+    mutation_model_t mut_model;
+    char value[21]; /* longest 64 bit int has 20 digits */
+    int64_t next_mutation_id;
+
+    CU_ASSERT_FATAL(rng != NULL);
+    ret = tsk_table_collection_init(&tables, 0);
+    CU_ASSERT_EQUAL_FATAL(ret, 0);
+    insert_single_tree(&tables, -1);
+
+    /* Trying to generate mutations that overflow raises an error */
+    ret = slim_v6_mutation_model_alloc(&mut_model, INT64_MAX, 0);
+    CU_ASSERT_EQUAL_FATAL(ret, 0);
+    ret = mutgen_alloc(&mutgen, rng, &tables, &mut_model, 0);
+    CU_ASSERT_EQUAL_FATAL(ret, 0);
+    ret = mutgen_set_rate(&mutgen, 2);
+    CU_ASSERT_EQUAL_FATAL(ret, 0);
+    ret = mutgen_generate(&mutgen, MSP_DISCRETE_SITES);
+    CU_ASSERT_EQUAL_FATAL(ret, MSP_ERR_MUTATION_ID_OVERFLOW);
+    mutgen_free(&mutgen);
+    mutation_model_free(&mut_model);
+
+    tsk_mutation_table_clear(&tables.mutations);
+    tsk_site_table_clear(&tables.sites);
+    /* Try out with a large value that doesn't hit the ceiling */
+    next_mutation_id = INT64_MAX - 100;
+    ret = slim_v6_mutation_model_alloc(&mut_model, next_mutation_id, 0);
+    CU_ASSERT_EQUAL_FATAL(ret, 0);
+    ret = mutgen_alloc(&mutgen, rng, &tables, &mut_model, 0);
+    CU_ASSERT_EQUAL_FATAL(ret, 0);
+    ret = mutgen_set_rate(&mutgen, 2);
+    CU_ASSERT_EQUAL_FATAL(ret, 0);
+    ret = mutgen_generate(&mutgen, MSP_DISCRETE_SITES);
+    CU_ASSERT_EQUAL_FATAL(ret, 0);
+    CU_ASSERT_TRUE(tables.mutations.num_rows > 10);
+    mutgen_free(&mutgen);
+    mutation_model_free(&mut_model);
+
+    ret = tsk_mutation_table_get_row(&tables.mutations, 0, &mut);
+    CU_ASSERT_EQUAL_FATAL(ret, 0);
+    CU_ASSERT_EQUAL(mut.derived_state_length, 19);
+    sprintf(value, "%" PRId64, next_mutation_id);
+    CU_ASSERT_NSTRING_EQUAL(value, mut.derived_state, mut.derived_state_length);
+
+    tsk_table_collection_free(&tables);
+    gsl_rng_free(rng);
+}
+
 static void
 test_mutgen_infinite_alleles(void)
 {
@@ -1306,6 +1574,33 @@ test_slim_mutation_model_properties(void)
     mutation_model_free(&model);
 }
 
+static void
+test_slim_v6_mutation_model_errors(void)
+{
+    int ret;
+    mutation_model_t model;
+    int64_t next_mutation_id = 0;
+
+    next_mutation_id--;
+    ret = slim_v6_mutation_model_alloc(&model, next_mutation_id, 0);
+    CU_ASSERT_EQUAL_FATAL(ret, MSP_ERR_BAD_SLIM_PARAMETERS);
+    mutation_model_free(&model);
+}
+
+static void
+test_slim_v6_mutation_model_properties(void)
+{
+    int ret;
+    mutation_model_t model;
+    int64_t next_mutation_id = 2;
+
+    ret = slim_v6_mutation_model_alloc(&model, next_mutation_id, 0);
+    CU_ASSERT_EQUAL_FATAL(ret, 0);
+    CU_ASSERT_EQUAL_FATAL(model.params.slim_mutator.next_mutation_id, 2);
+
+    mutation_model_free(&model);
+}
+
 int
 main(int argc, char **argv)
 {
@@ -1331,6 +1626,10 @@ main(int argc, char **argv)
         { "test_mutgen_slim_mutations", test_mutgen_slim_mutations },
         { "test_mutgen_slim_mutation_large_values",
             test_mutgen_slim_mutation_large_values },
+        { "test_mutgen_slim_v6_mutations", test_mutgen_slim_v6_mutations },
+        { "test_mutgen_slim_v6_mutations_keeps", test_mutgen_slim_v6_mutations_keeps },
+        { "test_mutgen_slim_v6_mutation_large_values",
+            test_mutgen_slim_v6_mutation_large_values },
         { "test_mutgen_infinite_alleles", test_mutgen_infinite_alleles },
         { "test_mutgen_infinite_alleles_large_values",
             test_mutgen_infinite_alleles_large_values },
@@ -1339,6 +1638,9 @@ main(int argc, char **argv)
             test_matrix_mutation_model_properties },
         { "test_slim_mutation_model_errors", test_slim_mutation_model_errors },
         { "test_slim_mutation_model_properties", test_slim_mutation_model_properties },
+        { "test_slim_v6_mutation_model_errors", test_slim_v6_mutation_model_errors },
+        { "test_slim_v6_mutation_model_properties",
+            test_slim_v6_mutation_model_properties },
         CU_TEST_INFO_NULL,
     };
 
diff --git a/msprime/__init__.py b/msprime/__init__.py
index 11acef934..1878d7661 100644
--- a/msprime/__init__.py
+++ b/msprime/__init__.py
@@ -81,6 +81,7 @@
     NUCLEOTIDES,
     PAM,
     SLiMMutationModel,
+    SLiMv6MutationModel,
     sim_mutations,
 )
 
@@ -151,6 +152,7 @@
     "RateMap",
     "RecombinationMap",
     "SLiMMutationModel",
+    "SLiMv6MutationModel",
     "Sample",
     "SampleSet",
     "SimpleBottleneck",
diff --git a/msprime/_msprimemodule.c b/msprime/_msprimemodule.c
index f3cf82d94..e9b38d673 100644
--- a/msprime/_msprimemodule.c
+++ b/msprime/_msprimemodule.c
@@ -65,6 +65,11 @@ typedef struct {
     mutation_model_t *mutation_model;
 } SLiMMutationModel;
 
+typedef struct {
+    PyObject_HEAD
+    mutation_model_t *mutation_model;
+} SLiMv6MutationModel;
+
 typedef struct {
     PyObject_HEAD
     mutation_model_t *mutation_model;
@@ -796,6 +801,98 @@ static PyTypeObject SLiMMutationModelType = {
     .tp_new = PyType_GenericNew,
 };
 
+/*===================================================================
+ * SLiM mutation model, v6
+ *===================================================================
+ */
+
+static int
+SLiMv6MutationModel_check_state(SLiMv6MutationModel *self)
+{
+    int ret = 0;
+    if (self->mutation_model == NULL) {
+        PyErr_SetString(PyExc_SystemError, "SLiMv6MutationModel not initialised");
+        ret = -1;
+    }
+    return ret;
+}
+
+static void
+SLiMv6MutationModel_dealloc(SLiMv6MutationModel *self)
+{
+    if (self->mutation_model != NULL) {
+        mutation_model_free(self->mutation_model);
+        PyMem_Free(self->mutation_model);
+        self->mutation_model = NULL;
+    }
+    Py_TYPE(self)->tp_free((PyObject *) self);
+}
+
+static int
+SLiMv6MutationModel_init(SLiMv6MutationModel *self, PyObject *args, PyObject *kwds)
+{
+    int ret = -1;
+    int err;
+    static char *kwlist[] = { "next_id", "block_size", NULL };
+    long long next_id = 0;
+    Py_ssize_t block_size = 0;
+
+    self->mutation_model = NULL;
+    if (!PyArg_ParseTupleAndKeywords(args, kwds, "|Ln", kwlist, &next_id, &block_size)) {
+        goto out;
+    }
+
+    /* Note: it's important we zero out mutation_model here because
+     * we can error before we can mutation_model_alloc, leaving the
+     * object in an uninitialised state */
+    self->mutation_model = PyMem_Calloc(1, sizeof(*self->mutation_model));
+    if (self->mutation_model == NULL) {
+        PyErr_NoMemory();
+        goto out;
+    }
+    err = slim_v6_mutation_model_alloc(
+        self->mutation_model, (int64_t) next_id, (size_t) block_size);
+    if (err != 0) {
+        handle_library_error(err);
+        goto out;
+    }
+    ret = 0;
+out:
+    return ret;
+}
+
+static PyObject *
+SLiMv6MutationModel_get_next_id(SLiMv6MutationModel *self, void *closure)
+{
+    slim_mutator_t *params;
+    PyObject *ret = NULL;
+
+    if (SLiMv6MutationModel_check_state(self) != 0) {
+        goto out;
+    }
+    params = &self->mutation_model->params.slim_mutator;
+    ret = Py_BuildValue("L", params->next_mutation_id);
+out:
+    return ret;
+}
+
+static PyGetSetDef SLiMv6MutationModel_getsetters[] = {
+    { "next_id", (getter) SLiMv6MutationModel_get_next_id, NULL,
+        "Return the next mutation id" },
+    { NULL } /* Sentinel */
+};
+
+static PyTypeObject SLiMv6MutationModelType = {
+    .tp_name = "_msprime.SLiMv6MutationModel",
+    .tp_basicsize = sizeof(SLiMv6MutationModel),
+    .tp_dealloc = (destructor) SLiMv6MutationModel_dealloc,
+    .tp_flags = Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE,
+    .tp_doc = "SLiMv6MutationModel objects",
+    .tp_getset = SLiMv6MutationModel_getsetters,
+    .tp_init = (initproc) SLiMv6MutationModel_init,
+    .tp_new = PyType_GenericNew,
+};
+
 /*===================================================================
  * Infinite alleles mutation model
  *===================================================================
@@ -3014,6 +3111,7 @@ parse_mutation_model(PyObject *py_model)
     mutation_model_t *model = NULL;
     MatrixMutationModel *matrix_mutation_model = NULL;
     SLiMMutationModel *slim_mutation_model = NULL;
+    SLiMv6MutationModel *slim_v6_mutation_model = NULL;
     InfiniteAllelesMutationModel *infinite_alleles_model = NULL;
 
     if (PyObject_TypeCheck(py_model, &MatrixMutationModelType)) {
@@ -3028,6 +3126,12 @@ parse_mutation_model(PyObject *py_model)
             goto out;
         }
         model = slim_mutation_model->mutation_model;
+    } else if (PyObject_TypeCheck(py_model, &SLiMv6MutationModelType)) {
+        slim_v6_mutation_model = (SLiMv6MutationModel *) py_model;
+        if (SLiMv6MutationModel_check_state(slim_v6_mutation_model) != 0) {
+            goto out;
+        }
+        model = slim_v6_mutation_model->mutation_model;
     } else if (PyObject_TypeCheck(py_model, &InfiniteAllelesMutationModelType)) {
         infinite_alleles_model = (InfiniteAllelesMutationModel *) py_model;
         if (InfiniteAllelesMutationModel_check_state(infinite_alleles_model) != 0) {
@@ -3037,7 +3141,7 @@ parse_mutation_model(PyObject *py_model)
     } else {
         PyErr_SetString(PyExc_TypeError,
             "model must be an instance of MatrixMutationModel, "
-            "SLiMMutationModel or InfiniteAllelesMutationModel.");
+            "SLiMMutationModel, SLiMv6MutationModel or InfiniteAllelesMutationModel.");
         goto out;
     }
 out:
@@ -3256,6 +3360,15 @@ PyInit__msprime(void)
     Py_INCREF(&SLiMMutationModelType);
     PyModule_AddObject(module, "SLiMMutationModel", (PyObject *) &SLiMMutationModelType);
 
+    /* SLiMv6MutationModel type */
+    SLiMv6MutationModelType.tp_base = &BaseMutationModelType;
+    if (PyType_Ready(&SLiMv6MutationModelType) < 0) {
+        return NULL;
+    }
+    Py_INCREF(&SLiMv6MutationModelType);
+    PyModule_AddObject(
+        module, "SLiMv6MutationModel", (PyObject *) &SLiMv6MutationModelType);
+
     /* InfiniteAllelesMutationModel type */
     InfiniteAllelesMutationModelType.tp_base = &BaseMutationModelType;
     if (PyType_Ready(&InfiniteAllelesMutationModelType) < 0) {
diff --git a/msprime/mutations.py b/msprime/mutations.py
index 55c58305d..4d7fbdfa8 100644
--- a/msprime/mutations.py
+++ b/msprime/mutations.py
@@ -108,7 +108,9 @@ def __str__(self):
 
 class SLiMMutationModel(_msprime.SLiMMutationModel, MutationModel):
     """
-    An infinite-alleles model of mutation producing "SLiM-style" mutations.
+    An infinite-alleles model of mutation producing "SLiM-style" mutations
+    for versions of SLiM before 6.0. For SLiM version 6.0 and beyond,
+    use {class}`.SLiMv6MutationModel`.
 
     To agree with mutations produced by SLiM, the ancestral state of each new
     site is set to the empty string, and each derived state is produced by
@@ -144,6 +146,48 @@ def __str__(self):
         )
 
 
+class SLiMv6MutationModel(_msprime.SLiMv6MutationModel, MutationModel):
+    """
+    An infinite-alleles model of mutation producing "SLiM-style" mutations
+    for versions of SLiM 6.0 and beyond. For older versions of SLiM,
+    use {class}`.SLiMMutationModel`.
+
+    To agree with mutations produced by SLiM, the ancestral state of each new
+    site is set to the empty string, and each derived state is produced by
+    appending the "next allele" to the previous state.  The result is that each
+    allele is a comma-separated string of all mutations that have occurred up
+    to the root. Alleles are numeric IDs, starting with ``next_id``.
+
+    The same list of numeric IDs is stored in metadata, in binary, as int64s.
+    It is actually these IDs in metadata that SLiM uses when reading in a tree
+    sequence, not the text-based version in the derived state.
+
+    This model works as following: for each mutation, if the next ID is ``n``,
+    then it sets the derived state to be the string ``"n"`` appended to its
+    inherited state (which is the derived state of the parent mutation, if any,
+    or the ancestral state otherwise), with an intervening ``","`` if the inherited
+    state is nonempty. Also, it sets the metadata to be the 64-bit integer n
+    appended to the bytes of the inherited metadata (which similarly comes from
+    the parent mutation, if any, or the site).  So, if this mutation model is
+    applied to a tree sequence with existing mutations of a different sort (especially
+    if the existing sites have nonempty metadata) the result may be surprising:
+    see {ref}`sec_mutations_mutation_slim_mutations_stacking`.
+
+    To update your code to use this model instead of {class}`.SLiMMutationModel`,
+    remove the ``type`` and ``slim_generation`` arguments, and follow
+    up {func}`.sim_mutations` with a call to ``pyslim.add_mutation_metadata``.
+
+    :param int next_id: The nonnegative integer to start assigning alleles from.
+        (default: 0)
+    :param int block_size: The block size for allocating derived states.
+        You do not need to change this unless you get an "out of memory" error
+        due to a very large number of stacked mutations.
+    """
+
+    def __str__(self):
+        return f"Mutation model for SLiM mutations, v6+. Next ID: {self.next_id}\n"
+
+
 # NOTE: we use a hacky workaround here for documenting the next_allele instance
 # variable because of the nasty sphinx bug:
 # https://github.com/sphinx-doc/sphinx/issues/2549
@@ -1512,7 +1556,8 @@ def mutate(
     "jc69": JC69,
     "blosum62": BLOSUM62,
     "pam": PAM,
-    # "slim": SLiMMutationModel(), Needs type argument so can't be string init'd
+    # "slim": SLiMMutationModel(), or SLiMv6MutationModel(), not clear
+    # and maybe needs the next_id argument
     # "hky": HKY(), Needs kappa argument
     # "f84": F84(), Needs kappa argument
     # "gtr": GTR(), Needs relative_rates argument
diff --git a/tests/test_mutations.py b/tests/test_mutations.py
index 569b5d9a2..1deae6608 100644
--- a/tests/test_mutations.py
+++ b/tests/test_mutations.py
@@ -1551,7 +1551,57 @@ class SlimMetadata:
     nucleotide: int
 
 
-class TestSLiMMutationModel:
+class SLiMModelMixin:
+    """
+    Common tests for both SLiM mutation models
+    """
+
+    def test_binary_n_4_low_rate(self):
+        ts = msprime.sim_ancestry(4, sequence_length=10, random_seed=5)
+        mts = self.run_mutate(ts, rate=0.1, random_seed=23)
+        assert mts.num_mutations > 1
+        self.validate_slim_mutations(mts)
+
+    def test_binary_n_4_high_rate(self):
+        ts = msprime.sim_ancestry(4, sequence_length=2, random_seed=5)
+        mts = self.run_mutate(ts, rate=2.0, random_seed=23)
+        assert mts.num_mutations > 10
+        self.validate_slim_mutations(mts)
+
+    def test_binary_n_8_low_rate(self):
+        ts = msprime.sim_ancestry(8, sequence_length=10, random_seed=50)
+        mts = self.run_mutate(ts, rate=0.1, random_seed=342)
+        assert mts.num_mutations > 1
+        self.validate_slim_mutations(mts)
+
+    def test_binary_n_8_high_rate(self):
+        ts = msprime.sim_ancestry(8, sequence_length=10, random_seed=5)
+        mts = self.run_mutate(ts, rate=2.0, random_seed=23)
+        assert mts.num_mutations > 10
+        self.validate_slim_mutations(mts)
+
+    def test_binary_incomplete_trees(self):
+        ts = msprime.sim_ancestry(8, sequence_length=5, random_seed=50, end_time=0.1)
+        assert ts.first().num_roots > 1
+        mts = self.run_mutate(ts, rate=2.0, random_seed=23)
+        assert mts.num_mutations > 10
+        self.validate_slim_mutations(mts)
+
+    def test_binary_many_trees(self):
+        ts = msprime.sim_ancestry(
+            8,
+            sequence_length=5,
+            recombination_rate=5,
+            random_seed=50,
+            discrete_genome=False,
+        )
+        assert ts.num_trees > 20
+        mts = self.run_mutate(ts, rate=2.0, random_seed=23)
+        assert mts.num_mutations > 10
+        self.validate_slim_mutations(mts)
+
+
+class TestSLiMMutationModel(SLiMModelMixin):
     """
     Tests for the SLiM mutation generator.
     """
@@ -1619,18 +1669,10 @@ def run_mutate(
 
         t1 = mts1.dump_tables()
         t2 = mts2.dump_tables()
-        assert t1.sites == t2.sites
         # Drop the mutation metadata - we're validating that elsewhere and
         # it's not worth complicating the Python generator with it.
-        t1.mutations.set_columns(
-            site=t1.mutations.site,
-            node=t1.mutations.node,
-            parent=t1.mutations.parent,
-            time=t1.mutations.time,
-            derived_state=t1.mutations.derived_state,
-            derived_state_offset=t1.mutations.derived_state_offset,
-        )
-        assert t1.mutations == t2.mutations
+        t1.sites.assert_equals(t2.sites)
+        t1.mutations.assert_equals(t2.mutations, ignore_metadata=True)
         return mts1
 
     def test_slim_mutation_type(self):
@@ -1651,49 +1693,194 @@ def test_slim_generation(self):
             assert mts.num_mutations > 10
             self.validate_slim_mutations(mts, slim_generation=slim_generation)
 
-    def test_binary_n_4_low_rate(self):
-        ts = msprime.sim_ancestry(4, sequence_length=10, random_seed=5)
-        mts = self.run_mutate(ts, rate=0.1, random_seed=23)
-        assert mts.num_mutations > 1
-        self.validate_slim_mutations(mts)
 
-    def test_binary_n_4_high_rate(self):
-        ts = msprime.sim_ancestry(4, sequence_length=2, random_seed=5)
-        mts = self.run_mutate(ts, rate=2.0, random_seed=23)
-        assert mts.num_mutations > 10
-        self.validate_slim_mutations(mts)
+class TestSLiMv6MutationModel(SLiMModelMixin):
+    """
+    Tests for the SLiMv6 mutation generator.
+    """
 
-    def test_binary_n_8_low_rate(self):
-        ts = msprime.sim_ancestry(8, sequence_length=10, random_seed=50)
-        mts = self.run_mutate(ts, rate=0.1, random_seed=342)
-        assert mts.num_mutations > 1
-        self.validate_slim_mutations(mts)
+    metadata_schema = tskit.MetadataSchema(
+        {
+            "codec": "struct",
+            "type": "object",
+            "properties": {
+                "derived_states": {
+                    "items": {
+                        "binaryFormat": "q",
+                        "type": "number",
+                    },
+                    "noLengthEncodingExhaustBuffer": True,
+                    "type": "array",
+                }
+            },
+            "required": ["derived_states"],
+            "additionalProperties": False,
+        }
+    )
 
-    def test_binary_n_8_high_rate(self):
-        ts = msprime.sim_ancestry(8, sequence_length=10, random_seed=5)
-        mts = self.run_mutate(ts, rate=2.0, random_seed=23)
-        assert mts.num_mutations > 10
-        self.validate_slim_mutations(mts)
+    def validate_slim_mutations(self, ts):
+        # slim alleles should be lists of integers
+        # and ancestral states the empty string
+        # while metadata should be lists of integers
+        t = ts.dump_tables()
+        t.mutations.metadata_schema = self.metadata_schema
+        ts = t.tree_sequence()
+        for site in ts.sites():
+            assert site.ancestral_state == ""
+            alleles = {}
+            for mutation in site.mutations:
+                a = list(map(int, mutation.derived_state.split(",")))
+                md = mutation.metadata["derived_states"]
+                assert len(a) == len(md)
+                for x, y in zip(a, md):
+                    assert x == y
+                alleles[mutation.id] = a
+                if mutation.parent == tskit.NULL:
+                    assert len(a) == 1
+                else:
+                    parent_allele = alleles[mutation.parent]
+                    assert a[:-1] == parent_allele
 
-    def test_binary_incomplete_trees(self):
-        ts = msprime.sim_ancestry(8, sequence_length=5, random_seed=50, end_time=0.1)
-        assert ts.first().num_roots > 1
-        mts = self.run_mutate(ts, rate=2.0, random_seed=23)
-        assert mts.num_mutations > 10
-        self.validate_slim_mutations(mts)
+    def run_mutate(
+        self,
+        ts,
+        rate=1,
+        random_seed=42,
+        mutation_id=0,
+        keep=True,
+    ):
+        model = msprime.SLiMv6MutationModel(next_id=mutation_id)
+        mts1 = msprime.sim_mutations(
+            ts,
+            rate=rate,
+            random_seed=random_seed,
+            model=model,
+            discrete_genome=True,
+            keep=keep,
+        )
+        if ts.num_mutations == 0 or not keep:
+            assert mts1.num_mutations == model.next_id
 
-    def test_binary_many_trees(self):
-        ts = msprime.sim_ancestry(
-            8,
-            sequence_length=5,
-            recombination_rate=5,
-            random_seed=50,
-            discrete_genome=False,
+        model = PythonSLiMv6MutationModel(next_id=mutation_id)
+        mts2 = py_sim_mutations(
+            ts,
+            rate=rate,
+            random_seed=random_seed,
+            model=model,
+            discrete_genome=True,
+            keep=keep,
         )
-        assert ts.num_trees > 20
+
+        t1 = mts1.dump_tables()
+        t2 = mts2.dump_tables()
+        # Drop the mutation metadata - we're validating that elsewhere and
+        # it's not worth complicating the Python generator with it.
+        t1.sites.assert_equals(t2.sites)
+        t1.mutations.assert_equals(t2.mutations, ignore_metadata=True)
+        return mts1
+
+    metadata_schema_with_extra = tskit.MetadataSchema(
+        {
+            "codec": "struct",
+            "type": "object",
+            "properties": {
+                "tag": {
+                    "index": 0,
+                    "type": "string",
+                    "binaryFormat": "24s",
+                    "nullTerminated": True,
+                },
+                "derived_states": {
+                    "index": 1,
+                    "items": {
+                        "binaryFormat": "q",
+                        "type": "number",
+                    },
+                    "noLengthEncodingExhaustBuffer": True,
+                    "type": "array",
+                },
+            },
+            "required": ["tag", "derived_states"],
+            "additionalProperties": False,
+        }
+    )
+
+    def test_keep(self):
+        ts = msprime.sim_ancestry(8, sequence_length=10, random_seed=5)
         mts = self.run_mutate(ts, rate=2.0, random_seed=23)
-        assert mts.num_mutations > 10
-        self.validate_slim_mutations(mts)
+        mt = mts.dump_tables()
+        # there can be weird interactions with site metadata, so add some
+        mt.sites.clear()
+        mt.sites.metadata_schema = self.metadata_schema_with_extra
+        for s in mts.sites():
+            md = {
+                "tag": f"lmnop{s.id}",
+                "derived_states": [],
+            }
+            mt.sites.append(s.replace(metadata=md))
+        mt.mutations.clear()
+        mt.mutations.metadata_schema = self.metadata_schema_with_extra
+        for m in mts.mutations():
+            md = {
+                "tag": f"lmnop{m.site}",
+                "derived_states": [int(x) for x in m.derived_state.split(",")],
+            }
+            mt.mutations.append(m.replace(metadata=md))
+        mts = mt.tree_sequence()
+        num_muts = {s.position: len(s.mutations) for s in mts.sites()}
+        mut_id_start = mts.num_mutations + 1
+        model = msprime.SLiMv6MutationModel(next_id=mut_id_start)
+        mmts = msprime.sim_mutations(
+            mts,
+            rate=2.0,
+            random_seed=1234,
+            model=model,
+            discrete_genome=True,
+            keep=True,
+        )
+        # check we've stacked on top of some previous ones
+        more_hits = False
+        for s in mmts.sites():
+            if s.position in num_muts:
+                if len(s.mutations) > num_muts[s.position]:
+                    more_hits = True
+                    break
+        assert more_hits
+        for mut in mmts.mutations():
+            site = mmts.site(mut.site)
+            mp = mut.parent
+            # mutations we put down the first time won't have inherited the state
+            # from mutations we put down the second time, so traverse up until
+            # we find the correct one we should have inherited from
+            mut_id = mut.metadata["derived_states"][-1]
+            if mut_id < mut_id_start:
+                while (
+                    mp >= 0
+                    and mmts.mutation(mp).metadata["derived_states"][-1] >= mut_id_start
+                ):
+                    mp = mmts.mutation(mp).parent
+            if mp == -1:
+                parent_allele = site.ancestral_state
+                parent_metadata = site.metadata
+            else:
+                parent_allele = mmts.mutation(mp).derived_state
+                parent_metadata = mmts.mutation(mp).metadata
+            n = len(parent_allele)
+            assert n < len(mut.derived_state)
+            assert parent_allele == mut.derived_state[:n]
+            if n > 0:
+                assert mut.derived_state[n] == ","
+                n += 1
+            assert mut_id == int(mut.derived_state[n:])
+            assert parent_metadata["tag"] == mut.metadata["tag"]
+            n = len(parent_metadata["derived_states"])
+            assert n + 1 == len(mut.metadata["derived_states"])
+            if n > 0:
+                assert (
+                    parent_metadata["derived_states"]
+                    == mut.metadata["derived_states"][:n]
+                )
+            assert mut.metadata["derived_states"][-1] == mut_id
 
 
 class TestInfiniteAllelesMutationModel:
@@ -1984,8 +2171,7 @@ def transition_allele(self, rng, current_allele):
 
 
 @dataclasses.dataclass
-class PythonSLiMMutationModel(PythonMutationModel):
-    mutation_type: int = 0
+class PythonSLiMv6MutationModel(PythonMutationModel):
     next_id: int = 0
 
     def root_allele(self, rng):
@@ -2000,6 +2186,11 @@ def transition_allele(self, rng, current_allele):
         return out
 
 
+@dataclasses.dataclass
+class PythonSLiMMutationModel(PythonSLiMv6MutationModel):
+    mutation_type: int = 0
+
+
 @dataclasses.dataclass
 class PythonInfiniteAllelesMutationModel(PythonMutationModel):
     start_allele: int = 0
@@ -2302,6 +2493,7 @@ def test_bad_models(self):
     def test_returns_mutation_model_instances_without_copying(self):
         models = [
             msprime.SLiMMutationModel(0, 0),
+            msprime.SLiMv6MutationModel(0),
             msprime.InfiniteAlleles(),
             msprime.BinaryMutationModel(),
             msprime.JC69(),
@@ -2324,6 +2516,11 @@ def test_slim(self):
         assert m.slim_generation == 9
         assert str(m) == "Mutation model for SLiM mutations of type m1\n  next ID: 2\n"
 
+    def test_slim_v6(self):
+        m = msprime.SLiMv6MutationModel(next_id=2)
+        assert m.next_id == 2
+        assert str(m) == "Mutation model for SLiM mutations, v6+. Next ID: 2\n"
+
     def test_infinite_alleles(self):
         m = msprime.InfiniteAlleles(start_allele=1)
         assert (
diff --git a/tests/test_python_c.py b/tests/test_python_c.py
index 5f12e5e7e..105b35c09 100644
--- a/tests/test_python_c.py
+++ b/tests/test_python_c.py
@@ -2808,6 +2808,30 @@ def test_generation(self):
         assert model.slim_generation == 1
 
 
+class TestSLiMv6MutationModel:
+    """
+    Tests for the slim v6 mutation model class.
+    """
+
+    def test_constructor_errors(self):
+        for bad_id in ["sdr", 0.222, None]:
+            with pytest.raises(TypeError):
+                _msprime.SLiMv6MutationModel(next_id=bad_id)
+
+        with pytest.raises(_msprime.LibraryError):
+            _msprime.SLiMv6MutationModel(next_id=-1)
+
+    def test_uninitialised(self):
+        model = _msprime.SLiMv6MutationModel.__new__(_msprime.SLiMv6MutationModel)
+        with pytest.raises(SystemError):
+            _ = model.next_id
+
+    def test_next_id(self):
+        for next_id in [0, 10, 2**63 - 1]:
+            model = _msprime.SLiMv6MutationModel(next_id=next_id)
+            assert model.next_id == next_id
+
+
 class TestInfiniteAllelesMutationModel:
     """
     Tests for the infinite alleles mutation model class.
@@ -2984,6 +3008,7 @@ def test_model(self):
                 _msprime.sim_mutations(tables, rng, rate_map=imap, model=bad_type)
         model_classes = [
             _msprime.SLiMMutationModel,
+            _msprime.SLiMv6MutationModel,
             _msprime.InfiniteAllelesMutationModel,
             _msprime.MatrixMutationModel,
         ]
@@ -3002,6 +3027,14 @@ def test_slim_model(self):
         _msprime.sim_mutations(tables, rng, imap, model)
         assert model.next_id == 5678
 
+    def test_slim_v6_model(self):
+        rng = _msprime.RandomGenerator(1)
+        imap = uniform_rate_map(1)
+        tables = _msprime.LightweightTableCollection(1.0)
+        model = _msprime.SLiMv6MutationModel(5678)
+        _msprime.sim_mutations(tables, rng, imap, model)
+        assert model.next_id == 5678
+
     def test_infinite_alleles_model(self):
         rng = _msprime.RandomGenerator(1)
         imap = uniform_rate_map(1, 1)