diff --git a/fwdpy11/_evolvets.py b/fwdpy11/_evolvets.py
index a94cc4829..9830d35d5 100644
--- a/fwdpy11/_evolvets.py
+++ b/fwdpy11/_evolvets.py
@@ -155,11 +155,19 @@ def evolvets(
                 assert isinstance(r, fwdpy11.mvDES), f"{type(r)}"
                 assert isinstance(r.des, fwdpy11.Sregion)  # type: ignore
                 nd = demographic_model.number_of_demes()
-                if nd > params.gvalue.ndemes:
-                    msg = f"maxmimum number of demes in model ({nd})"
-                    msg += " is not compatible with genetic value"
-                    msg += f" number of demes {params.gvalue.ndemes}"
-                    raise ValueError(msg)
+                if isinstance(params.gvalue, list):
+                    for g in params.gvalue:
+                        if nd > g.ndemes:
+                            msg = f"maxmimum number of demes in model ({nd})"
+                            msg += " is not compatible with genetic value"
+                            msg += f" number of demes {params.gvalue.ndemes}"
+                            raise ValueError(msg)
+                else:
+                    if nd > params.gvalue.ndemes:
+                        msg = f"maxmimum number of demes in model ({nd})"
+                        msg += " is not compatible with genetic value"
+                        msg += f" number of demes {params.gvalue.ndemes}"
+                        raise ValueError(msg)
                 if r.des.beg < 0:
                     raise ValueError(f"{r} has begin value < 0.0")
                 if r.des.end > pop.tables.genome_length:
diff --git a/tests/test_tree_sequences_different_des_in_different_demes.py b/tests/test_tree_sequences_different_des_in_different_demes.py
index 2767ff656..638ceb187 100644
--- a/tests/test_tree_sequences_different_des_in_different_demes.py
+++ b/tests/test_tree_sequences_different_des_in_different_demes.py
@@ -278,5 +278,84 @@ def test_invalid_model(mvDES, gvalue):
         fwdpy11.evolvets(rng, pop, params, 100)
 
 
+def test_two_demes_divergent_optima():
+    N = 500
+    rho = 0
+
+    yaml = f"""
+    time_units: generations
+    demes:
+      - name: alpha
+        epochs:
+          - start_size: {N}
+      - name: beta
+        epochs:
+          - start_size: {N}
+    migrations:
+      - demes: [alpha, beta]
+        rate: 1e-2
+    """
+
+    demography = fwdpy11.ForwardDemesGraph.from_demes(
+        yaml, burnin=N + 50, burnin_is_exact=True)
+
+    moving_optimum_deme_0 = fwdpy11.GaussianStabilizingSelection(
+        optima=[
+            fwdpy11.Optimum(when=0, optimum=0.0, VS=1.0),
+            fwdpy11.Optimum(when=N, optimum=1.0, VS=1.0),
+        ],
+        is_single_trait=True,
+    )
+    moving_optimum_deme_1 = fwdpy11.GaussianStabilizingSelection(
+        optima=[
+            fwdpy11.Optimum(when=0, optimum=0.0, VS=1.0),
+            fwdpy11.Optimum(when=N, optimum=-1.0, VS=1.0),
+            ],
+        is_single_trait=True,
+    )
+
+    # The covariance matrix for effect sizes.
+    # The marginal Gaussians will have mean zero and sd = 0.1
+    # The deviates will be highly positively correlated.
+    sd = 0.1
+    covariance_matrix = np.array([0.999] * 4).reshape((2, 2))
+    np.fill_diagonal(covariance_matrix, 1)
+    covariance_matrix *= np.power(sd, 2)
+
+    pdict = {
+        "nregions": [],
+        "sregions": [
+            # Multivariate Gaussian distribution of effect sizes
+            fwdpy11.mvDES(
+                fwdpy11.MultivariateGaussianEffects(
+                    beg=0, end=1, weight=1, h=1, cov_matrix=covariance_matrix
+                ),
+                # Means of zero for each marginal Gaussian
+                np.zeros(2),
+            )
+        ],
+        "recregions": [fwdpy11.PoissonInterval(0, 1, rho / (4 * N))],
+        "rates": (0, 1e-3, None),
+        "demography": demography,
+        "simlen": 10 * N + 200,  # 200 gens past optimum shift
+        # Specify one gvalue object per deme
+        "gvalue": [
+            fwdpy11.Additive(
+                ndemes=2,  # Number of demes
+                scaling=2,  # 0, 0+sh, 0+2s for AA, Aa, and aa, respectively
+                # Mapping of trait (genetic) value to fitness
+                gvalue_to_fitness=moving_optimum_deme_0,
+            ),
+            fwdpy11.Additive(ndemes=2, scaling=2,
+                             gvalue_to_fitness=moving_optimum_deme_1),
+        ],
+    }
+
+    params = fwdpy11.ModelParams(**pdict)
+    pop = fwdpy11.DiploidPopulation([N, N], 1.0)
+    rng = fwdpy11.GSLrng(1010)
+    fwdpy11.evolvets(rng, pop, params, 100)
+
+
 if __name__ == "__main__":
     unittest.main()