diff --git a/metagraph/src/annotation/annotation_converters.cpp b/metagraph/src/annotation/annotation_converters.cpp
index dc72f5248c..9505a33cae 100644
--- a/metagraph/src/annotation/annotation_converters.cpp
+++ b/metagraph/src/annotation/annotation_converters.cpp
@@ -1157,7 +1157,8 @@ void convert_to_row_diff(const std::vector<std::string> &files,
 
     if (construction_stage == RowDiffStage::CONVERT) {
         assign_anchors(graph_fname, graph_fname, out_dir, max_path_length,
-                       ".row_reduction", get_num_threads());
+                       ".row_reduction", get_num_threads(),
+                       with_values || with_coordinates);
 
         const std::string anchors_fname = graph_fname + kRowDiffAnchorExt;
         if (!fs::exists(anchors_fname)) {
@@ -1172,6 +1173,20 @@ void convert_to_row_diff(const std::vector<std::string> &files,
             return;
         }
         mem_bytes -= anchor_size;
+
+        const std::string rd_succ_fname = graph_fname + kRowDiffForkSuccExt;
+        if (!fs::exists(rd_succ_fname)) {
+            logger->error("Can't find row-diff successor bitmap at {}", rd_succ_fname);
+            exit(1);
+        }
+        uint64_t rd_succ_size = fs::file_size(rd_succ_fname);
+        if (rd_succ_size > mem_bytes) {
+            logger->warn("row-diff successor bitmap ({} MiB) is larger than"
+                         " the memory allocated ({} MiB). Reserve more RAM.",
+                         rd_succ_size >> 20, mem_bytes >> 20);
+            return;
+        }
+        mem_bytes -= rd_succ_size;
     }
 
     if (!files.size())
diff --git a/metagraph/src/annotation/binary_matrix/row_diff/row_diff.cpp b/metagraph/src/annotation/binary_matrix/row_diff/row_diff.cpp
index f4a47255b2..5e22027c1d 100644
--- a/metagraph/src/annotation/binary_matrix/row_diff/row_diff.cpp
+++ b/metagraph/src/annotation/binary_matrix/row_diff/row_diff.cpp
@@ -37,13 +37,15 @@ void IRowDiff::load_fork_succ(const std::string &filename) {
     fork_succ_.load(f);
 }
 
-std::pair<std::vector<BinaryMatrix::Row>, std::vector<std::vector<size_t>>>
+std::pair<std::vector<BinaryMatrix::Row>, std::vector<std::vector<std::pair<size_t, size_t>>>>
 IRowDiff::get_rd_ids(const std::vector<BinaryMatrix::Row> &row_ids) const {
     assert(graph_ && "graph must be loaded");
-    assert(!fork_succ_.size() || fork_succ_.size() == graph_->get_boss().get_last().size());
+    assert(!fork_succ_.size() || fork_succ_.size() == graph_->num_nodes() + 1);
 
     using Row = BinaryMatrix::Row;
 
+    const size_t RD_PATH_RESERVE_SIZE = 2;
+
     // diff rows annotating nodes along the row-diff paths
     std::vector<Row> rd_ids;
     rd_ids.reserve(row_ids.size() * RD_PATH_RESERVE_SIZE);
@@ -55,38 +57,50 @@ IRowDiff::get_rd_ids(const std::vector<BinaryMatrix::Row> &row_ids) const {
     // Truncated row-diff paths, indexes to |rd_rows|.
     // The last index in each path points to an anchor or to a row which had
     // been reached before, and thus, will be reconstructed before this one.
-    std::vector<std::vector<size_t>> rd_paths_trunc(row_ids.size());
-
-    const graph::boss::BOSS &boss = graph_->get_boss();
-    const bit_vector &rd_succ = fork_succ_.size() ? fork_succ_ : boss.get_last();
+    std::vector<std::vector<std::pair<size_t, size_t>>> rd_paths_trunc(row_ids.size());
 
     for (size_t i = 0; i < row_ids.size(); ++i) {
-        Row row = row_ids[i];
-
-        graph::boss::BOSS::edge_index boss_edge = graph_->kmer_to_boss_index(
-                graph::AnnotatedSequenceGraph::anno_to_graph_index(row));
-
-        while (true) {
-            row = graph::AnnotatedSequenceGraph::graph_to_anno_index(
-                    graph_->boss_to_kmer_index(boss_edge));
-
+        std::vector<std::pair<size_t, size_t>> &rd_path = rd_paths_trunc[i];
+
+        std::vector<size_t> path;
+        Vector<std::pair<size_t, Row>> queue;
+        queue.emplace_back(0, row_ids[i]);
+
+        while (queue.size()) {
+            size_t depth = queue.back().first;
+            Row row = queue.back().second;
+            queue.pop_back();
+            while (depth < path.size()) {
+                assert(path.size() > 1);
+                rd_path.emplace_back(*(path.rbegin() + 1), *path.rbegin());
+                path.pop_back();
+            }
             auto [it, is_new] = node_to_rd.try_emplace(row, rd_ids.size());
-            rd_paths_trunc[i].push_back(it.value());
-
+            path.push_back(it.value());
             // If a node had been reached before, we interrupt the diff path.
             // The annotation for that node will have been reconstructed earlier
             // than for other nodes in this path as well. Thus, we will start
             // reconstruction from that node and don't need its successors.
             if (!is_new)
-                break;
+                continue;
 
             rd_ids.push_back(row);
 
             if (anchor_[row])
-                break;
+                continue;
+
+            auto node = graph::AnnotatedSequenceGraph::anno_to_graph_index(row);
+            graph_->call_row_diff_successors(node, fork_succ_, [&](auto succ) {
+                queue.emplace_back(depth + 1, graph::AnnotatedSequenceGraph::graph_to_anno_index(succ));
+            });
+        }
 
-            boss_edge = boss.row_diff_successor(boss_edge, rd_succ);
+        while (path.size() > 1) {
+            rd_path.emplace_back(*(path.rbegin() + 1), *path.rbegin());
+            path.pop_back();
         }
+        assert(path.size());
+        rd_path.emplace_back(-1, path[0]);
     }
 
     return std::make_pair(std::move(rd_ids), std::move(rd_paths_trunc));
diff --git a/metagraph/src/annotation/binary_matrix/row_diff/row_diff.hpp b/metagraph/src/annotation/binary_matrix/row_diff/row_diff.hpp
index d75e206d0a..74aecb0d62 100644
--- a/metagraph/src/annotation/binary_matrix/row_diff/row_diff.hpp
+++ b/metagraph/src/annotation/binary_matrix/row_diff/row_diff.hpp
@@ -13,7 +13,6 @@
 #include "common/logger.hpp"
 #include "common/utils/template_utils.hpp"
 #include "graph/annotated_dbg.hpp"
-#include "graph/representation/succinct/boss.hpp"
 #include "graph/representation/succinct/dbg_succinct.hpp"
 
 
@@ -44,7 +43,7 @@ class IRowDiff {
 
   protected:
     // get row-diff paths starting at |row_ids|
-    std::pair<std::vector<BinaryMatrix::Row>, std::vector<std::vector<size_t>>>
+    std::pair<std::vector<BinaryMatrix::Row>, std::vector<std::vector<std::pair<size_t, size_t>>>>
     get_rd_ids(const std::vector<BinaryMatrix::Row> &row_ids) const;
 
     const graph::DBGSuccinct *graph_ = nullptr;
@@ -115,7 +114,7 @@ template <class BaseMatrix>
 bool RowDiff<BaseMatrix>::get(Row row, Column column) const {
     assert(graph_ && "graph must be loaded");
     assert(anchor_.size() == diffs_.num_rows() && "anchors must be loaded");
-    assert(!fork_succ_.size() || fork_succ_.size() == graph_->get_boss().get_last().size());
+    assert(!fork_succ_.size() || fork_succ_.size() == graph_->num_nodes() + 1);
 
     SetBitPositions set_bits = get_row(row);
     SetBitPositions::iterator v = std::lower_bound(set_bits.begin(), set_bits.end(), column);
@@ -130,9 +129,9 @@ template <class BaseMatrix>
 std::vector<BinaryMatrix::Row> RowDiff<BaseMatrix>::get_column(Column column) const {
     assert(graph_ && "graph must be loaded");
     assert(anchor_.size() == diffs_.num_rows() && "anchors must be loaded");
+    assert(!fork_succ_.size() || fork_succ_.size() == graph_->num_nodes() + 1);
 
     const graph::boss::BOSS &boss = graph_->get_boss();
-    assert(!fork_succ_.size() || fork_succ_.size() == boss.get_last().size());
 
     std::vector<Row> result;
     // TODO: implement a more efficient algorithm
@@ -151,21 +150,16 @@ template <class BaseMatrix>
 BinaryMatrix::SetBitPositions RowDiff<BaseMatrix>::get_row(Row row) const {
     assert(graph_ && "graph must be loaded");
     assert(anchor_.size() == diffs_.num_rows() && "anchors must be loaded");
-    assert(!fork_succ_.size() || fork_succ_.size() == graph_->get_boss().get_last().size());
+    assert(!fork_succ_.size() || fork_succ_.size() == graph_->num_nodes() + 1);
 
     Vector<uint64_t> result = diffs_.get_row(row);
     std::sort(result.begin(), result.end());
 
-    uint64_t boss_edge = graph_->kmer_to_boss_index(
-            graph::AnnotatedSequenceGraph::anno_to_graph_index(row));
-    const graph::boss::BOSS &boss = graph_->get_boss();
-    const bit_vector &rd_succ = fork_succ_.size() ? fork_succ_ : boss.get_last();
+    auto node = graph::AnnotatedSequenceGraph::anno_to_graph_index(row);
 
     while (!anchor_[row]) {
-        boss_edge = boss.row_diff_successor(boss_edge, rd_succ);
-
-        row = graph::AnnotatedSequenceGraph::graph_to_anno_index(
-                graph_->boss_to_kmer_index(boss_edge));
+        node = graph_->row_diff_successor(node, fork_succ_);
+        row = graph::AnnotatedSequenceGraph::graph_to_anno_index(node);
 
         auto diff_row = diffs_.get_row(row);
         std::sort(diff_row.begin(), diff_row.end());
@@ -180,7 +174,7 @@ std::vector<BinaryMatrix::SetBitPositions>
 RowDiff<BaseMatrix>::get_rows(const std::vector<Row> &row_ids) const {
     assert(graph_ && "graph must be loaded");
     assert(anchor_.size() == diffs_.num_rows() && "anchors must be loaded");
-    assert(!fork_succ_.size() || fork_succ_.size() == graph_->get_boss().get_last().size());
+    assert(!fork_succ_.size() || fork_succ_.size() == graph_->num_nodes() + 1);
 
     // diff rows annotating nodes along the row-diff paths
     std::vector<Row> rd_ids;
@@ -199,19 +193,10 @@ RowDiff<BaseMatrix>::get_rows(const std::vector<Row> &row_ids) const {
     // been reached before, and thus, will be reconstructed before this one.
     std::vector<std::vector<size_t>> rd_paths_trunc(row_ids.size());
 
-    const graph::boss::BOSS &boss = graph_->get_boss();
-    const bit_vector &rd_succ = fork_succ_.size() ? fork_succ_ : boss.get_last();
-
     for (size_t i = 0; i < row_ids.size(); ++i) {
         Row row = row_ids[i];
 
-        graph::boss::BOSS::edge_index boss_edge = graph_->kmer_to_boss_index(
-                graph::AnnotatedSequenceGraph::anno_to_graph_index(row));
-
         while (true) {
-            row = graph::AnnotatedSequenceGraph::graph_to_anno_index(
-                    graph_->boss_to_kmer_index(boss_edge));
-
             auto [it, is_new] = node_to_rd.try_emplace(row, rd_ids.size());
             rd_paths_trunc[i].push_back(it.value());
 
@@ -230,7 +215,9 @@ RowDiff<BaseMatrix>::get_rows(const std::vector<Row> &row_ids) const {
             if (anchor_[row])
                 break;
 
-            boss_edge = boss.row_diff_successor(boss_edge, rd_succ);
+            auto node = graph::AnnotatedSequenceGraph::anno_to_graph_index(row);
+            node = graph_->row_diff_successor(node, fork_succ_);
+            row = graph::AnnotatedSequenceGraph::graph_to_anno_index(node);
         }
     }
 
diff --git a/metagraph/src/annotation/int_matrix/row_diff/int_row_diff.hpp b/metagraph/src/annotation/int_matrix/row_diff/int_row_diff.hpp
index 978c48959b..fa5da63d3b 100644
--- a/metagraph/src/annotation/int_matrix/row_diff/int_row_diff.hpp
+++ b/metagraph/src/annotation/int_matrix/row_diff/int_row_diff.hpp
@@ -13,7 +13,6 @@
 #include "common/logger.hpp"
 #include "common/utils/template_utils.hpp"
 #include "graph/annotated_dbg.hpp"
-#include "graph/representation/succinct/boss.hpp"
 #include "graph/representation/succinct/dbg_succinct.hpp"
 #include "annotation/binary_matrix/row_diff/row_diff.hpp"
 #include "annotation/int_matrix/base/int_matrix.hpp"
@@ -86,9 +85,9 @@ template <class BaseMatrix>
 std::vector<IntMatrix::Row> IntRowDiff<BaseMatrix>::get_column(Column j) const {
     assert(graph_ && "graph must be loaded");
     assert(anchor_.size() == diffs_.num_rows() && "anchors must be loaded");
+    assert(!fork_succ_.size() || fork_succ_.size() == graph_->num_nodes() + 1);
 
     const graph::boss::BOSS &boss = graph_->get_boss();
-    assert(!fork_succ_.size() || fork_succ_.size() == boss.get_last().size());
 
     // TODO: implement a more efficient algorithm
     std::vector<Row> result;
@@ -113,7 +112,7 @@ std::vector<IntMatrix::RowValues>
 IntRowDiff<BaseMatrix>::get_row_values(const std::vector<Row> &row_ids) const {
     assert(graph_ && "graph must be loaded");
     assert(anchor_.size() == diffs_.num_rows() && "anchors must be loaded");
-    assert(!fork_succ_.size() || fork_succ_.size() == graph_->get_boss().get_last().size());
+    assert(!fork_succ_.size() || fork_succ_.size() == graph_->num_nodes() + 1);
 
     // get row-diff paths
     auto [rd_ids, rd_paths_trunc] = get_rd_ids(row_ids);
@@ -121,6 +120,7 @@ IntRowDiff<BaseMatrix>::get_row_values(const std::vector<Row> &row_ids) const {
     std::vector<RowValues> rd_rows = diffs_.get_row_values(rd_ids);
     for (auto &row : rd_rows) {
         decode_diffs(&row);
+        std::sort(row.begin(), row.end());
     }
 
     rd_ids = std::vector<Row>();
@@ -129,17 +129,15 @@ IntRowDiff<BaseMatrix>::get_row_values(const std::vector<Row> &row_ids) const {
     std::vector<RowValues> rows(row_ids.size());
 
     for (size_t i = 0; i < row_ids.size(); ++i) {
-        RowValues &result = rows[i];
+        const auto &rd_path = rd_paths_trunc[i];
         // propagate back and reconstruct full annotations for predecessors
-        for (auto it = rd_paths_trunc[i].rbegin(); it != rd_paths_trunc[i].rend(); ++it) {
-            std::sort(rd_rows[*it].begin(), rd_rows[*it].end());
-            add_diff(rd_rows[*it], &result);
-            // replace diff row with full reconstructed annotation
-            rd_rows[*it] = result;
+        for (size_t j = 0; j + 1 < rd_path.size(); ++j) {
+            auto [node, succ] = rd_path[j];
+            // reconstruct annotation by adding the diff (full succ + diff)
+            add_diff(rd_rows[succ], &rd_rows[node]);
         }
-        assert(std::all_of(result.begin(), result.end(),
-                           [](auto &p) { return p.second; }));
-        assert(std::all_of(result.begin(), result.end(),
+        rows[i] = rd_rows[rd_path.back().second];
+        assert(std::all_of(rows[i].begin(), rows[i].end(),
                            [](auto &p) { return (int64_t)p.second > 0; }));
     }
 
diff --git a/metagraph/src/annotation/int_matrix/row_diff/tuple_row_diff.hpp b/metagraph/src/annotation/int_matrix/row_diff/tuple_row_diff.hpp
index 107ef3ee74..45e310754e 100644
--- a/metagraph/src/annotation/int_matrix/row_diff/tuple_row_diff.hpp
+++ b/metagraph/src/annotation/int_matrix/row_diff/tuple_row_diff.hpp
@@ -13,7 +13,6 @@
 #include "common/logger.hpp"
 #include "common/utils/template_utils.hpp"
 #include "graph/annotated_dbg.hpp"
-#include "graph/representation/succinct/boss.hpp"
 #include "graph/representation/succinct/dbg_succinct.hpp"
 #include "annotation/binary_matrix/row_diff/row_diff.hpp"
 #include "annotation/int_matrix/base/int_matrix.hpp"
@@ -52,6 +51,7 @@ class TupleRowDiff : public binmat::IRowDiff, public MultiIntMatrix {
 
   private:
     static void decode_diffs(RowTuples *diffs);
+    static void shift_coords(RowTuples *diffs);
     static void add_diff(const RowTuples &diff, RowTuples *row);
 
     BaseMatrix diffs_;
@@ -69,9 +69,9 @@ template <class BaseMatrix>
 std::vector<MultiIntMatrix::Row> TupleRowDiff<BaseMatrix>::get_column(Column j) const {
     assert(graph_ && "graph must be loaded");
     assert(anchor_.size() == diffs_.num_rows() && "anchors must be loaded");
+    assert(!fork_succ_.size() || fork_succ_.size() == graph_->num_nodes() + 1);
 
     const graph::boss::BOSS &boss = graph_->get_boss();
-    assert(!fork_succ_.size() || fork_succ_.size() == boss.get_last().size());
 
     // TODO: implement a more efficient algorithm
     std::vector<Row> result;
@@ -96,7 +96,7 @@ std::vector<MultiIntMatrix::RowTuples>
 TupleRowDiff<BaseMatrix>::get_row_tuples(const std::vector<Row> &row_ids) const {
     assert(graph_ && "graph must be loaded");
     assert(anchor_.size() == diffs_.num_rows() && "anchors must be loaded");
-    assert(!fork_succ_.size() || fork_succ_.size() == graph_->get_boss().get_last().size());
+    assert(!fork_succ_.size() || fork_succ_.size() == graph_->num_nodes() + 1);
 
     // get row-diff paths
     auto [rd_ids, rd_paths_trunc] = get_rd_ids(row_ids);
@@ -104,6 +104,7 @@ TupleRowDiff<BaseMatrix>::get_row_tuples(const std::vector<Row> &row_ids) const
     std::vector<RowTuples> rd_rows = diffs_.get_row_tuples(rd_ids);
     for (auto &row : rd_rows) {
         decode_diffs(&row);
+        std::sort(row.begin(), row.end());
     }
 
     rd_ids = std::vector<Row>();
@@ -112,19 +113,30 @@ TupleRowDiff<BaseMatrix>::get_row_tuples(const std::vector<Row> &row_ids) const
     std::vector<RowTuples> rows(row_ids.size());
 
     for (size_t i = 0; i < row_ids.size(); ++i) {
-        RowTuples &result = rows[i];
-
-        auto it = rd_paths_trunc[i].rbegin();
-        std::sort(rd_rows[*it].begin(), rd_rows[*it].end());
-        result = rd_rows[*it];
+        const auto &rd_path = rd_paths_trunc[i];
+        uint64_t last_node = -1;
         // propagate back and reconstruct full annotations for predecessors
-        for (++it ; it != rd_paths_trunc[i].rend(); ++it) {
-            std::sort(rd_rows[*it].begin(), rd_rows[*it].end());
-            add_diff(rd_rows[*it], &result);
-            // replace diff row with full reconstructed annotation
-            rd_rows[*it] = result;
+        for (size_t j = 0; j + 1 < rd_path.size(); ++j) {
+            auto [node, succ] = rd_path[j];
+            if (last_node == succ)
+                shift_coords(&rd_rows[last_node]);
+            // reconstruct annotation by adding the diff (full succ + diff)
+            add_diff(rd_rows[succ], &rd_rows[node]);
+            last_node = node;
+        }
+        if (last_node != (uint64_t)-1) {
+            assert(last_node == rd_path.back().second);
+            shift_coords(&rd_rows[last_node]);
+        }
+        auto &row = rd_rows[rd_path.back().second];
+        rows[i].reserve(std::count_if(row.begin(), row.end(),
+                                      [](const auto &v) { return !v.second.empty(); }));
+        for (auto&& v : row) {
+            if (v.second.size())
+                rows[i].push_back(std::move(v));
         }
-        assert(std::all_of(result.begin(), result.end(),
+        row = rows[i];
+        assert(std::all_of(rows[i].begin(), rows[i].end(),
                            [](auto &p) { return p.second.size(); }));
     }
 
@@ -152,6 +164,16 @@ void TupleRowDiff<BaseMatrix>::decode_diffs(RowTuples *diffs) {
     // no encoding
 }
 
+template <class BaseMatrix>
+void TupleRowDiff<BaseMatrix>::shift_coords(RowTuples *diffs) {
+    for (auto &[j, tuple] : *diffs) {
+        assert(std::is_sorted(tuple.begin(), tuple.end()));
+        for (uint64_t &c : tuple) {
+            c -= SHIFT;
+        }
+    }
+}
+
 template <class BaseMatrix>
 void TupleRowDiff<BaseMatrix>::add_diff(const RowTuples &diff, RowTuples *row) {
     assert(std::is_sorted(row->begin(), row->end()));
@@ -168,32 +190,34 @@ void TupleRowDiff<BaseMatrix>::add_diff(const RowTuples &diff, RowTuples *row) {
                 result.push_back(*it);
                 ++it;
             } else if (it->first > it2->first) {
-                result.push_back(*it2);
+                if (it2->second.size())
+                    result.push_back(*it2);
                 ++it2;
             } else {
-                if (it2->second.size()) {
-                    result.emplace_back(it->first, Tuple{});
+                result.emplace_back(it->first, Tuple{});
+                if (it->second.size()) {
+                    assert(std::is_sorted(it->second.begin(), it->second.end()));
+                    assert(std::is_sorted(it2->second.begin(), it2->second.end()));
                     std::set_symmetric_difference(it->second.begin(), it->second.end(),
                                                   it2->second.begin(), it2->second.end(),
                                                   std::back_inserter(result.back().second));
+                    if (result.back().second.empty())
+                        result.pop_back();
                 }
                 ++it;
                 ++it2;
             }
         }
         std::copy(it, row->end(), std::back_inserter(result));
-        std::copy(it2, diff.end(), std::back_inserter(result));
+        for ( ; it2 != diff.end(); ++it2) {
+            if (it2->second.size())
+                result.push_back(*it2);
+        }
 
         row->swap(result);
     }
 
     assert(std::is_sorted(row->begin(), row->end()));
-    for (auto &[j, tuple] : *row) {
-        assert(std::is_sorted(tuple.begin(), tuple.end()));
-        for (uint64_t &c : tuple) {
-            c -= SHIFT;
-        }
-    }
 }
 
 } // namespace matrix
diff --git a/metagraph/src/annotation/row_diff_builder.cpp b/metagraph/src/annotation/row_diff_builder.cpp
index 000cf8601d..a925cf4f30 100644
--- a/metagraph/src/annotation/row_diff_builder.cpp
+++ b/metagraph/src/annotation/row_diff_builder.cpp
@@ -258,13 +258,14 @@ void sum_and_call_counts(const fs::path &dir,
     }
 }
 
-rd_succ_bv_type route_at_forks(const graph::DBGSuccinct &graph,
-                               const std::string &rd_succ_filename,
-                               const std::string &count_vectors_dir,
-                               const std::string &row_count_extension) {
+void route_at_forks(const graph::DBGSuccinct &graph,
+                    const std::string &rd_succ_fname,
+                    const std::string &count_vectors_dir,
+                    const std::string &row_count_extension) {
     logger->trace("Assigning row-diff successors at forks...");
 
-    rd_succ_bv_type rd_succ;
+    const bit_vector &last = graph.get_boss().get_last();
+    sdsl::bit_vector rd_succ_bv(graph.num_nodes() + 1, false);
 
     bool optimize_forks = false;
     for (const auto &p : fs::directory_iterator(count_vectors_dir)) {
@@ -276,13 +277,9 @@ rd_succ_bv_type route_at_forks(const graph::DBGSuccinct &graph,
         logger->trace("RowDiff successors will be set to the adjacent nodes with"
                       " the largest number of labels");
 
-        const bit_vector &last = graph.get_boss().get_last();
         graph::DeBruijnGraph::node_index graph_idx = to_node(0);
-
         std::vector<uint32_t> outgoing_counts;
 
-        sdsl::bit_vector rd_succ_bv(last.size(), false);
-
         sum_and_call_counts(count_vectors_dir, row_count_extension, "row counts",
             [&](int32_t count) {
                 // TODO: skip single outgoing
@@ -292,7 +289,7 @@ rd_succ_bv_type route_at_forks(const graph::DBGSuccinct &graph,
                     size_t max_pos = std::max_element(outgoing_counts.rbegin(),
                                                       outgoing_counts.rend())
                                      - outgoing_counts.rbegin();
-                    rd_succ_bv[graph.kmer_to_boss_index(graph_idx - max_pos)] = true;
+                    rd_succ_bv[graph_idx - max_pos] = true;
                     outgoing_counts.resize(0);
                 }
                 graph_idx++;
@@ -305,19 +302,23 @@ rd_succ_bv_type route_at_forks(const graph::DBGSuccinct &graph,
             exit(1);
         }
 
-        rd_succ = rd_succ_bv_type(std::move(rd_succ_bv));
-
     } else {
+        // TODO: remove this mode?
         logger->warn("No count vectors could be found in {}. The last outgoing"
                      " edges will be selected for assigning RowDiff successors",
                      count_vectors_dir);
+
+        last.call_ones([&](BOSS::edge_index i) {
+            rd_succ_bv[graph.boss_to_kmer_index(i)] = true;
+        });
+        // npos is never a successor
+        rd_succ_bv[0] = false;
     }
 
-    std::ofstream f(rd_succ_filename, ios::binary);
-    rd_succ.serialize(f);
+    std::ofstream f(rd_succ_fname, ios::binary);
+    rd_succ_bv.serialize(f);
     logger->trace("RowDiff successors are assigned for forks and written to {}",
-                  rd_succ_filename);
-    return rd_succ;
+                  rd_succ_fname);
 }
 
 void build_pred_succ(const std::string &graph_fname,
@@ -344,8 +345,8 @@ void build_pred_succ(const std::string &graph_fname,
     }
 
     // assign row-diff successors at forks
-    rd_succ_bv_type rd_succ = route_at_forks(graph, outfbase + kRowDiffForkSuccExt,
-                                             count_vectors_dir, row_count_extension);
+    route_at_forks(graph, outfbase + kRowDiffForkSuccExt,
+                   count_vectors_dir, row_count_extension);
 
     const BOSS &boss = graph.get_boss();
 
@@ -379,28 +380,22 @@ void build_pred_succ(const std::string &graph_fname,
                 BOSS::edge_index next = boss.fwd(boss_idx, d);
                 assert(next);
                 if (!dummy[next]) {
-                    while (rd_succ.size() && !rd_succ[next]) {
-                        next--;
-                        assert(!boss.get_last(next));
-                    }
-                    succ_buf.push_back(to_row(graph.boss_to_kmer_index(next)));
-                    succ_boundary_buf.push_back(0);
+                    do {
+                        succ_buf.push_back(to_row(graph.boss_to_kmer_index(next)));
+                        succ_boundary_buf.push_back(0);
+                    } while (!boss.get_last(--next));
                 }
-                // compute predecessors only for row-diff successors
-                if (rd_succ.size() ? rd_succ[boss_idx] : boss.get_last(boss_idx)) {
-                    BOSS::TAlphabet d = boss.get_node_last_value(boss_idx);
-                    BOSS::edge_index back_idx = boss.bwd(boss_idx);
-                    boss.call_incoming_to_target(back_idx, d,
-                        [&](BOSS::edge_index pred) {
-                            // dummy predecessors are ignored
-                            if (!dummy[pred]) {
-                                uint64_t node_index = graph.boss_to_kmer_index(pred);
-                                pred_buf.push_back(to_row(node_index));
-                                pred_boundary_buf.push_back(0);
-                            }
+                BOSS::edge_index back_idx = boss.bwd(boss_idx);
+                boss.call_incoming_to_target(back_idx, boss.get_node_last_value(boss_idx),
+                    [&](BOSS::edge_index pred) {
+                        // dummy predecessors are ignored
+                        if (!dummy[pred]) {
+                            uint64_t node_index = graph.boss_to_kmer_index(pred);
+                            pred_buf.push_back(to_row(node_index));
+                            pred_boundary_buf.push_back(0);
                         }
-                    );
-                }
+                    }
+                );
             }
             succ_boundary_buf.push_back(1);
             pred_boundary_buf.push_back(1);
@@ -425,7 +420,8 @@ void assign_anchors(const std::string &graph_fname,
                     const std::filesystem::path &count_vectors_dir,
                     uint32_t max_length,
                     const std::string &row_reduction_extension,
-                    uint32_t num_threads) {
+                    uint32_t num_threads,
+                    bool multiple_fork_successors) {
     std::string anchor_filename = outfbase + kRowDiffAnchorExt;
     if (fs::exists(anchor_filename)) {
         logger->trace("Using existing anchors {}", anchor_filename);
@@ -478,27 +474,43 @@ void assign_anchors(const std::string &graph_fname,
     // assign extra anchors and restrict the length of row-diff paths
     logger->trace("Assigning required anchors...");
     {
-        rd_succ_bv_type rd_succ;
+        sdsl::bit_vector rd_succ_bv;
         const std::string &rd_succ_fname = outfbase + kRowDiffForkSuccExt;
         std::ifstream f(rd_succ_fname, ios::binary);
-        if (!rd_succ.load(f)) {
+        try {
+            rd_succ_bv.load(f);
+        } catch (...) {
             logger->error("Couldn't load row-diff successor bitmap from {}", rd_succ_fname);
             exit(1);
         }
-
-        if (rd_succ.size()) {
-            logger->trace("Assigning anchors for RowDiff successors {}...", rd_succ_fname);
-            boss.row_diff_traverse(num_threads, max_length, rd_succ, &anchors_bv);
-        } else {
-            logger->warn("Assigning anchors without chosen RowDiff successors."
-                         " The last outgoing edges will be used for routing.");
-            boss.row_diff_traverse(num_threads, max_length, boss.get_last(), &anchors_bv);
+        if (rd_succ_bv.size() != graph.num_nodes() + 1) {
+            logger->error("Successor bitmap {} is incompatible with the graph."
+                          " Vector size: {}, number of nodes: {}",
+                          rd_succ_fname, rd_succ_bv.size(), graph.num_nodes());
+            exit(1);
         }
+
+        logger->trace("Assigning anchors for RowDiff successors {}...", rd_succ_fname);
+        boss.row_diff_traverse(num_threads, max_length,
+                               [&](BOSS::edge_index i) { return rd_succ_bv[graph.boss_to_kmer_index(i)]; },
+                               &anchors_bv);
+
+        logger->trace("Adding branching off forks to RowDiff successors {}...", rd_succ_fname);
+        const uint64_t num_fork_successors = sdsl::util::cnt_one_bits(rd_succ_bv);
+        if (multiple_fork_successors)
+            graph.add_rd_successors_at_forks(num_threads, anchors_bv, &rd_succ_bv, 10 * max_length);
+
+        rd_succ_bv_type rd_succ(std::move(rd_succ_bv));
+        logger->trace("Number of successors at forks increased from {} to {}",
+                      num_fork_successors, rd_succ.num_set_bits());
+        std::ofstream out(rd_succ_fname, ios::binary);
+        rd_succ.serialize(out);
+        logger->trace("Updated RowDiff successors {}", rd_succ_fname);
     }
 
     // anchors_bv uses BOSS edges as indices, so we need to map it to annotation indices
     {
-        sdsl::bit_vector anchors(num_rows, false);
+        sdsl::bit_vector anchors(num_rows, 0);
         for (BOSS::edge_index i = 1; i < anchors_bv.size(); ++i) {
             if (anchors_bv[i]) {
                 uint64_t graph_idx = graph.boss_to_kmer_index(i);
@@ -536,7 +548,7 @@ using CallOnes = std::function<void(const bit_vector &source_col,
                                     uint64_t row_idx_chunk,
                                     size_t source_idx,
                                     size_t col_idx,
-                                    const uint64_t *succ,
+                                    const uint64_t *next_begin, const uint64_t *next_end,
                                     const uint64_t *pred_begin, const uint64_t *pred_end)>;
 
 void read_next_block(sdsl::int_vector_buffer<>::iterator &it,
@@ -657,12 +669,10 @@ void traverse_anno_chunked(
                 source_col.call_ones_in_range(chunk, chunk + block_size,
                     [&](uint64_t i) {
                         assert(succ_chunk_idx[i - chunk + 1] >= succ_chunk_idx[i - chunk]);
-                        assert(succ_chunk_idx[i - chunk + 1] <= succ_chunk_idx[i - chunk] + 1);
-                        const uint64_t *succ = succ_chunk_idx[i - chunk + 1]
-                                                > succ_chunk_idx[i - chunk]
-                                                ? succ_chunk.data() + succ_chunk_idx[i - chunk]
-                                                : NULL;
-                        call_ones(source_col, i, i - chunk, l_idx, j, succ,
+                        assert(pred_chunk_idx[i - chunk + 1] >= pred_chunk_idx[i - chunk]);
+                        call_ones(source_col, i, i - chunk, l_idx, j,
+                                  succ_chunk.data() + succ_chunk_idx[i - chunk],
+                                  succ_chunk.data() + succ_chunk_idx[i - chunk + 1],
                                   pred_chunk.data() + pred_chunk_idx[i - chunk],
                                   pred_chunk.data() + pred_chunk_idx[i - chunk + 1]);
                     }
@@ -782,6 +792,7 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
     }
 
     anchor_bv_type anchor;
+    rd_succ_bv_type rd_succ;
     if (!compute_row_reduction) {
         const std::string anchors_fname = pred_succ_fprefix + kRowDiffAnchorExt;
         std::ifstream f(anchors_fname, std::ios::binary);
@@ -795,6 +806,19 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
                           anchors_fname, anchor.size(), num_rows);
             exit(1);
         }
+
+        const std::string &rd_succ_fname = pred_succ_fprefix + kRowDiffForkSuccExt;
+        f = std::ifstream(rd_succ_fname, ios::binary);
+        if (!rd_succ.load(f)) {
+            logger->error("Couldn't load row-diff successor bitmap from {}", rd_succ_fname);
+            exit(1);
+        }
+        if (rd_succ.size() != num_rows + 1) {
+            logger->error("Successor bitmap {} is incompatible with annotations."
+                          " Vector size: {}, number of rows: {}",
+                          rd_succ_fname, rd_succ.size(), num_rows);
+            exit(1);
+        }
     }
 
     const bool swap_disk = !swap_dir.empty();
@@ -825,7 +849,7 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
         v.resize(0);
     };
 
-    std::function<std::function<void(const std::function<void(uint64_t)> &)>(size_t, size_t)> call_diffs;
+    std::function<std::vector<T>(size_t, size_t)> get_diff_values;
     if (swap_disk) {
         uint64_t total_num_labels = 0;
         for (size_t s = 0; s < sources.size(); ++s) {
@@ -854,46 +878,61 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
             exit(1);
         }
 
-        call_diffs = [&,chunks_open_per_thread](size_t s, size_t j) {
-            return [&,s,j](const std::function<void(uint64_t)> &call) {
-                std::vector<std::string> filenames;
-                // for stage 1, fwd bits are already counted, so we skip that chunk
-                for (uint32_t chunk = compute_row_reduction ? 1 : 0;
-                                    chunk < num_chunks[s][j]; ++chunk) {
-                    filenames.push_back(tmp_file(s, j, chunk));
-                }
+        get_diff_values = [&,chunks_open_per_thread](size_t s, size_t j) {
+            std::vector<std::string> filenames;
+            // for stage 1, fwd bits are already counted, so we skip that chunk
+            for (uint32_t chunk = compute_row_reduction ? 1 : 0;
+                                chunk < num_chunks[s][j]; ++chunk) {
+                filenames.push_back(tmp_file(s, j, chunk));
+            }
 
-                const bool remove_chunks = true;
-                uint64_t r = 0;
+            std::vector<T> buffer;
+            const bool remove_chunks = true;
+            if constexpr(with_values) {
                 elias_fano::merge_files<T>(filenames, [&](T v) {
-                    call(utils::get_first(v));
-                    if constexpr(with_values) {
-                        assert(v.second && "zero diffs must have been skipped");
-                        values[s][j][r++] = matrix::encode_diff(v.second);
+                    assert(v.second && "zero diffs must have been skipped");
+                    if (buffer.size() && buffer.back().first == v.first) {
+                        buffer.back().second += v.second;
+                    } else {
+                        buffer.push_back(v);
                     }
                 }, remove_chunks, chunks_open_per_thread);
-            };
+            } else {
+                elias_fano::merge_files<T>(filenames, [&](T v) { buffer.push_back(v); },
+                                           remove_chunks, chunks_open_per_thread, true /* dedupe */);
+            }
+            return buffer;
         };
     } else {
         logger->info("Diff-transform in memory without disk swap");
 
-        call_diffs = [&](size_t s, size_t j) {
-            return [&,s,j](const std::function<void(uint64_t)> &call) {
-                auto &v = set_rows_fwd[s][j];
-                v.insert(v.end(), set_rows_bwd[s][j].begin(), set_rows_bwd[s][j].end());
-                set_rows_bwd[s][j] = {};
+        get_diff_values = [&](size_t s, size_t j) {
+            auto &buffer = set_rows_fwd[s][j];
+            assert(std::is_sorted(buffer.begin(), buffer.end()));
 
-                std::sort(v.begin(), v.end());
+            buffer.insert(buffer.end(), set_rows_bwd[s][j].begin(), set_rows_bwd[s][j].end());
+            set_rows_bwd[s][j] = {};
 
+            std::sort(buffer.begin(), buffer.end());
+
+            if constexpr(with_values) {
                 uint64_t r = 0;
-                for (size_t i = 0; i < v.size(); ++i) {
-                    call(utils::get_first(v[i]));
-                    if constexpr(with_values) {
-                        assert(v[i].second && "zero diffs must have been skipped");
-                        values[s][j][r++] = matrix::encode_diff(v[i].second);
+                for (size_t i = 1; i < buffer.size(); ++i) {
+                    assert(buffer[i].second && "zero diffs must have been skipped");
+                    if (buffer[i].first == buffer[r].first) {
+                        buffer[r].second += buffer[i].second;
+                    } else {
+                        buffer[++r] = buffer[i];
                     }
                 }
-            };
+                if (buffer.size())
+                    buffer.resize(r + 1);
+
+            } else {
+                buffer.erase(std::unique(buffer.begin(), buffer.end()), buffer.end());
+            }
+
+            return std::move(buffer);
         };
     }
 
@@ -962,15 +1001,22 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
 
     // get bit at position |i| or its value
     auto get_value = [&](const bit_vector &col,
-                         size_t s, size_t j, uint64_t i) -> uint64_t {
+                         size_t s, size_t j,
+                         const uint64_t *it, const uint64_t *end) -> uint64_t {
         if (!with_values) {
-            return col[i];
+            for ( ; it != end; ++it) {
+                if ((compute_row_reduction || rd_succ[to_node(*it)]) && col[*it])
+                    return true;
+            }
+            return false;
         } else {
-            if (uint64_t rk = col.conditional_rank1(i)) {
-                return values[s][j][rk - 1];
-            } else {
-                return 0;
+            uint64_t value = 0;
+            uint64_t rk;
+            for ( ; it != end; ++it) {
+                if ((compute_row_reduction || rd_succ[to_node(*it)]) && (rk = col.conditional_rank1(*it)))
+                    value += values[s][j][rk - 1];
             }
+            return value;
         }
     };
 
@@ -982,7 +1028,7 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
             },
             [&](const bit_vector &source_col, uint64_t row_idx, uint64_t chunk_idx,
                     size_t source_idx, size_t j,
-                    const uint64_t *succ,
+                    const uint64_t *succ_begin, const uint64_t *succ_end,
                     const uint64_t *pred_begin, const uint64_t *pred_end) {
 
                 // get bits for these positions (or values, hence uint64_t)
@@ -994,15 +1040,20 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
                 }
 
                 if (compute_row_reduction) {
-                    if (succ && curr_value == get_value(source_col, source_idx, j, *succ)) {
-                        // reduction (zero diff)
-                        __atomic_add_fetch(&row_nbits_block[chunk_idx], 1, __ATOMIC_RELAXED);
+                    if (succ_begin != succ_end) {
+                        const uint64_t succ_value = get_value(source_col, source_idx, j, succ_begin, succ_end);
+                        bool before = curr_value; // if anchor
+                        bool after = (curr_value != succ_value); // if diff
+                        // reduction
+                        __atomic_add_fetch(&row_nbits_block[chunk_idx],
+                                           (int)before - (int)after, __ATOMIC_RELAXED);
                     }
                 } else {
                     bool is_anchor = anchor[row_idx];
                     // add current bit if this node is an anchor
                     // or if the successor has zero diff
-                    uint64_t succ_value = is_anchor ? 0 : get_value(source_col, source_idx, j, *succ);
+                    uint64_t succ_value = is_anchor ? 0 : get_value(source_col, source_idx,
+                                                                    j, succ_begin, succ_end);
                     if (succ_value != curr_value) {
                         // no reduction, we must keep the bit
                         auto &v = set_rows_fwd[source_idx][j];
@@ -1019,9 +1070,12 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
                     }
                 }
 
+                if (!curr_value || !(compute_row_reduction || rd_succ[to_node(row_idx)]))
+                    return;
+
                 // check non-anchor predecessor nodes and add them if they are zero
                 for (const uint64_t *pred_p = pred_begin; pred_p < pred_end; ++pred_p) {
-                    if (curr_value && !source_col[*pred_p] && (compute_row_reduction || !anchor[*pred_p])) {
+                    if (!source_col[*pred_p] && (compute_row_reduction || !anchor[*pred_p])) {
                         auto &v = set_rows_bwd[source_idx][j];
                         if constexpr(with_values) {
                             v.emplace_back(*pred_p, -curr_value);
@@ -1104,14 +1158,22 @@ void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
         std::vector<std::unique_ptr<bit_vector>> columns(label_encoders[l_idx].size());
 
         for (size_t j = 0; j < label_encoders[l_idx].size(); ++j) {
+            std::vector<T> buffer = get_diff_values(l_idx, j);
             if constexpr(with_values) {
                 // diff values may be negative, hence we need wider integers
-                values[l_idx][j] = sdsl::int_vector<>(row_diff_bits[l_idx][j], 0,
+                values[l_idx][j] = sdsl::int_vector<>(buffer.size(), 0,
                                                       std::min(values[l_idx][j].width() + 1, 64));
+                for (size_t i = 0; i < buffer.size(); ++i) {
+                    values[l_idx][j][i] = matrix::encode_diff(buffer[i].second);
+                }
             }
-
-            columns[j] = std::make_unique<bit_vector_smart>(call_diffs(l_idx, j), num_rows,
-                                                            row_diff_bits[l_idx][j]);
+            columns[j] = std::make_unique<bit_vector_smart>(
+                            [&](auto call_index) {
+                                for (const auto &v : buffer) {
+                                    call_index(utils::get_first(v));
+                                }
+                            },
+                            num_rows, buffer.size());
         }
 
         if (compute_row_reduction) {
@@ -1222,15 +1284,34 @@ void convert_batch_to_row_diff_coord(const std::string &pred_succ_fprefix,
     );
     logger->trace("Done loading coordinates");
 
-    // get bit at position |i| or its value
-    auto get_value = [&](const bit_vector &col,
-                         size_t s, size_t j, uint64_t i) {
+    rd_succ_bv_type rd_succ;
+
+    // get coordinates at position |i|
+    auto get_value = [&](const bit_vector &col, size_t s, size_t j, uint64_t i) {
         std::vector<uint64_t> result;
         if (uint64_t rk = col.conditional_rank1(i)) {
             uint64_t t = delims[s][j].select1(rk);
             while (!delims[s][j][++t]) {
                 result.push_back(coords[s][j][t - rk]);
             }
+            assert(result.size());
+        }
+        assert(std::is_sorted(result.begin(), result.end()));
+        return result;
+    };
+
+    auto get_succ = [&](const bit_vector &col,
+                        size_t s, size_t j,
+                        const uint64_t *it, const uint64_t *end) {
+        std::vector<uint64_t> result;
+        uint64_t rk;
+        for ( ; it != end; ++it) {
+            if ((compute_row_reduction || rd_succ[to_node(*it)]) && (rk = col.conditional_rank1(*it))) {
+                uint64_t t = delims[s][j].select1(rk);
+                while (!delims[s][j][++t]) {
+                    result.push_back(coords[s][j][t - rk]);
+                }
+            }
         }
         std::sort(result.begin(), result.end());
         return result;
@@ -1297,17 +1378,17 @@ void convert_batch_to_row_diff_coord(const std::string &pred_succ_fprefix,
                 },
                 [&](const bit_vector &source_col, uint64_t row_idx, uint64_t chunk_idx,
                         size_t s, size_t j,
-                        const uint64_t *succ, const uint64_t *, const uint64_t *) {
-                    if (!succ)
+                        const uint64_t *succ_begin, const uint64_t *succ_end,
+                        const uint64_t *, const uint64_t *) {
+                    if (succ_begin == succ_end)
                         return;
 
                     // get annotated coordinates for this k-mer
                     const auto curr_value = get_value(source_col, s, j, row_idx);
-                    const auto diff = get_diff(curr_value, get_value(source_col, s, j, *succ));
+                    const auto diff = get_diff(curr_value, get_succ(source_col, s, j, succ_begin, succ_end));
                     // reduction (zero diff)
                     __atomic_add_fetch(&row_nbits_block[chunk_idx],
-                                       curr_value.size() - diff.size(),
-                                       __ATOMIC_RELAXED);
+                                       curr_value.size() - diff.size(), __ATOMIC_RELAXED);
                 },
                 [&](uint64_t block_begin) {
                     __atomic_thread_fence(__ATOMIC_ACQUIRE);
@@ -1331,9 +1412,22 @@ void convert_batch_to_row_diff_coord(const std::string &pred_succ_fprefix,
         return;
     }
 
+    const std::string &rd_succ_fname = pred_succ_fprefix + kRowDiffForkSuccExt;
+    std::ifstream f(rd_succ_fname, ios::binary);
+    if (!rd_succ.load(f)) {
+        logger->error("Couldn't load row-diff successor bitmap from {}", rd_succ_fname);
+        exit(1);
+    }
+    if (rd_succ.size() != num_rows + 1) {
+        logger->error("Successor bitmap {} is incompatible with annotations."
+                      " Vector size: {}, number of rows: {}",
+                      rd_succ_fname, rd_succ.size(), num_rows);
+        exit(1);
+    }
+
     anchor_bv_type anchor;
     const std::string anchors_fname = pred_succ_fprefix + kRowDiffAnchorExt;
-    std::ifstream f(anchors_fname, std::ios::binary);
+    f = std::ifstream(anchors_fname, std::ios::binary);
     if (!anchor.load(f)) {
         logger->error("Can't load anchors from {}", anchors_fname);
         exit(1);
@@ -1351,7 +1445,6 @@ void convert_batch_to_row_diff_coord(const std::string &pred_succ_fprefix,
     // edges, for each of the sources, per chunk. set_rows_fwd is already sorted
     std::vector<std::vector<std::vector<T>>> set_rows_bwd(sources.size());
     std::vector<std::vector<std::vector<T>>> set_rows_fwd(sources.size());
-    std::vector<std::vector<uint64_t>> row_diff_bits(sources.size());
     std::vector<std::vector<uint64_t>> row_diff_coords(sources.size());
     std::vector<std::vector<uint64_t>> num_coords_anchored(sources.size());
     std::vector<std::vector<uint64_t>> num_chunks(sources.size());
@@ -1372,7 +1465,6 @@ void convert_batch_to_row_diff_coord(const std::string &pred_succ_fprefix,
     for (size_t s = 0; s < sources.size(); ++s) {
         set_rows_fwd[s].resize(sources[s].num_labels());
         set_rows_bwd[s].resize(sources[s].num_labels());
-        row_diff_bits[s].assign(sources[s].num_labels(), 0);
         row_diff_coords[s].assign(sources[s].num_labels(), 0);
         num_coords_anchored[s].assign(sources[s].num_labels(), 0);
         // The first chunk will contain forward bits, all sorted.
@@ -1392,14 +1484,14 @@ void convert_batch_to_row_diff_coord(const std::string &pred_succ_fprefix,
 
     // We use this dummy index for an optimization where we don't store diff
     // for non-anchor k-mers with no coordinates.
-    uint64_t DUMMY_COORD = -1;
+    const uint64_t DUMMY_COORD = -1;
 
     traverse_anno_chunked(
             num_rows, pred_succ_fprefix, sources,
             [&](uint64_t) {},
             [&](const bit_vector &source_col, uint64_t row_idx, uint64_t,
                     size_t s, size_t j,
-                    const uint64_t *succ,
+                    const uint64_t *succ_begin, const uint64_t *succ_end,
                     const uint64_t *pred_begin, const uint64_t *pred_end) {
                 // get annotated coordinates for this k-mer
                 const auto curr_value = get_value(source_col, s, j, row_idx);
@@ -1409,15 +1501,15 @@ void convert_batch_to_row_diff_coord(const std::string &pred_succ_fprefix,
                 if (is_anchor)
                     num_coords_anchored[s][j] += curr_value.size();
 
-                const auto diff = is_anchor
+                const auto value = is_anchor
                     ? curr_value
-                    : get_diff(curr_value, get_value(source_col, s, j, *succ));
+                    : get_diff(curr_value, get_succ(source_col, s, j, succ_begin, succ_end));
 
-                if (diff.size()) {
+                if (value.size()) {
                     // must write the coordinates/diff
                     auto &v = set_rows_fwd[s][j];
 
-                    for (uint64_t coord : diff) {
+                    for (uint64_t coord : value) {
                         assert((!v.size() || v.back() != std::make_pair(row_idx, coord))
                                && "coordinates must be unique and can't repeat");
                         v.emplace_back(row_idx, coord);
@@ -1429,16 +1521,22 @@ void convert_batch_to_row_diff_coord(const std::string &pred_succ_fprefix,
                             v.resize(0);
                         }
                     }
-                    row_diff_bits[s][j]++;
                 }
 
+                if (!curr_value.size() || !rd_succ[to_node(row_idx)])
+                    return;
+
                 // check non-anchor predecessor nodes and add them if they are zero
                 for (const uint64_t *pred_p = pred_begin; pred_p < pred_end; ++pred_p) {
-                    if (curr_value.size() && !source_col[*pred_p] && !anchor[*pred_p]) {
+                    if (!source_col[*pred_p] && !anchor[*pred_p]
+                            && (!num_coords_per_seq
+                                || std::any_of(curr_value.begin(), curr_value.end(),
+                                               [&](uint64_t c) { return c % num_coords_per_seq; }))) {
                         auto &v = set_rows_bwd[s][j];
                         // indicate that there are no coordinates for the predecessor
+                        // FYI: in case of multiple successors at forks, this adds duplicate
+                        //      DUMMY_COORD from each fork-successor the empty predecessor.
                         v.emplace_back(*pred_p, DUMMY_COORD);
-                        row_diff_bits[s][j]++;
 
                         if (v.size() == v.capacity()) {
                             std::sort(v.begin(), v.end());
@@ -1511,10 +1609,9 @@ void convert_batch_to_row_diff_coord(const std::string &pred_succ_fprefix,
             // diff values may be negative, hence we need wider integers
             sdsl::int_vector<> diff_coords(row_diff_coords[l_idx][j]);
             auto coords_it = diff_coords.begin();
-            sdsl::bit_vector diff_delims(diff_coords.size() + row_diff_bits[l_idx][j] + 1, 0);
-            auto delims_it = diff_delims.begin();
-
-            auto call_ones = [&](const std::function<void(uint64_t)> &call) {
+            std::vector<uint64_t> ids;
+            std::vector<bool> diff_delims;
+            {
                 std::vector<std::string> filenames;
                 // skip chunk with fwd bits which have already been counted if stage 1
                 for (uint32_t chunk = 0; chunk < num_chunks[l_idx][j]; ++chunk) {
@@ -1525,22 +1622,23 @@ void convert_batch_to_row_diff_coord(const std::string &pred_succ_fprefix,
                 elias_fano::merge_files<T>(filenames, [&](T pair) {
                     const auto &[i, coord] = pair;
                     if (i != last) {
-                        call(i);
+                        ids.push_back(i);
                         last = i;
-                        *delims_it++ = 1;
+                        diff_delims.push_back(1);
                     }
                     if (coord != DUMMY_COORD) {
                         *coords_it++ = coord;
-                        ++delims_it;
+                        diff_delims.push_back(0);
                     }
-                }, true, chunks_open_per_thread);
-                *delims_it++ = 1;
+                }, true, chunks_open_per_thread, false);
+                diff_delims.push_back(1);
                 assert(coords_it == diff_coords.end());
-                assert(delims_it == diff_delims.end());
-            };
-            columns[j] = std::make_unique<bit_vector_smart>(call_ones, num_rows,
-                                                            row_diff_bits[l_idx][j]);
-            bit_vector_smart(std::move(diff_delims)).serialize(out_coord);
+                assert(diff_delims.size() == diff_coords.size() + ids.size() + 1);
+            }
+            columns[j] = std::make_unique<bit_vector_smart>(
+                    [&](auto callback) { std::for_each(ids.begin(), ids.end(), callback); },
+                    num_rows, ids.size());
+            bit_vector_smart(to_sdsl(diff_delims)).serialize(out_coord);
             sdsl::util::bit_compress(diff_coords);
             diff_coords.serialize(out_coord);
         }
diff --git a/metagraph/src/annotation/row_diff_builder.hpp b/metagraph/src/annotation/row_diff_builder.hpp
index f57fe4c38c..d9ae8dd1ed 100644
--- a/metagraph/src/annotation/row_diff_builder.hpp
+++ b/metagraph/src/annotation/row_diff_builder.hpp
@@ -27,7 +27,8 @@ void assign_anchors(const std::string &graph_filename,
                     const std::filesystem::path &dest_dir,
                     uint32_t max_length,
                     const std::string &row_reduction_extension,
-                    uint32_t num_threads);
+                    uint32_t num_threads,
+                    bool multiple_fork_successors);
 
 void convert_batch_to_row_diff(const std::string &pred_succ_fprefix,
                                const std::vector<std::string> &source_files,
diff --git a/metagraph/src/common/elias_fano/elias_fano_merger.hpp b/metagraph/src/common/elias_fano/elias_fano_merger.hpp
index f42298c470..939aa64535 100644
--- a/metagraph/src/common/elias_fano/elias_fano_merger.hpp
+++ b/metagraph/src/common/elias_fano/elias_fano_merger.hpp
@@ -188,7 +188,8 @@ template <typename T>
 void merge_files(std::vector<std::string> sources,
                  const std::function<void(const T &)> &on_new_item,
                  bool remove_sources = true,
-                 size_t max_sources_open = -1) {
+                 size_t max_sources_open = -1,
+                 bool deduplicate = true) {
     if (!sources.size())
         return;
 
@@ -223,7 +224,7 @@ void merge_files(std::vector<std::string> sources,
                     EliasFanoEncoderBuffered<T>::append_block(buf, new_chunks.back());
                     buf.resize(0);
                 }
-            });
+            }, remove_sources, -1, deduplicate);
             if (buf.size()) {
                 EliasFanoEncoderBuffered<T>::append_block(buf, new_chunks.back());
             }
@@ -238,7 +239,7 @@ void merge_files(std::vector<std::string> sources,
     T last = decoder.pop();
     while (!decoder.empty()) {
         T curr = decoder.pop();
-        if (curr != last) {
+        if (!deduplicate || curr != last) {
             on_new_item(last);
             last = curr;
         }
diff --git a/metagraph/src/graph/representation/succinct/boss.cpp b/metagraph/src/graph/representation/succinct/boss.cpp
index dd112cf390..58ffb017e2 100644
--- a/metagraph/src/graph/representation/succinct/boss.cpp
+++ b/metagraph/src/graph/representation/succinct/boss.cpp
@@ -2496,7 +2496,7 @@ void update_terminal_bits(size_t max_length,
  * if the row-diff successor in a fork has already been visited.
  */
 void traverse_rd_path(const BOSS &boss,
-                      const bit_vector &rd_succ,
+                      const std::function<bool(edge_index)> &is_successor,
                       edge_index edge,
                       size_t max_length,
                       sdsl::bit_vector *visited,
@@ -2524,7 +2524,7 @@ void traverse_rd_path(const BOSS &boss,
 
         path.push_back(edge);
 
-        edge = boss.row_diff_successor(edge, rd_succ);
+        edge = boss.row_diff_successor(edge, is_successor);
     }
 
     // mark terminal and near terminal nodes
@@ -2720,12 +2720,12 @@ void BOSS::call_sequences(Call<std::string&&, std::vector<edge_index>&&> callbac
 // Reach all k-mers that merge into anchor |edge| by following their diff paths.
 template <typename T>
 void traverse_rd_path_backward(const BOSS &boss,
-                               const bit_vector &rd_succ,
+                               const std::function<bool(edge_index)> &is_successor,
                                edge_index edge,
                                T max_length,
                                sdsl::bit_vector *visited,
                                sdsl::bit_vector *terminal,
-                               sdsl::bit_vector *dummy,
+                               const sdsl::bit_vector &dummy,
                                ProgressBar &progress_bar) {
     constexpr bool async = true;
 
@@ -2740,7 +2740,7 @@ void traverse_rd_path_backward(const BOSS &boss,
         std::tie(edge, dist_to_anchor) = queue.back();
         queue.pop_back();
 
-        assert(boss.get_W(edge) && !fetch_bit(dummy->data(), edge, async));
+        assert(boss.get_W(edge) && !dummy[edge]);
 
         // mark as visited
         // also check if the node had already been visited (in case of loop)
@@ -2761,7 +2761,7 @@ void traverse_rd_path_backward(const BOSS &boss,
         // AAAX - AAX$
         // ^^^^
         // AAAY - ****
-        if (!rd_succ[edge])
+        if (!is_successor(edge))
             continue;
 
         // |edge| is the row-diff successor. Thus, it is part of a diff
@@ -2769,7 +2769,7 @@ void traverse_rd_path_backward(const BOSS &boss,
         // So, we propagate the diff path backward.
         boss.call_incoming_to_target(boss.bwd(edge), boss.get_node_last_value(edge),
             [&](edge_index pred) {
-                if (!fetch_bit(dummy->data(), pred, async))
+                if (!dummy[pred])
                     queue.emplace_back(pred, dist_to_anchor + 1);
             }
         );
@@ -2778,7 +2778,7 @@ void traverse_rd_path_backward(const BOSS &boss,
 
 void BOSS::row_diff_traverse(size_t num_threads,
                              size_t max_length,
-                             const bit_vector &rd_succ,
+                             const std::function<bool(edge_index)> &is_successor,
                              sdsl::bit_vector *terminal) const {
     // TODO: can we do it without using this extra |dummy| bitmap?
     sdsl::bit_vector dummy(W_->size(), false);
@@ -2848,8 +2848,8 @@ void BOSS::row_diff_traverse(size_t num_threads,
 
     // backward traversal
     traverse_path = [&](edge_index anchor) {
-        traverse_rd_path_backward(*this, rd_succ, anchor, max_length, &visited,
-                                  terminal, &dummy, progress_bar);
+        traverse_rd_path_backward(*this, is_successor, anchor, max_length, &visited,
+                                  terminal, dummy, progress_bar);
     };
 
     // run backward traversal from every anchor
@@ -2860,7 +2860,7 @@ void BOSS::row_diff_traverse(size_t num_threads,
 
     // forward traversal
     traverse_path = [&](edge_index start) {
-        traverse_rd_path(*this, rd_succ, start, max_length, &visited,
+        traverse_rd_path(*this, is_successor, start, max_length, &visited,
                          terminal, &near_terminal, progress_bar);
     };
     // start traversal from the dummy source edges first ($X...X)
diff --git a/metagraph/src/graph/representation/succinct/boss.hpp b/metagraph/src/graph/representation/succinct/boss.hpp
index 6b2ec5fe25..eed3fc9215 100644
--- a/metagraph/src/graph/representation/succinct/boss.hpp
+++ b/metagraph/src/graph/representation/succinct/boss.hpp
@@ -153,17 +153,18 @@ class BOSS {
      */
     void row_diff_traverse(size_t num_threads,
                            size_t max_length,
-                           const bit_vector &rd_succ,
+                           const std::function<bool(edge_index)> &is_successor,
                            sdsl::bit_vector *terminal) const;
 
-    edge_index row_diff_successor(edge_index edge, const bit_vector &rd_succ) const {
+    // TODO: move all row-diff related things to DBGSuccinct and remove this function
+    edge_index row_diff_successor(edge_index edge, const std::function<bool(edge_index)> &is_successor) const {
         TAlphabet d = get_W(edge) % alph_size;
         assert(d != kSentinelCode && "sinks have no row-diff successors");
         // make one traversal step
         edge = fwd(edge, d);
         // pick the row-diff successor
         if (!get_last(edge - 1)) {
-            while (!rd_succ[edge]) {
+            while (!is_successor(edge)) {
                 edge--;
                 assert(!get_last(edge) && "a row-diff successor must exist");
             }
diff --git a/metagraph/src/graph/representation/succinct/dbg_succinct.cpp b/metagraph/src/graph/representation/succinct/dbg_succinct.cpp
index 3f401b80be..fd9025c19a 100644
--- a/metagraph/src/graph/representation/succinct/dbg_succinct.cpp
+++ b/metagraph/src/graph/representation/succinct/dbg_succinct.cpp
@@ -6,6 +6,8 @@
 #include <string>
 #include <filesystem>
 
+#include <progress_bar.hpp>
+
 #include "common/seq_tools/reverse_complement.hpp"
 #include "common/serialization.hpp"
 #include "common/logger.hpp"
@@ -1030,5 +1032,62 @@ void DBGSuccinct::print(std::ostream &out) const {
     }
 }
 
+void DBGSuccinct::add_rd_successors_at_forks(size_t num_threads,
+                                             const sdsl::bit_vector &anchor,
+                                             sdsl::bit_vector *rd_succ,
+                                             size_t max_length) const {
+    const auto &W = boss_graph_->get_W();
+
+    ProgressBar progress_bar(W.size(), "Adding fork successors",
+                             std::cerr, !common::get_verbose());
+
+    constexpr bool async = true;
+    sdsl::bit_vector excluded(W.size(), false);
+
+    // start from 0 and go with blocks of 1024 bits to avoid race conditions
+    #pragma omp parallel for num_threads(num_threads) schedule(static, 1024)
+    for (BOSS::edge_index i = 0; i < W.size(); ++i) {
+        ++progress_bar;
+        // skip edge if it's not a fork or it's already selected
+        if (i < 2 || boss_graph_->is_single_outgoing(i) || (*rd_succ)[boss_to_kmer_index(i)])
+            continue;
+
+        // // make nodes with multiple incoming edges have only one successor
+        // if (indegree(boss_to_kmer_index(i)) > 1)
+        //     continue;
+
+        // (*rd_succ)[i] = true;
+        // continue;
+        // TODO: test this
+        // if (!is_single_incoming(i, get_W(i)))
+        //     continue;
+
+        std::vector<BOSS::edge_index> queue = { i };
+
+        // make branching edge a successor if it reaches anchors not too deep
+        for (size_t depth = 0; depth < max_length && queue.size(); ++depth) {
+            BOSS::edge_index edge = queue.back();
+            queue.pop_back();
+            // stop branch if anchor is reached
+            if (anchor[edge])
+                continue;
+
+            assert(W[edge] % boss_graph_->alph_size);
+
+            boss_graph_->call_outgoing(boss_graph_->fwd(edge, W[edge] % boss_graph_->alph_size),
+                [&](BOSS::edge_index next) {
+                    if (!fetch_bit(excluded.data(), next, async))
+                        queue.push_back(next);
+                }
+            );
+        }
+        if (queue.empty()) {
+            (*rd_succ)[boss_to_kmer_index(i)] = true;
+        } else {
+            set_bit(excluded.data(), i, async);
+        }
+    }
+}
+
 } // namespace graph
 } // namespace mtg
diff --git a/metagraph/src/graph/representation/succinct/dbg_succinct.hpp b/metagraph/src/graph/representation/succinct/dbg_succinct.hpp
index 698678504c..f80ecffdf0 100644
--- a/metagraph/src/graph/representation/succinct/dbg_succinct.hpp
+++ b/metagraph/src/graph/representation/succinct/dbg_succinct.hpp
@@ -172,6 +172,15 @@ class DBGSuccinct : public DeBruijnGraph {
 
     virtual void call_source_nodes(const std::function<void(node_index)> &callback) const override final;
 
+    node_index row_diff_successor(node_index node, const bit_vector &rd_succ) const;
+    template <class Callback>
+    void call_row_diff_successors(node_index node, const bit_vector &rd_succ, const Callback &callback) const;
+
+    void add_rd_successors_at_forks(size_t num_threads,
+                                    const sdsl::bit_vector &anchors,
+                                    sdsl::bit_vector *rd_succ,
+                                    size_t max_length) const;
+
     uint64_t kmer_to_boss_index(node_index kmer_index) const;
     node_index boss_to_kmer_index(uint64_t boss_index) const;
 
@@ -197,6 +206,53 @@ class DBGSuccinct : public DeBruijnGraph {
     std::unique_ptr<mtg::kmer::KmerBloomFilter<>> bloom_filter_;
 };
 
+inline DBGSuccinct::node_index
+DBGSuccinct::row_diff_successor(node_index node, const bit_vector &rd_succ) const {
+    const boss::BOSS &boss = *boss_graph_;
+    boss::BOSS::edge_index edge = kmer_to_boss_index(node);
+    boss::BOSS::TAlphabet d = boss.get_W(edge) % boss.alph_size;
+    assert(d != boss::BOSS::kSentinelCode && "sinks have no row-diff successors");
+    // make one traversal step
+    edge = boss.fwd(edge, d);
+    node = boss_to_kmer_index(edge);
+
+    if (!rd_succ.size() || boss.get_last(edge - 1))
+        return node;
+
+    // pick the row-diff successor
+    while (!rd_succ[node]) {
+        node--;
+        edge--;
+        assert(!boss.get_last(edge) && "a row-diff successor must exist");
+    }
+    return node;
+}
+
+template <class Callback>
+inline void DBGSuccinct::call_row_diff_successors(node_index node,
+                                                  const bit_vector &rd_succ,
+                                                  const Callback &callback) const {
+    const boss::BOSS &boss = *boss_graph_;
+    boss::BOSS::edge_index edge = kmer_to_boss_index(node);
+    boss::BOSS::TAlphabet d = boss.get_W(edge) % boss.alph_size;
+    assert(d != boss::BOSS::kSentinelCode && "sinks have no row-diff successors");
+    // make one traversal step
+    edge = boss.fwd(edge, d);
+    node = boss_to_kmer_index(edge);
+
+    if (!rd_succ.size() || boss.get_last(edge - 1)) {
+        callback(node);
+        return;
+    }
+
+    // pick the row-diff successor
+    do {
+        if (rd_succ[node])
+            callback(node);
+        node--;
+    } while (!boss.get_last(--edge));
+}
+
 } // namespace graph
 } // namespace mtg
 
diff --git a/metagraph/tests/graph/succinct/test_boss.cpp b/metagraph/tests/graph/succinct/test_boss.cpp
index 39da899143..23d42ab266 100644
--- a/metagraph/tests/graph/succinct/test_boss.cpp
+++ b/metagraph/tests/graph/succinct/test_boss.cpp
@@ -952,7 +952,7 @@ TEST(BOSS, CallSequencesRowDiff_EmptyGraph) {
             BOSS empty(k);
 
             sdsl::bit_vector terminal(empty.get_last().size(), false);
-            empty.row_diff_traverse(num_threads, 1, empty.get_last(), &terminal);
+            empty.row_diff_traverse(num_threads, 1, [&](auto i) { return empty.get_last(i); }, &terminal);
             ASSERT_EQ(sdsl::bit_vector(2, false), terminal)
                 << "Empty graph must have no anchors";
         }
@@ -1067,7 +1067,7 @@ TEST(BOSS, CallSequenceRowDiff_TwoLoops) {
             ASSERT_EQ(2u, graph.num_edges());
 
             sdsl::bit_vector terminal(graph.get_last().size(), false);
-            graph.row_diff_traverse(num_threads, 1, graph.get_last(), &terminal);
+            graph.row_diff_traverse(num_threads, 1, [&](auto i) { return graph.get_last(i); }, &terminal);
             ASSERT_EQ(graph.num_edges() + 1, terminal.size());
             ASSERT_EQ(sdsl::bit_vector({ 0, 0, 1 }), terminal);
         }
@@ -1128,11 +1128,11 @@ TEST(BOSS, CallSequenceRowDiff_TwoBigLoops) {
         BOSS graph(&constructor);
 
         sdsl::bit_vector terminal(graph.get_last().size(), false);
-        graph.row_diff_traverse(num_threads, 100, graph.get_last(), &terminal);
+        graph.row_diff_traverse(num_threads, 100, [&](auto i) { return graph.get_last(i); }, &terminal);
         ASSERT_EQ(graph.num_edges() + 1, terminal.size());
         ASSERT_EQ(2, std::accumulate(terminal.begin() + 1, terminal.end(), 0U));
 
-        graph.row_diff_traverse(num_threads, 1, graph.get_last(), &terminal);
+        graph.row_diff_traverse(num_threads, 1, [&](auto i) { return graph.get_last(i); }, &terminal);
         ASSERT_EQ(graph.num_edges() + 1, terminal.size());
         ASSERT_EQ(104, std::accumulate(terminal.begin() + 1, terminal.end(), 0U));
     }
@@ -1204,7 +1204,7 @@ TEST(BOSS, CallSequenceRowDiff_FourLoops) {
             BOSS graph(&constructor);
 
             sdsl::bit_vector terminal(graph.get_last().size(), false);
-            graph.row_diff_traverse(num_threads, 1, graph.get_last(), &terminal);
+            graph.row_diff_traverse(num_threads, 1, [&](auto i) { return graph.get_last(i); }, &terminal);
             ASSERT_EQ(graph.num_edges() + 1, terminal.size());
             ASSERT_EQ(4, std::accumulate(terminal.begin() + 1, terminal.end(), 0U));
         }
@@ -1220,11 +1220,11 @@ TEST(BOSS, CallSequenceRowDiff_FourPaths) {
     BOSS graph(&constructor);
     for (size_t num_threads : { 1, 4 }) {
         sdsl::bit_vector terminal(graph.get_last().size(), false);
-        graph.row_diff_traverse(num_threads, 20, graph.get_last(), &terminal);
+        graph.row_diff_traverse(num_threads, 20, [&](auto i) { return graph.get_last(i); }, &terminal);
         ASSERT_EQ(graph.num_edges() + 1, terminal.size());
         ASSERT_EQ(4, std::accumulate(terminal.begin() + 1, terminal.end(), 0U));
 
-        graph.row_diff_traverse(num_threads, 1, graph.get_last(), &terminal);
+        graph.row_diff_traverse(num_threads, 1, [&](auto i) { return graph.get_last(i); }, &terminal);
         ASSERT_EQ(graph.num_edges() + 1, terminal.size());
         ASSERT_EQ(19, std::accumulate(terminal.begin() + 1, terminal.end(), 0U));
     }