diff --git a/cgp/genome.py b/cgp/genome.py
index bf4d3240..fde198e0 100644
--- a/cgp/genome.py
+++ b/cgp/genome.py
@@ -349,7 +349,7 @@ def _is_active_input_gene(self, gene_idx: int) -> bool:
         else:
             assert False  # should never be reached
 
-    def _select_gene_indices_for_mutation(self, mutation_rate, len_dna):
+    def _select_gene_indices_for_mutation(self, mutation_rate, len_dna, rng):
         """Selects the gene indices for mutations
 
         Parameters
@@ -358,6 +358,8 @@ def _select_gene_indices_for_mutation(self, mutation_rate, len_dna):
             Proportion of genes to be mutated, between 0 and 1.
         len_dna : int
             Length of the genome dna.
+        rng : numpy.random.RandomState
+            Random number generator instance to use for selecting the indices.
 
         Returns
         ----------
@@ -365,7 +367,7 @@ def _select_gene_indices_for_mutation(self, mutation_rate, len_dna):
             indices of the genes selected for mutation.
         """
 
-        selected_gene_indices = np.nonzero(np.random.rand(len_dna) < mutation_rate)[0]
+        selected_gene_indices = np.nonzero(rng.rand(len_dna) < mutation_rate)[0]
         return selected_gene_indices
 
     def mutate(self, mutation_rate: float, rng: np.random.RandomState):
@@ -389,7 +391,9 @@ def mutate(self, mutation_rate: float, rng: np.random.RandomState):
         dna = list(self._dna)
         only_silent_mutations = True
 
-        selected_gene_indices = self._select_gene_indices_for_mutation(mutation_rate, len(dna))
+        selected_gene_indices = self._select_gene_indices_for_mutation(
+            mutation_rate, len(dna), rng
+        )
 
         for (gene_idx, allele) in zip(selected_gene_indices, np.array(dna)[selected_gene_indices]):
 
diff --git a/test/test_genome.py b/test/test_genome.py
index 876e1bf6..8c65b155 100644
--- a/test/test_genome.py
+++ b/test/test_genome.py
@@ -449,11 +449,12 @@ def test_only_silent_mutations(genome_params, mutation_rate, rng_seed):
         3 * length_per_region
     )  # operator gene of 2nd internal unit, a silent unit in this dna
 
-    def select_gene_indices_silent(mutation_rate, dna):
+    def select_gene_indices_silent(mutation_rate, dna, rng):
         selected_gene_indices = [gene_to_be_mutated_non_active]
         return selected_gene_indices
 
-    genome._select_gene_indices_for_mutation = select_gene_indices_silent  # monkey patch the selection of indices to selecting a silent gene
+    genome._select_gene_indices_for_mutation = \
+        select_gene_indices_silent  # monkey patch the selection of indices to selecting a silent gene
     genome.dna = dna_fixed
     only_silent_mutations = genome.mutate(mutation_rate, rng)
     assert only_silent_mutations is True
@@ -462,7 +463,7 @@ def select_gene_indices_silent(mutation_rate, dna):
         active_regions[-1] * length_per_region
     )  # function gene of the 1st active internal gene, should always be mutable
 
-    def select_gene_indices_non_silent(mutation_rate, dna):
+    def select_gene_indices_non_silent(mutation_rate, dna, rng):
         selected_gene_indices = [gene_to_be_mutated_active]
         return selected_gene_indices
 
@@ -487,7 +488,8 @@ def test_permissible_values_hidden(rng_seed):
     gene_idx = 6  # function gene of first hidden region
     allele = 0  #
     region_idx = None  # can be any number, since not touched for function gene
-    permissible_function_gene_values = genome._determine_alternative_permissible_values_hidden_gene(
+    permissible_function_gene_values = \
+        genome._determine_alternative_permissible_values_hidden_gene(
         gene_idx, allele, region_idx
     )
     assert permissible_function_gene_values == [
@@ -501,7 +503,8 @@ def test_permissible_values_hidden(rng_seed):
     region_idx = gene_idx // (
         genome.primitives.max_arity + 1
     )  # function gene + input gene addresses
-    permissible_hidden_input_gene_values = genome._determine_alternative_permissible_values_hidden_gene(
+    permissible_hidden_input_gene_values = \
+        genome._determine_alternative_permissible_values_hidden_gene(
         gene_idx, allele, region_idx
     )