diff --git a/pyslim/methods.py b/pyslim/methods.py
index cf84634..674c22f 100644
--- a/pyslim/methods.py
+++ b/pyslim/methods.py
@@ -450,28 +450,34 @@ def add_mutation_metadata_tables(tables, mutation_type=0, remove_unused=False):
         for mut in tables.mutations
         for j in mut.derived_state.split(",")
     ]
-    mut_ids.sort()
-    mut_ids = np.array(mut_ids, dtype="int")  # floors times
-    # remove duplicate IDs, keeping the last (most recent)
-    keep = np.full(len(mut_ids), True, dtype="bool")
-    keep[np.where(np.diff(mut_ids[:, 0]) == 0)[0]] = False
-    mut_ids = mut_ids[keep, :]
-    mut_ids[:, 1] = slim_time(
-        tables, mut_ids[:, 1], stage="late", ts_metadata=ts_metadata
-    )
-    # this assumes mutations were added in late(), which is what SLiM does
-    ts_metadata["SLiM_mutation_list"].extend(
-        [
-            default_slim_metadata(
-                "mutation_list_entry",
-                mutation_id=int(j),
-                mutation_type=mutation_type,
-                slim_time=int(t),
-            )
-            for j, t in mut_ids
-            if j not in existing_muts
-        ]
-    )
+
+    # Only do things if there are mutations
+    # avoids an edge case where numpy indexing is out of bounds
+    # because empty lists create numpy arrays with only one dimension
+    if len(mut_ids) > 0:
+        mut_ids.sort()
+        mut_ids = np.array(mut_ids, dtype="int")  # floors times
+        # remove duplicate IDs, keeping the last (most recent)
+        keep = np.full(len(mut_ids), True, dtype="bool")
+        keep[np.where(np.diff(mut_ids[:, 0]) == 0)[0]] = False
+        mut_ids = mut_ids[keep, :]
+        mut_ids[:, 1] = slim_time(
+            tables, mut_ids[:, 1], stage="late", ts_metadata=ts_metadata
+        )
+        # this assumes mutations were added in late(), which is what SLiM does
+        ts_metadata["SLiM_mutation_list"].extend(
+            [
+                default_slim_metadata(
+                    "mutation_list_entry",
+                    num_traits=len(ts_metadata["SLiM"]["traits"]),
+                    mutation_id=int(j),
+                    mutation_type=mutation_type,
+                    slim_time=int(t),
+                )
+                for j, t in mut_ids
+                if j not in existing_muts
+            ]
+        )
     if remove_unused and len(mut_ids) < len(ts_metadata["SLiM_mutation_list"]):
         ts_metadata["SLiM_mutation_list"] = [
             x for x in ts_metadata["SLiM_mutation_list"] if x["mutation_id"] in mut_ids
@@ -486,7 +492,7 @@ def convert_alleles(ts):
     have "" (the empty string) for the ancestral state at each site; this method
     will replace this with the corresponding nucleotide from the reference sequence.
     For mutations, SLiM records the 'derived state' as a SLiM mutation ID; this
-    method will this with the nucleotide from the mutation's metadata.
+    method will replace this with the nucleotide from the mutation's metadata.
 
     This operation is not reversible: since SLiM mutation IDs are lost, the tree
     sequence will not be able to be read back into SLiM.