diff --git a/docs/metadata.md b/docs/metadata.md
index 9f6ff4b4..8e66ad6c 100644
--- a/docs/metadata.md
+++ b/docs/metadata.md
@@ -67,10 +67,10 @@ and `ts.metadata["SLiM"]` contains information about the simulation:
     * `index`: the index of the trait in SLiM
     * `name`: the name in SLiM for the trait
     * `type`: `"additive"`, `"multiplicative"`, or `"logistic"`
-    * `baselineOffsetFromUser`: a value added to all traits
-    * `baselineOffsetFromSubstitutions`: the total effect of all substitutions on the traits
+    * `baselineOffsetF`, `baselineOffsetM`, `baselineOffsetH`: values added to all individuals, depending on their sex
+    * `substitutionOffsetF`, `substitutionOffsetM`, `substitutionOffsetH`: total effect of all substitutions on traits, by sex,
       at the time the tree sequence was saved
-    * `baselineAccumulation`: whether the effect of substitutions accumulate in that value
+    * `substitutionAccumulation`: whether the effect of substitutions accumulate in that value
     * `directFitnessEffect`: whether the trait has a direct effect on fitness
     * `individualOffsetMean`, `individualOffsetSD`: parameters governing the distribution of individual-level offsets
       (i.e., "environment" effects)
@@ -144,17 +144,14 @@ and refer to that object, since otherwise you can incur runtime penalties
 for decoding and copying the metadata every time you call `ts.metadata`.
 This can be substantial, given the amount of mutation information
 in top-level metadata.
-For instance, to subtract baseline offsets from individual's trait values,
+For instance, to the top-level "tick" value from each node's time,
 we might do:
 ```{code-cell}
-md = ts.metadata
-traits = md["SLiM"]["traits"]
-values = [
-    [x['phenotype'] - y["baselineOffsetFromUser"] for x, y in zip(ind.metadata['per_trait'], traits)]
-    for ind in ts.individuals()
-]
+ts_metadata = ts.metadata
+t = ts_metadata["SLiM"]["tick"]
+adjusted_ticks = [n.time - t for n in ts.nodes()]
 ```
-If we instead inserted ``ts.metadata["SLiM"]["traits"]`` directly into the loop,
+If we instead inserted ``ts.metadata["SLiM"]["tick"]`` directly into the loop,
 this would become infeasibly slow.
 
 In some more detail:
diff --git a/docs/phenotypes.md b/docs/phenotypes.md
index 4685ad9b..0dcc3ee9 100644
--- a/docs/phenotypes.md
+++ b/docs/phenotypes.md
@@ -226,17 +226,23 @@ The final component is offsets.
 To match SLiM we need to add in both the individual's offset
 (which is stored in their metadata)
 and the global ("baseline") offset.
-Within SLiM, the `baselineOffset` property is a sum of the value provided by the user
-on initialization of the trait, and the accumulated value of any substitutions
+Within SLiM, the `baselineOffsetX` properties
+(where `X` is `M`, `F`, or `H`) values provided by the user
+on initialization of the trait;
+within SLiM this is added to the `substitutionOffsetX` property
 (if various options do not modify this behavior; see the SLiM manual).
 Since the tree sequence does not distinguish substitutions from other mutations,
-SLiM stores these two components separately, as `baselineOffsetFromUser`
-and `baselineOffsetFromSubstitutions`. So, we just need to add in `baselineOffsetFromUser`.
-
+it is important that SLiM stores these two components separately,
+and here we do nothing with the `substitutionOffset`.
 
 ```{code-cell}
 def additive_offset(ts_metadata, ind):
-    out = np.array([x['baselineOffsetFromUser'] for x in ts_metadata['SLiM']['traits']])
+    k = {
+        pyslim.INDIVIDUAL_TYPE_HERMAPHRODITE : "baselineOffsetH",
+        pyslim.INDIVIDUAL_TYPE_FEMALE : "baselineOffsetF",
+        pyslim.INDIVIDUAL_TYPE_MALE : "baselineOffsetM",
+    }[ind.metadata["sex"]]
+    out = np.array([x[k] for x in ts_metadata['SLiM']['traits']])
     out += [x['offset'] for x in ind.metadata['per_trait']]
     return out
 
@@ -254,7 +260,7 @@ print(f"SLiM: {np.array([x['phenotype'] for x in ind.metadata['per_trait']])}")
 ```
 
 That looks good! To compare we need to account for floating-point error
-(actually the numbers stored by SLiM and computed by us may differ by $10^{-10}$).
+(actually the numbers stored by SLiM and computed by us may differ by `1e-10` or so).
 
 
 ```{code-cell}
@@ -370,7 +376,12 @@ Offsets:
 
 ```{code-cell}
 def multiplicative_offset(ts_metadata, ind):
-    out = np.array([x['baselineOffsetFromUser'] for x in ts_metadata['SLiM']['traits']])
+    k = {
+        pyslim.INDIVIDUAL_TYPE_HERMAPHRODITE : "baselineOffsetH",
+        pyslim.INDIVIDUAL_TYPE_FEMALE : "baselineOffsetF",
+        pyslim.INDIVIDUAL_TYPE_MALE : "baselineOffsetM",
+    }[ind.metadata["sex"]]
+    out = np.array([x[k] for x in ts_metadata['SLiM']['traits']])
     out *= [x['offset'] for x in ind.metadata['per_trait']]
     return out
 
