diff --git a/causallearn/utils/PCUtils/SkeletonDiscovery.py b/causallearn/utils/PCUtils/SkeletonDiscovery.py
index 6cfaa492..428d2cee 100644
--- a/causallearn/utils/PCUtils/SkeletonDiscovery.py
+++ b/causallearn/utils/PCUtils/SkeletonDiscovery.py
@@ -120,8 +120,9 @@ def skeleton_discovery(
                     else:
                         if verbose:
                             print('%d dep %d | %s with p-value %f\n' % (x, y, S, p))
-                append_value(cg.sepset, x, y, tuple(sepsets))
-                append_value(cg.sepset, y, x, tuple(sepsets))
+                if (x, y) in edge_removal or not cg.G.get_edge(cg.G.nodes[x], cg.G.nodes[y]):
+                    append_value(cg.sepset, x, y, tuple(sepsets))
+                    append_value(cg.sepset, y, x, tuple(sepsets))
 
         if show_progress:
             pbar.refresh()
diff --git a/tests/TestSkeletonDiscovery.py b/tests/TestSkeletonDiscovery.py
new file mode 100644
index 00000000..10b0a5a2
--- /dev/null
+++ b/tests/TestSkeletonDiscovery.py
@@ -0,0 +1,18 @@
+from unittest import TestCase
+import numpy as np
+from causallearn.search.ConstraintBased.PC import pc
+import networkx as nx
+from causallearn.utils.cit import chisq, fisherz, gsq, kci, mv_fisherz, d_separation
+
+
+class TestSkeletonDiscovery(TestCase):
+    def test_sepset(self):
+        truth_DAG_directed_edges = {(0, 2), (1, 2), (2, 3), (2, 4)}
+
+        true_dag_netx = nx.DiGraph()
+        true_dag_netx.add_nodes_from(list(range(5)))
+        true_dag_netx.add_edges_from(truth_DAG_directed_edges)
+
+        data = np.zeros((100, len(true_dag_netx.nodes)))  # just a placeholder
+        cg = pc(data, 0.05, d_separation, True, 0, -1, true_dag=true_dag_netx)
+        assert cg.sepset[0, 2] is None
\ No newline at end of file