@@ -408,3 +419,72 @@ print(f"SLiM: {np.array([x['phenotype'] for x in ind.metadata['per_trait']])}")
 
 More elegant code would pull the trait types out of top-level metadata
 and use additive or multiplicative effects accordingly, etcetera.
+
+(sec_phenotypes_technical_details)=
+
+## Technical details
+
+You're better off using SLiM to calculate phenotypes
+than doing it yourself in python.
+There's lots of corner cases, and
+phenotypes can even be modified by the SLiM script directly.
+But, there's some situations where it's important to know about those corner cases.
+
+First: phenotypes are determined from the various contributions by
+either a sum ("additive" or "logistic") or a product ("multiplicative").
+Logistic traits are then transformed.
+The contributions come from the global `baselineOffsetX` values
+(which one depends on the sex of the individual),
+from individual offsets, and cumulative mutation effects.
+
+The contribution of 0, 1, or 2 copies of a SLiM mutation to a diploid
+is either 0, 2hs, or 2s (additive) or 1, 1+hs, 1+s (multiplicative),
+where h is the dominance coefficient and s is the effect size.
+For a hemizygous individual, the contributions are the same (for 0 and 1 copies),
+but using the hemizygous dominance coefficient for h.
+For an individual without the chromosome at all, there is no effect (obviously).
+
+Usually, that's all we need to know.
+However, there's some more complications that affect the default trait
+(i.e., the trait you get - called `simT` - if you don't explicitly declare any traits),
+or if you have set up a trait without baseline accumulation
+and with mutations that convert to substitutions.
+This can also be important to understand if you remove fixed mutations in python,
+and then want to read the file back into SLiM:
+to keep phenotypes the same, you have to include the effects of any removed mutations
+that count towards the trait in the `baselineOffsetX` value.
+
+Which mutations count?
+Within SLiM, when a Mutation fixes, it *might* convert into a Substitution.
+Whether this occurs is determined by
+the mutation type's `convertToSubstitution` property,
+which is True by default in WF models and False by default in nonWF models.
+Substitutions are no longer used by SLiM for calculating individual traits,
+and so to retain the effects of fixed mutations,
+sometimes scripts will set `convertToSubstitution=F` if it was not already set.
+Another mechanism to retain these effects is the Trait property `substitutionAccumulation`,
+which defaults to True.
+If a trait has `substitutionAccumulation=T`, then the effects of any Substitutions
+are accumulated in the "substitution offset", and thus contribute to every individual's
+phenotype. The net effect of this is that the contributions of any mutations that fix
+will be included in the phenotype values of all individuals
+as long as `convertToSubstitution=F` for their mutation type,
+or `substitutionAccumulation=T` for the trait, or both.
+
+However, the tree sequence does not record whether a given mutation was converted to
+a Substitution or not, and also does not record anything about mutation types,
+including their `convertToSubstitution` property.
+(This is because genome structure is set up before .trees files can be loaded.)
+So, when we find a fixed mutation in the tree sequence, we don't know *a priori*
+if this is a Substitution or not.
+If `substitutionAccumulation=T` for the trait
+(information that *is* in the top-level metadata),
+then it doesn't matter, mutations will contribute either way.
+The potentially missing information is in the script: are these mutation types
+converting to substitutions or not?
+
+Finally, a note on the `substitutionOffsetX` properties (where `X` is `M`, `F`, or `H`).
+This is where SLiM stores the effect of any substitutions on the traits of
+males, females, and heterozygotes, respectively.
+These are distinct because of the differing effects of fixed effects on individuals
+with 0, 1, or two copies of the chromosome in question.
diff --git a/pyslim/slim_metadata.py b/pyslim/slim_metadata.py
index 22b7fc87..dbebed8a 100644
--- a/pyslim/slim_metadata.py
+++ b/pyslim/slim_metadata.py
@@ -184,17 +184,17 @@ def is_vacant_num_bytes(num_chromosomes):
                             "description": "The traits defined for this tree sequence; each mutation and individual will have per-trait metadata.",
                             "items": {
                                 "properties": {
-                                    "baselineAccumulation": {
-                                        "description": "Whether the baseline offset includes accumulated effects from fixed (substituted) mutations.",
-                                        "type": "boolean",
+                                    "baselineOffsetH": {
+                                        "type": "number",
+                                        "description": "The baseline offset of the trait, for hermaphrodites.",
                                     },
-                                    "baselineOffsetFromUser": {
+                                    "baselineOffsetM": {
                                         "type": "number",
-                                        "description": "The from-user component of the baseline offset of the trait.",
+                                        "description": "The baseline offset of the trait, for males.",
                                     },
-                                    "baselineOffsetFromSubstitutions": {
+                                    "baselineOffsetF": {
                                         "type": "number",
-                                        "description": "The from-substitutions component of the baseline offset of the trait.",
+                                        "description": "The baseline offset of the trait, for females.",
                                     },
                                     "directFitnessEffect": {
                                         "description": "Whether the trait's effects are used directly as fitness effects.",
@@ -216,6 +216,22 @@ def is_vacant_num_bytes(num_chromosomes):
                                         "description": "The string name for the trait.",
                                         "type": "string",
                                     },
+                                    "substitutionAccumulation": {
+                                        "type": "boolean",
+                                        "description": "Whether the substitution offset accumulates effects from fixed (substituted) mutations.",
+                                    },
+                                    "substitutionOffsetH": {
+                                        "type": "number",
+                                        "description": "The substitution offset of the trait, for hermaphrodites.",
+                                    },
+                                    "substitutionOffsetM": {
+                                        "type": "number",
+                                        "description": "The substitution offset of the trait, for males.",
+                                    },
+                                    "substitutionOffsetF": {
+                                        "type": "number",
+                                        "description": "The substitution offset of the trait, for females.",
+                                    },
                                     "type": {
                                         "description": "The type of the trait; this must be 'additive', 'multiplicative', or 'logistic'.",
                                         "enum": [
diff --git a/tests/recipe_specs.py b/tests/recipe_specs.py
index fbe39a2d..8bd8e31e 100644
--- a/tests/recipe_specs.py
+++ b/tests/recipe_specs.py
@@ -27,10 +27,23 @@
         "pedigree": True,
         "record_mutations": True,
         "minimal": True,
+        "traits": True,  # just the default
+    },
+    "recipe_nonWF_X.slim": {
+        "nonWF": True,
+        "pedigree": True,
+        "X": True,
+        "minimal": True,
+        "traits": True,
+    },
+    "recipe_nonWF_Y.slim": {"nonWF": True, "pedigree": True, "Y": True, "traits": True},
+    "recipe_nonWF_H.slim": {
+        "nonWF": True,
+        "pedigree": True,
+        "H": True,
+        "minimal": True,
+        "traits": True,
     },
-    "recipe_nonWF_X.slim": {"nonWF": True, "pedigree": True, "X": True, "minimal": True},
-    "recipe_nonWF_Y.slim": {"nonWF": True, "pedigree": True, "Y": True},
-    "recipe_nonWF_H.slim": {"nonWF": True, "pedigree": True, "H": True, "minimal": True},
     "recipe_WF_X.slim": {
         "WF": True,
         "pedigree": True,
@@ -219,11 +232,24 @@
         "multichrom": True,
         "H-": True,
     },
-    "recipe_with_traits.slim": {
+    "recipe_with_traits_single_chrom.slim": {
         "WF": True,
         "traits": True,
+        "begun_late": True,
+    },
+    "recipe_with_traits_simple.slim": {
+        "WF": True,
+        "traits": True,
+        "begun_late": True,
         "multichrom": True,
+        "X": True,
+        "H": True,
+    },
+    "recipe_with_traits.slim": {
+        "WF": True,
+        "traits": True,
         "begun_late": True,
+        "multichrom": True,
         "X": True,
         "Y": True,
         "H": True,
diff --git a/tests/test_metadata.py b/tests/test_metadata.py
index bbd0288b..a299177d 100644
--- a/tests/test_metadata.py
+++ b/tests/test_metadata.py
@@ -2,10 +2,13 @@
 Test cases for the metadata reading/writing of pyslim.
 """
 
+import dataclasses
 import sys
+from collections import Counter
 
 import numpy as np
 import pytest
+import scipy.sparse as sparse
 import tskit
 
 import pyslim
@@ -424,3 +427,250 @@ def test_chromosome_types(self, recipe):
         for chrom in chroms.keys():
             assert chrom in chrom_list_d
             assert chrom_list_d[chrom] == chrom_info[chrom]
+
+
+class TraitCalculator:
+    """
+    Calculating phenotypes by hand from metadata takes a good bit of
+    internal state, especially if it's from a multi-chromosome simulation,
+    so here's a class that does it.  See below for how to use this.
+    """
+
+    def __init__(self):
+        self.ts_metadata = None
+        self.slim_phenotypes = None
+        self.phenotypes = None
+        self.offsets = None
+
+    def add_ts(self, ts):
+        ts_metadata = ts.metadata
+        if self.ts_metadata is not None:
+            assert ts_metadata["SLiM"]["traits"] == self.ts_metadata["SLiM"]["traits"]
+        self.ts_metadata = ts_metadata
+        self.ts = ts
+        chrom_type = self.ts_metadata["SLiM"]["this_chromosome"]["type"]
+        self.haploid_chromosome = chrom_type in ("H", "HF", "HM")
+        self.mut_metadata = pyslim.mutation_metadata(ts)
+        self.nodes_vacant = pyslim.nodes_vacant(ts)
+        self.sexes = self.ts.tables.individuals.metadata_vector("sex")
+        self.num_traits = len(self.ts_metadata["SLiM"]["traits"])
+        self.types = [m["type"] for m in self.ts_metadata["SLiM"]["traits"]]
+        self.accumulates = [
+            m["substitutionAccumulation"] for m in self.ts_metadata["SLiM"]["traits"]
+        ]
+        # to properly do traits in the case where substitutionAccumulation=F
+        # and there are mutation types with convertToSubstitution=T then we
+        # need to store this in metadata: (so, this code won't work if this is
+        # not present for such traits, which are probably only the default trait!)
+        if (
+            "user_metadata" in ts_metadata["SLiM"]
+            and "mutation_types" in ts_metadata["SLiM"]["user_metadata"]
+        ):
+            self.mutation_types = {
+                k: v[0]["convertToSubstitution"][0]
+                for k, v in ts_metadata["SLiM"]["user_metadata"]["mutation_types"][
+                    0
+                ].items()
+            }
+        else:
+            # guess T if there just a single trait without substitution accumulation
+            # whose name is simT and this is a WF model
+            # since that describes "default trait"
+            convert = (
+                self.ts_metadata["SLiM"]["model_type"] == "WF"
+                and len(self.ts_metadata["SLiM"]["traits"]) == 1
+                and self.ts_metadata["SLiM"]["traits"][0]["name"] == "simT"
+                and self.ts_metadata["SLiM"]["traits"][0]["substitutionAccumulation"]
+                == False
+            )
+            mtypes = set(md["mutation_type"] for md in self.mut_metadata.values())
+            self.mutation_types = {f"m{k}": convert for k in mtypes}
+        sp = np.column_stack(
+            [
+                self.ts.tables.individuals.metadata_vector(["per_trait", j, "phenotype"])
+                for j in range(self.num_traits)
+            ]
+        )
+        if self.slim_phenotypes is None:
+            self.slim_phenotypes = sp
+        assert np.all(self.slim_phenotypes == sp)
+        if self.phenotypes is None:
+            # initialize phenotypes where we'll accumulate changes
+            self.phenotypes = np.full(
+                self.slim_phenotypes.shape,
+                [1.0 if t == "multiplicative" else 0.0 for t in self.types],
+                dtype="float",
+            )
+        offsets = np.column_stack(
+            [
+                self.ts.tables.individuals.metadata_vector(["per_trait", j, "offset"])
+                for j in range(self.num_traits)
+            ]
+        )
+        if self.offsets is None:
+            self.offsets = offsets
+            # only do offsets if this is the first ts!
+            self.do_offsets()
+        else:
+            assert np.all(offsets == self.offsets)
+        self.get_frequencies()
+        self.do_effects()
+
+    def transform(self):
+        for k, t in enumerate(self.types):
+            if t == "logistic":
+                self.phenotypes[:, k] = 1 / (1 + np.exp(-self.phenotypes[:, k]))
+
+    def is_converted(self, a):
+        # is allele a from a convert-to-substitution mutation type?
+        # if we haven't recorded whether or not it's converted, assume T
+        md = self.mut_metadata[int(a)]
+        return self.mutation_types[f"m{md['mutation_type']}"]
+
+    def get_frequencies(self):
+        # Fixed mutations count unless they are converted to substitutions
+        # and substitution accumulation is off, so we need to look up both these things:
+        # (frequency, whether they're converted to substitutions) for each mutation.
+        self.frequencies = {}
+        for v in self.ts.variants(isolated_as_missing=False):
+            freqs = {}
+            for ds, f in v.counts().items():
+                for a in ds.split(","):
+                    if a != "":
+                        if a not in freqs:
+                            convert = self.is_converted(a)
+                            freqs[a] = [0, convert]
+                        freqs[a][0] += f
+            self.frequencies[v.site.id] = freqs
+
+    def do_offsets(self):
+        okeys = {
+            pyslim.INDIVIDUAL_TYPE_HERMAPHRODITE: "baselineOffsetH",
+            pyslim.INDIVIDUAL_TYPE_FEMALE: "baselineOffsetF",
+            pyslim.INDIVIDUAL_TYPE_MALE: "baselineOffsetM",
+        }
+        for k, t in enumerate(self.types):
+            bO = np.array(
+                [self.ts_metadata["SLiM"]["traits"][k][okeys[s]] for s in self.sexes]
+            )
+            if t == "multiplicative":
+                self.phenotypes[:, k] *= bO * self.offsets[:, k]
+            else:
+                self.phenotypes[:, k] += bO + self.offsets[:, k]
+
+    def do_effects(self):
+        for ind in self.ts.individuals():
+            ploidy = np.sum(~self.nodes_vacant[ind.nodes])
+            if ploidy == 0:
+                continue
+            else:
+                hemizygous = ploidy == 1
+                for k, t in enumerate(self.types):
+                    if t == "multiplicative":
+                        g = self.multiplicative_effect(ind, k, hemizygous)
+                        self.phenotypes[ind.id, k] *= g
+                    else:
+                        g = self.additive_effect(ind, k, hemizygous)
+                        self.phenotypes[ind.id, k] += g
+
+    def skip_mut(self, m, freqs, trait_id):
+        if m == "":
+            out = True
+        else:
+            f, convert = freqs[m]
+            out = f == 0 or (
+                f == self.ts.num_samples and convert and not self.accumulates[trait_id]
+            )
+        return out
+
+    def additive_effect(self, ind, trait_id, hemizygous):
+        out = 0.0
+        for v in self.ts.variants(samples=ind.nodes, isolated_as_missing=False):
+            freqs = self.frequencies[v.site.id]
+            a = ",".join([v.alleles[g] for g in v.genotypes])
+            muts = Counter(a.split(","))
+            for m in muts:
+                if self.skip_mut(m, freqs, trait_id):
+                    continue
+                md = self.mut_metadata[int(m)]["per_trait"][trait_id]
+                s = md["effect_size"]
+                if muts[m] == 2:
+                    h = 1
+                else:
+                    assert muts[m] == 1
+                    if self.haploid_chromosome:
+                        h = 1
+                    else:
+                        if hemizygous:
+                            h = md["hemizygous_dominance"]
+                        else:
+                            h = md["dominance"]
+                        if np.isnan(h):
+                            h = 1 / 2
+                out += 2 * h * s
+        return out
+
+    def multiplicative_effect(self, ind, trait_id, hemizygous):
+        out = 1.0
+        for v in self.ts.variants(samples=ind.nodes, isolated_as_missing=False):
+            freqs = self.frequencies[v.site.id]
+            a = ",".join([v.alleles[g] for g in v.genotypes])
+            muts = Counter(a.split(","))
+            for m in muts:
+                if self.skip_mut(m, freqs, trait_id):
+                    continue
+                assert muts[m] > 0 and muts[m] <= 2
+                md = self.mut_metadata[int(m)]["per_trait"][trait_id]
+                s = md["effect_size"]
+                if muts[m] == 2:
+                    h = 1
+                else:
+                    assert muts[m] == 1
+                    if self.haploid_chromosome:
+                        h = 1
+                    elif hemizygous:
+                        h = md["hemizygous_dominance"]
+                    else:
+                        h = md["dominance"]
+                    if np.isnan(h):
+                        # "independent dominance occurs when (1+hs)(1+hs) equals 1+s,
+                        # which occurs when h=(sqrt(1+s)−1)/s"
+                        h = (np.sqrt(1 + s) - 1) / s if s != 0 else 0
+                out *= 1 + h * s
+        return out
+
+
+class TestTraits(tests.PyslimTestCase):
+    @pytest.mark.parametrize("recipe", recipe_eq("traits"), indirect=True)
+    def test_traits_consistency(self, recipe):
+        trait_md = None
+        offsets = None
+        slim_phenotypes = None
+        tc = TraitCalculator()
+        for ts in recipe["ts"].values():
+            assert ts.num_mutations > 0
+            ts_metadata = ts.metadata
+            # top-level info about traits
+            if trait_md is None:
+                trait_md = ts_metadata["SLiM"]["traits"]
+            assert trait_md == ts_metadata["SLiM"]["traits"]
+            # now compute them
+            tc.add_ts(ts)
+        # phenotypes as computed by us
+        tc.transform()
+        print([ind.metadata["sex"] for ind in ts.individuals()])
+        for k in range(len(tc.types)):
+            print(k, tc.types[k])
+            print(
+                np.column_stack(
+                    [
+                        tc.phenotypes[:, k],
+                        tc.slim_phenotypes[:, k],
+                        tc.phenotypes[:, k] - tc.slim_phenotypes[:, k],
+                    ]
+                )
+            )
+        print(np.max(tc.phenotypes - tc.slim_phenotypes))
+        print(np.min(tc.phenotypes - tc.slim_phenotypes))
+        print(np.sum(tc.phenotypes - tc.slim_phenotypes))
+        assert np.allclose(tc.phenotypes, tc.slim_phenotypes)
diff --git a/tests/test_recipes/recipe_WF_H.slim b/tests/test_recipes/recipe_WF_H.slim
index 99065111..93a8aa04 100644
--- a/tests/test_recipes/recipe_WF_H.slim
+++ b/tests/test_recipes/recipe_WF_H.slim
@@ -29,6 +29,7 @@ initialize()
 }
 
 10 late() {
+    sim.demandPhenotype(sim.subpopulations);
     sim.treeSeqOutput(TREES_FILE);
     catn("Done.");
     sim.simulationFinished();
diff --git a/tests/test_recipes/recipe_WF_HF.slim b/tests/test_recipes/recipe_WF_HF.slim
index 7e71e21f..13f712b1 100644
--- a/tests/test_recipes/recipe_WF_HF.slim
+++ b/tests/test_recipes/recipe_WF_HF.slim
@@ -30,6 +30,7 @@ initialize()
 }
 
 10 late() {
+    sim.demandPhenotype(sim.subpopulations);
     sim.treeSeqOutput(TREES_FILE);
     catn("Done.");
     sim.simulationFinished();
diff --git a/tests/test_recipes/recipe_WF_HM.slim b/tests/test_recipes/recipe_WF_HM.slim
index 12ca7be4..5604f8a2 100644
--- a/tests/test_recipes/recipe_WF_HM.slim
+++ b/tests/test_recipes/recipe_WF_HM.slim
@@ -30,6 +30,7 @@ initialize()
 }
 
 10 late() {
+    sim.demandPhenotype(sim.subpopulations);
     sim.treeSeqOutput(TREES_FILE);
     catn("Done.");
     sim.simulationFinished();
diff --git a/tests/test_recipes/recipe_WF_W.slim b/tests/test_recipes/recipe_WF_W.slim
index ec7d9fb5..3b3c4078 100644
--- a/tests/test_recipes/recipe_WF_W.slim
+++ b/tests/test_recipes/recipe_WF_W.slim
@@ -30,6 +30,7 @@ initialize()
 }
 
 10 late() {
+    sim.demandPhenotype(sim.subpopulations);
     sim.treeSeqOutput(TREES_FILE);
     catn("Done.");
     sim.simulationFinished();
diff --git a/tests/test_recipes/recipe_WF_X.slim b/tests/test_recipes/recipe_WF_X.slim
index 5f40dba2..f12f0c9c 100644
--- a/tests/test_recipes/recipe_WF_X.slim
+++ b/tests/test_recipes/recipe_WF_X.slim
@@ -30,6 +30,7 @@ initialize()
 }
 
 10 late() {
+    sim.demandPhenotype(sim.subpopulations);
     sim.treeSeqOutput(TREES_FILE);
     catn("Done.");
     sim.simulationFinished();
diff --git a/tests/test_recipes/recipe_WF_Y.slim b/tests/test_recipes/recipe_WF_Y.slim
index 50172112..022cbc14 100644
--- a/tests/test_recipes/recipe_WF_Y.slim
+++ b/tests/test_recipes/recipe_WF_Y.slim
@@ -30,6 +30,7 @@ initialize()
 }
 
 10 late() {
+    sim.demandPhenotype(sim.subpopulations);
     sim.treeSeqOutput(TREES_FILE);
     catn("Done.");
     sim.simulationFinished();
diff --git a/tests/test_recipes/recipe_WF_Z.slim b/tests/test_recipes/recipe_WF_Z.slim
index c680f145..30c3d4fe 100644
--- a/tests/test_recipes/recipe_WF_Z.slim
+++ b/tests/test_recipes/recipe_WF_Z.slim
@@ -30,6 +30,7 @@ initialize()
 }
 
 10 late() {
+    sim.demandPhenotype(sim.subpopulations);
     sim.treeSeqOutput(TREES_FILE);
     catn("Done.");
     sim.simulationFinished();
diff --git a/tests/test_recipes/recipe_nonWF.slim b/tests/test_recipes/recipe_nonWF.slim
index 7174a853..044b53c2 100644
--- a/tests/test_recipes/recipe_nonWF.slim
+++ b/tests/test_recipes/recipe_nonWF.slim
@@ -45,7 +45,7 @@ mutation() {
 	muts.setValue(asString(mut.id), m);
 	return T;
 }
-10 late() {
+20 late() {
 	subs = Dictionary();
 	for (mut in sim.substitutions) {
 		nuc = mut.mutationType.nucleotideBased ? mut.nucleotide else "N";
@@ -62,7 +62,7 @@ mutation() {
 }
 
 // OUTPUT/FINISH
-10 late() {
+20 late() {
     sim.treeSeqOutput(TREES_FILE, metadata=MD);
     catn("Done.");
     sim.simulationFinished();
diff --git a/tests/test_recipes/recipe_nonWF_H.slim b/tests/test_recipes/recipe_nonWF_H.slim
index 2baaa231..b1467d01 100644
--- a/tests/test_recipes/recipe_nonWF_H.slim
+++ b/tests/test_recipes/recipe_nonWF_H.slim
@@ -28,7 +28,7 @@ early() {
     p1.fitnessScaling = K / p1.individualCount;
 }
 
-10 late() {
+30 late() {
     sim.treeSeqOutput(TREES_FILE);
     catn("Done.");
     sim.simulationFinished();
diff --git a/tests/test_recipes/recipe_nonWF_X.slim b/tests/test_recipes/recipe_nonWF_X.slim
index 5250e81a..932f49bd 100644
--- a/tests/test_recipes/recipe_nonWF_X.slim
+++ b/tests/test_recipes/recipe_nonWF_X.slim
@@ -29,7 +29,7 @@ early() {
     p1.fitnessScaling = K / p1.individualCount;
 }
 
-10 late() {
+30 late() {
     sim.treeSeqOutput(TREES_FILE);
     catn("Done.");
     sim.simulationFinished();
diff --git a/tests/test_recipes/recipe_nonWF_Y.slim b/tests/test_recipes/recipe_nonWF_Y.slim
index 01e6e32f..4e0e1137 100644
--- a/tests/test_recipes/recipe_nonWF_Y.slim
+++ b/tests/test_recipes/recipe_nonWF_Y.slim
@@ -9,7 +9,7 @@ initialize()
     initializeSex();
     initializeChromosome(1, 100, "Y", "Y");
     initializeMutationRate(1e-2);
-    initializeMutationType("m1", 0.5, "f", -0.1);
+    initializeMutationType("m1", 0.5, "f", -0.01);
     initializeGenomicElementType("g1", m1, 1.0);
     initializeGenomicElement(g1, 0, 99);
     initializeRecombinationRate(1e-2);
@@ -29,7 +29,7 @@ early() {
     p1.fitnessScaling = K / p1.individualCount;
 }
 
-10 late() {
+30 late() {
     sim.treeSeqOutput(TREES_FILE);
     catn("Done.");
     sim.simulationFinished();
diff --git a/tests/test_recipes/recipe_nonWF_selfing.slim b/tests/test_recipes/recipe_nonWF_selfing.slim
deleted file mode 100644
index 379b8bdd..00000000
--- a/tests/test_recipes/recipe_nonWF_selfing.slim
+++ /dev/null
@@ -1,84 +0,0 @@
-initialize()
-{
-    setSeed(23);
-    if (!exists("TREES_FILE")) defineGlobal("TREES_FILE", "out.trees"); 
-    if (!exists("PEDIGREE_FILE")) defineGlobal("PEDIGREE_FILE", "out.pedigree"); 
-    initializeSLiMModelType("nonWF");
-    initializeSLiMOptions(keepPedigrees=T);
-    initializeTreeSeq();
-    initializeMutationRate(1e-2);
-    initializeMutationType("m1", 0.5, "f", -0.1);
-    initializeGenomicElementType("g1", m1, 1.0);
-    initializeGenomicElement(g1, 0, 99);
-    initializeRecombinationRate(1e-2);
-    defineConstant("K", 10);
-}
-
-reproduction() {
-    if (runif(1) < 0.25) {
-        mate = subpop.sampleIndividuals(1);
-    } else {
-        mate = individual;
-    }
-    for (_ in seqLen(rpois(1, 2))) {
-      subpop.addCrossed(individual, mate);
-    }
-}
-
-1 early() {
-    sim.addSubpop("p1", 10);
-}
-
-early() {
-    p1.fitnessScaling = K / p1.individualCount;
-}
-
-10 late() {
-    sim.treeSeqOutput(TREES_FILE);
-    catn("Done.");
-    sim.simulationFinished();
-}
-
-// PEDIGREE OUTPUT
-1 first() {
-    writeFile(PEDIGREE_FILE,
-       paste(c("generation", "stage", "individual", "age", "parent1", "parent2"),
-             sep="\t"));
-}
-
-1: first() {
-    for (pop in sim.subpopulations) {
-        for (ind in pop.individuals) {
-            age = community.modelType == "nonWF" ? ind.age else -1;
-            writeFile(PEDIGREE_FILE,
-             paste(c(sim.cycle, community.cycleStage, ind.pedigreeID, age, ind.pedigreeParentIDs),
-                   sep='\t'),
-             append=T);
-        }
-    }
-}
-
-1: early() {
-    for (pop in sim.subpopulations) {
-        for (ind in pop.individuals) {
-            age = community.modelType == "nonWF" ? ind.age else -1;
-            writeFile(PEDIGREE_FILE,
-             paste(c(sim.cycle, community.cycleStage, ind.pedigreeID, age, ind.pedigreeParentIDs),
-                   sep='\t'),
-             append=T);
-        }
-    }
-}
-
-1: late() {
-    for (pop in sim.subpopulations) {
-        for (ind in pop.individuals) {
-            age = community.modelType == "nonWF" ? ind.age else -1;
-            writeFile(PEDIGREE_FILE,
-             paste(c(sim.cycle, community.cycleStage, ind.pedigreeID, age, ind.pedigreeParentIDs),
-                   sep='\t'),
-             append=T);
-        }
-    }
-}
-
diff --git a/tests/test_recipes/recipe_with_traits.slim b/tests/test_recipes/recipe_with_traits.slim
index 72602d26..63b661b9 100644
--- a/tests/test_recipes/recipe_with_traits.slim
+++ b/tests/test_recipes/recipe_with_traits.slim
@@ -12,16 +12,15 @@ initialize() {
 	
 	initializeSex();
 	
-	popgen1T = initializeTrait("popgen1T", "m", 1.0, 0.0, 0.01, directFitnessEffect=T);
-	popgen2T = initializeTrait("popgen2T", "m", 1.0, 0.0, 0.01, directFitnessEffect=T);
+	popgen1T = initializeTrait("popgen1T", "m", c(1.0, 1.1), 0.0, 0.01, directFitnessEffect=T);
+	popgen2T = initializeTrait("popgen2T", "m", c(0.9, 1.0), 0.0, 0.01, directFitnessEffect=T);
 	n1T = initializeTrait("n1T", "m", directFitnessEffect=T);
 	n2T = initializeTrait("n2T", "m", directFitnessEffect=F);
 	
 	quant1T = initializeTrait("quant1T", "a", I1, 0.0, 0.01, directFitnessEffect=F);
 	quant2T = initializeTrait("quant2T", "a", I2, 0.0, 0.01, directFitnessEffect=F);
-	n3T = initializeTrait("n3T", "a", directFitnessEffect=F, baselineAccumulation=F);
+	n3T = initializeTrait("n3T", "a", directFitnessEffect=F);
 	
-
 	logistic1T = initializeTrait("logistic1T", "l", 0.0, 0.01, 0.01, directFitnessEffect=T);	
 	initializeMutationType("m1", 0.4, "f", 0.0);
 	initializeMutationType("m2", 0.4, "e", 0.05);
@@ -45,20 +44,20 @@ initialize() {
 	
 	ids = 1:5;
 	symbols = c(1, 2, "X", "Y", "MT");
-	lengths = rdunif(5, 1e7, 2e7);
+	lengths = rdunif(5, 1e3, 2e3);
 	types = c("A", "A", "X", "Y", "H");
 	names = c("A1", "A2", "X", "Y", "MT");
 	
 	for (id in ids, symbol in symbols, length in lengths, type in types, name in names)
 	{
 		initializeChromosome(id, length, type, symbol, name);
-		initializeMutationRate(1e-7);
-		initializeRecombinationRate(1e-8);
+		initializeMutationRate(1e-3);
+		initializeRecombinationRate(1e-4);
 		
 		if (id == 1)
 			initializeGenomicElement(g1); // autosome 1 is pure-neutral, using only m1
 		else
-			initializeGenomicElement(g2); // autosome 2 is a mix, using m1 / m2 / m3
+			initializeGenomicElement(g2); // others are a mix, using m1 / m2 / m3
 	}
 }
 
@@ -95,6 +94,7 @@ mutation(m3) {
 
 
 10 late() {
+    sim.demandPhenotype(sim.subpopulations);
     sim.treeSeqOutput(TREES_FILE);
     catn("Done.");
     sim.simulationFinished();
diff --git a/tests/test_recipes/recipe_with_traits_simple.slim b/tests/test_recipes/recipe_with_traits_simple.slim
new file mode 100644
index 00000000..22ca6204
--- /dev/null
+++ b/tests/test_recipes/recipe_with_traits_simple.slim
@@ -0,0 +1,45 @@
+initialize() {
+   setSeed(23);
+   if (!exists("TREES_FILE")) defineGlobal("TREES_FILE", "out.trees"); 
+   initializeSLiMOptions(keepPedigrees=T);
+   initializeTreeSeq(timeUnit="generations");
+
+   initializeSex();
+
+   popgen1T = initializeTrait("popgen1T", "l", c(-1.0, 1.1), 0.0, 0.01, directFitnessEffect=T);
+   initializeMutationType("m1", 0.4, "e", 0.05);
+   
+   initializeGenomicElementType("g1", m1, 1.0);
+   
+   symbols = c('1', '2', 'X', 'H');
+   types = c("A", 'A', 'X', 'H');
+   names = c("A1", 'A2', 'X', 'H');
+   ids = seqAlong(names);
+   lengths = rdunif(size(ids), 1e3, 2e3);
+   
+   for (id in ids, symbol in symbols, length in lengths, type in types, name in names)
+   {
+      initializeChromosome(id, length, type, symbol, name);
+      initializeMutationRate(1e-3);
+      initializeRecombinationRate(1e-4);
+      
+      initializeGenomicElement(g1);
+   }
+}
+
+mutation(m1) {
+	// set random dominance effects
+	mut.popgen1TDominance = runif(1);
+	return T;
+}
+
+1 late() {
+   sim.addSubpop("p1", 20);
+}
+
+10 late() {
+    sim.demandPhenotype(sim.subpopulations);
+    sim.treeSeqOutput(TREES_FILE);
+    catn("Done.");
+    sim.simulationFinished();
+}
diff --git a/tests/test_recipes/recipe_with_traits_single_chrom.slim b/tests/test_recipes/recipe_with_traits_single_chrom.slim
new file mode 100644
index 00000000..a748aa7b
--- /dev/null
+++ b/tests/test_recipes/recipe_with_traits_single_chrom.slim
@@ -0,0 +1,87 @@
+initialize() {
+   setSeed(23);
+   if (!exists("TREES_FILE")) defineGlobal("TREES_FILE", "out.trees"); 
+   initializeSLiMOptions(keepPedigrees=T);
+   initializeTreeSeq(timeUnit="generations");
+	defineConstant("I1", 5.0);
+	defineConstant("I2", -5.0);
+	defineConstant("OPT1", 10.0);
+	defineConstant("OPT2", 10.0);
+	defineConstant("SD1", 2.0);
+	defineConstant("SD2", 2.0);
+	
+	initializeSex();
+	
+	popgen1T = initializeTrait("popgen1T", "m", 1.0, 0.0, 0.01, directFitnessEffect=T);
+	popgen2T = initializeTrait("popgen2T", "m", 1.0, 0.0, 0.01, directFitnessEffect=T);
+	n1T = initializeTrait("n1T", "m", directFitnessEffect=T);
+	n2T = initializeTrait("n2T", "m", directFitnessEffect=F);
+	
+	quant1T = initializeTrait("quant1T", "a", I1, 0.0, 0.01, directFitnessEffect=F);
+	quant2T = initializeTrait("quant2T", "a", I2, 0.0, 0.01, directFitnessEffect=F);
+	n3T = initializeTrait("n3T", "a", directFitnessEffect=F, substitutionAccumulation=F);
+	
+
+	logistic1T = initializeTrait("logistic1T", "l", 0.0, 0.01, 0.01, directFitnessEffect=T);	
+	initializeMutationType("m1", 0.4, "f", 0.0);
+	initializeMutationType("m2", 0.4, "e", 0.05);
+	m2.setEffectSizeDistributionForTrait(c(n1T, n2T), "f", 0.0);
+	m2.setEffectSizeDistributionForTrait(c(quant1T, quant2T), "n", 0.0, 0.1);
+	m2.setEffectSizeDistributionForTrait(c(logistic1T), "n", -0.05, 0.1);
+	
+	initializeMutationType("m3", 0.4, "g", -0.05, 1.0);
+	m3.setEffectSizeDistributionForTrait(c(n1T, n2T), "f", 0.0);
+	m3.setEffectSizeDistributionForTrait(c(quant1T, quant2T), "n", 0.0, 0.1);
+	m3.setEffectSizeDistributionForTrait(c(logistic1T), "n", -0.05, 0.1);
+	
+	c(m2,m3).setEffectSizeDistributionForTrait(n3T, "n", -5.0, 0.5);
+	
+	c(m2,m3).setDefaultDominanceForTrait(c(popgen2T, quant2T), NAN);
+	
+	c(m2,m3).logMutationData(T, trait=NULL, effectSize=T, dominance=T);
+	
+   initializeMutationRate(1e-3);
+	initializeGenomicElementType("g1", m1, 1.0);
+	initializeGenomicElementType("g2", 1:3, c(3, 1, 2));
+   initializeGenomicElement(g1, 0, 1000);
+   initializeRecombinationRate(1e-4);
+}
+
+mutation(m2) {
+	// set random dominance effects for the popgen1T and quant1T and logistic1TDominance traits
+	// other effects are generated as specified by the mutation type DES
+	mut.popgen1TDominance = runif(1);
+	mut.quant1TDominance = runif(1);
+	mut.logistic1TDominance = runif(1);
+	return T;
+}
+mutation(m3) {
+	// set random dominance effects for the popgen1T and quant1T and logistic1TDominance traits
+	// other effects are generated as specified by the mutation type DES
+	mut.popgen1TDominance = runif(1);
+	mut.quant1TDominance = runif(1);
+	mut.logistic1TDominance = runif(1);
+	return T;
+}
+
+1 late() {
+	sim.addSubpop("p1", 20);
+}
+
+1: late() {
+	inds = sim.subpopulations.individuals;
+   sim.demandPhenotype(NULL, c(sim.quant1T, sim.quant2T));
+	phenotypes_q1 = inds.quant1T;
+	phenotypes_q2 = inds.quant2T;
+	fitnessEffect_q1 = dnorm(phenotypes_q1, OPT1, SD1) / dnorm(0.0, 0.0, SD1);
+	fitnessEffect_q2 = dnorm(phenotypes_q2, OPT2, SD2) / dnorm(0.0, 0.0, SD2);
+	inds.fitnessScaling = fitnessEffect_q1 * fitnessEffect_q2;
+}
+
+
+10 late() {
+    sim.demandPhenotype(sim.subpopulations);
+    sim.treeSeqOutput(TREES_FILE);
+    catn("Done.");
+    sim.simulationFinished();
+}