diff --git a/src/sage/algebras/fusion_rings/f_matrix.py b/src/sage/algebras/fusion_rings/f_matrix.py
index a8dc22901c8..e178b567a3e 100644
--- a/src/sage/algebras/fusion_rings/f_matrix.py
+++ b/src/sage/algebras/fusion_rings/f_matrix.py
@@ -1726,11 +1726,11 @@ def _partition_eqns(self, eqns=None, verbose=True):
         if eqns is None:
             eqns = self.ideal_basis
         graph = self.equations_graph(eqns)
-        partition = {tuple(c): [] for c in graph.connected_components()}
+        partition = {tuple(c): [] for c in graph.connected_components(sort=True)}
         for eq_tup in eqns:
-            partition[tuple(graph.connected_component_containing_vertex(variables(eq_tup)[0]))].append(eq_tup)
+            partition[tuple(graph.connected_component_containing_vertex(variables(eq_tup)[0], sort=True))].append(eq_tup)
         if verbose:
-            print("Partitioned {} equations into {} components of size:".format(len(eqns), len(graph.connected_components())))
+            print("Partitioned {} equations into {} components of size:".format(len(eqns), graph.connected_components_number()))
             print(graph.connected_components_sizes())
         return partition
 
diff --git a/src/sage/categories/complex_reflection_or_generalized_coxeter_groups.py b/src/sage/categories/complex_reflection_or_generalized_coxeter_groups.py
index 3bbe54cefc0..3dff3b983a6 100644
--- a/src/sage/categories/complex_reflection_or_generalized_coxeter_groups.py
+++ b/src/sage/categories/complex_reflection_or_generalized_coxeter_groups.py
@@ -766,7 +766,7 @@ def irreducible_component_index_sets(self):
                         for i,j in itertools.combinations(I,2)
                         if s[i]*s[j] != s[j]*s[i] ]],
                       format="vertices_and_edges")
-            return G.connected_components()
+            return G.connected_components(sort=False)
 
         @abstract_method(optional=True)
         def irreducible_components(self):
diff --git a/src/sage/categories/loop_crystals.py b/src/sage/categories/loop_crystals.py
index 2d5f1ca1914..83c403323cd 100644
--- a/src/sage/categories/loop_crystals.py
+++ b/src/sage/categories/loop_crystals.py
@@ -868,7 +868,7 @@ def one_dimensional_configuration_sum(self, q=None, group_components=True):
                 B = P0.algebra(q.parent())
                 if group_components:
                     G = self.digraph(index_set=self.cartan_type().classical().index_set())
-                    C = G.connected_components()
+                    C = G.connected_components(sort=False)
                     return B.sum(q**(c[0].energy_function())*B.sum(B(P0(b.weight())) for b in c)
                                  for c in C)
                 return B.sum(q**(b.energy_function())*B(P0(b.weight())) for b in self)
diff --git a/src/sage/combinat/cluster_algebra_quiver/mutation_type.py b/src/sage/combinat/cluster_algebra_quiver/mutation_type.py
index fc414b8f8fb..8707841b137 100644
--- a/src/sage/combinat/cluster_algebra_quiver/mutation_type.py
+++ b/src/sage/combinat/cluster_algebra_quiver/mutation_type.py
@@ -896,7 +896,7 @@ def _connected_mutation_type_AAtildeD(dg, ret_conn_vert=False):
                     dg_tmp = DiGraph( dg )
                     dg_tmp.delete_vertices( c1 )
 
-                    components = dg_tmp.connected_components()
+                    components = dg_tmp.connected_components(sort=False)
                     # if not len(components) == 2:
                     if len(components) != 2:
                         return _false_return(4)
@@ -938,7 +938,7 @@ def _connected_mutation_type_AAtildeD(dg, ret_conn_vert=False):
                     else:
                         c2.reverse()
                         dg_tmp.delete_edge( tuple( c2 ) )
-                    components = dg_tmp.connected_components()
+                    components = dg_tmp.connected_components(sort=False)
                     if len(components) != 2:
                         return _false_return(7)
                     else:
@@ -1190,7 +1190,7 @@ def _connected_mutation_type_AAtildeD(dg, ret_conn_vert=False):
             edge = long_cycle[0][0]
             sg = DiGraph( dg )
             sg. delete_vertices(edge)
-            connected_components = sg.connected_components()
+            connected_components = sg.connected_components(sort=False)
             cycle = []
             if connected_components:
                 cycle.append( ( edge[0], edge[1], len( connected_components[0] ) + 1 ) )
@@ -1200,7 +1200,7 @@ def _connected_mutation_type_AAtildeD(dg, ret_conn_vert=False):
             for edge in tmp:
                 sg = DiGraph( dg )
                 sg. delete_vertices(edge)
-                connected_components = sg.connected_components()
+                connected_components = sg.connected_components(sort=False)
                 if len( connected_components ) == 2:
                     #if len( list_intersection( [ connected_components[0], list_substract( long_cycle[0], [edge] )[0] ] ) ) > 0:
                     if len( set(connected_components[0]).intersection( set(long_cycle[0]).difference([edge]).pop() ) ) > 0:
diff --git a/src/sage/combinat/cluster_algebra_quiver/quiver.py b/src/sage/combinat/cluster_algebra_quiver/quiver.py
index 7b2ad09b1b2..dff95dcef24 100644
--- a/src/sage/combinat/cluster_algebra_quiver/quiver.py
+++ b/src/sage/combinat/cluster_algebra_quiver/quiver.py
@@ -1007,7 +1007,7 @@ def mutation_type(self):
 
             # checking the type for each connected component
             mutation_type = []
-            connected_components = sorted(dg.connected_components())
+            connected_components = sorted(dg.connected_components(sort=False))
             for component in connected_components:
                 # constructing the digraph for this component
                 dg_component = dg.subgraph( component )
@@ -1157,7 +1157,7 @@ def canonical_label(self, certificate=False):
             if dg.is_connected():
                 Q._mutation_type = self._mutation_type
             else:
-                CC = sorted( self._digraph.connected_components() )
+                CC = sorted(self._digraph.connected_components(sort=False))
                 CC_new = sorted(zip([sorted(iso[i] for i in L) for L in CC],
                                     range(len(CC))))
                 comp_iso = [L[1] for L in CC_new]
diff --git a/src/sage/combinat/constellation.py b/src/sage/combinat/constellation.py
index 0f5350204aa..8f7274cc256 100644
--- a/src/sage/combinat/constellation.py
+++ b/src/sage/combinat/constellation.py
@@ -475,7 +475,7 @@ def connected_components(self):
         G.add_vertices(list(range(self.degree())))
         for p in self._g:
             G.add_edges(enumerate(p.domain()), loops=False)
-        m = G.connected_components()
+        m = G.connected_components(sort=False)
         if len(m) == 1:
             return [self]
         for mm in m:
diff --git a/src/sage/combinat/crystals/littelmann_path.py b/src/sage/combinat/crystals/littelmann_path.py
index 1af0fed2dcb..c384a818e56 100644
--- a/src/sage/combinat/crystals/littelmann_path.py
+++ b/src/sage/combinat/crystals/littelmann_path.py
@@ -809,7 +809,7 @@ def weight(x):
             return P0.sum(int(c)*P0.basis()[i] for i,c in w if i in P0.index_set())
         if group_components:
             G = self.digraph(index_set=self.cartan_type().classical().index_set())
-            C = G.connected_components()
+            C = G.connected_components(sort=False)
             return sum(q**(c[0].energy_function())*B.sum(B(weight(b)) for b in c) for c in C)
         return B.sum(q**(b.energy_function())*B(weight(b)) for b in self)
 
@@ -1081,7 +1081,7 @@ def energy_function(self):
                 sage: La = R.weight_space().basis()
                 sage: LS = crystals.ProjectedLevelZeroLSPaths(2*La[1]+La[2])
                 sage: G = LS.digraph(index_set=[1,2])
-                sage: C = G.connected_components()
+                sage: C = G.connected_components(sort=False)
                 sage: [all(c[0].energy_function()==a.energy_function() for a in c) for c in C]
                 [True, True, True, True]
 
@@ -1093,7 +1093,7 @@ def energy_function(self):
                 sage: [(x.weight(), x.energy_function()) for x in hw]
                 [(-2*Lambda[0] + Lambda[2], 0), (-2*Lambda[0] + Lambda[1], 1), (0, 2)]
                 sage: G = LS.digraph(index_set=J)
-                sage: C = G.connected_components()
+                sage: C = G.connected_components(sort=False)
                 sage: [all(c[0].energy_function()==a.energy_function() for a in c) for c in C]
                 [True, True, True]
 
@@ -1101,7 +1101,7 @@ def energy_function(self):
                 sage: La = R.weight_space().basis()
                 sage: LS = crystals.ProjectedLevelZeroLSPaths(La[1]+La[2])
                 sage: G = LS.digraph(index_set=[1,2])
-                sage: C = G.connected_components()
+                sage: C = G.connected_components(sort=False)
                 sage: [all(c[0].energy_function()==a.energy_function() for a in c) for c in C] # long time
                 [True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True]
 
@@ -1110,7 +1110,7 @@ def energy_function(self):
                 sage: La = R.weight_space().basis()
                 sage: LS = crystals.ProjectedLevelZeroLSPaths(2*La[1]+La[2])
                 sage: G = LS.digraph(index_set=R.cartan_type().classical().index_set())
-                sage: C = G.connected_components()
+                sage: C = G.connected_components(sort=False)
                 sage: [all(c[0].energy_function()==a.energy_function() for a in c) for c in C] # long time
                 [True, True, True, True, True, True, True, True, True, True, True]
 
@@ -1118,7 +1118,7 @@ def energy_function(self):
                 sage: La = R.weight_space().basis()
                 sage: LS = crystals.ProjectedLevelZeroLSPaths(2*La[1]+La[2])
                 sage: G = LS.digraph(index_set=R.cartan_type().classical().index_set())
-                sage: C = G.connected_components()
+                sage: C = G.connected_components(sort=False)
                 sage: [all(c[0].energy_function()==a.energy_function() for a in c) for c in C] # long time
                 [True, True, True, True, True, True, True, True, True, True, True, True, True, True, True,
                 True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True]
diff --git a/src/sage/combinat/partition_algebra.py b/src/sage/combinat/partition_algebra.py
index 9136b682c59..0b7eb06c3a5 100644
--- a/src/sage/combinat/partition_algebra.py
+++ b/src/sage/combinat/partition_algebra.py
@@ -1974,7 +1974,7 @@ def set_partition_composition(sp1, sp2):
         True
     """
     g = pair_to_graph(sp1, sp2)
-    connected_components = g.connected_components()
+    connected_components = g.connected_components(sort=False)
 
     res = []
     total_removed = 0
diff --git a/src/sage/combinat/posets/mobile.py b/src/sage/combinat/posets/mobile.py
index d09eb370fa7..f291400ce5f 100644
--- a/src/sage/combinat/posets/mobile.py
+++ b/src/sage/combinat/posets/mobile.py
@@ -123,7 +123,7 @@ def _is_valid_ribbon(self, ribbon):
                     continue
 
                 G_un.delete_edge(lc, r)
-                P = Poset(G.subgraph(G_un.connected_component_containing_vertex(lc)))
+                P = Poset(G.subgraph(G_un.connected_component_containing_vertex(lc, sort=False)))
                 if P.top() != lc or not P.is_d_complete():
                     return False
                 G_un.add_edge(lc, r)
diff --git a/src/sage/combinat/posets/posets.py b/src/sage/combinat/posets/posets.py
index 65ad12be6bf..d37e9676d7c 100644
--- a/src/sage/combinat/posets/posets.py
+++ b/src/sage/combinat/posets/posets.py
@@ -1368,11 +1368,11 @@ def _latex_(self):
         EXAMPLES::
 
             sage: P = Poset(([1,2], [[1,2]]), cover_relations = True)
-            sage: print(P._latex_()) #optional - dot2tex graphviz
+            sage: print(P._latex_())  # optional - dot2tex graphviz
             \begin{tikzpicture}[>=latex,line join=bevel,]
             %%
-            \node (node_...) at (6.0...bp,...bp) [draw,draw=none] {$...$};
-              \node (node_...) at (6.0...bp,...bp) [draw,draw=none] {$...$};
+            \node (node_...) at (5...bp,...bp) [draw,draw=none] {$...$};
+              \node (node_...) at (5...bp,...bp) [draw,draw=none] {$...$};
               \draw [black,->] (node_...) ..controls (...bp,...bp) and (...bp,...bp)  .. (node_...);
             %
             \end{tikzpicture}
@@ -5330,7 +5330,7 @@ def edge_color(va, vb):
 
         fusion = fusion.transitive_closure()
         resu = []
-        for s in fusion.connected_components():
+        for s in fusion.connected_components(sort=False):
             subg = [x for x in prod_dg if all(x[i] == v0[i] for i in factors_range
                                               if i not in s)]
             resu.append(Poset(prod_dg.subgraph(subg)))
diff --git a/src/sage/combinat/root_system/ambient_space.py b/src/sage/combinat/root_system/ambient_space.py
index b5e4a90e6ad..7c7159a444c 100644
--- a/src/sage/combinat/root_system/ambient_space.py
+++ b/src/sage/combinat/root_system/ambient_space.py
@@ -119,7 +119,7 @@ def _test_norm_of_simple_roots(self, **options):
         T = self.cartan_type()
         D = T.symmetrizer()
         alpha = self.simple_roots()
-        for C in T.dynkin_diagram().connected_components():
+        for C in T.dynkin_diagram().connected_components(sort=False):
             tester.assertEqual(len( set( alpha[i].scalar(alpha[i]) / D[i] for i in C ) ), 1)
 
     # FIXME: attribute or method?
diff --git a/src/sage/combinat/root_system/cartan_type.py b/src/sage/combinat/root_system/cartan_type.py
index d8a37f350ce..013ceca8a02 100644
--- a/src/sage/combinat/root_system/cartan_type.py
+++ b/src/sage/combinat/root_system/cartan_type.py
@@ -1739,7 +1739,7 @@ def symmetrizer(self):
             M[i, n * i + j] = m[i,j]
             M[j, n * i + j] -= m[j,i]
         kern = M.integer_kernel()
-        c = len(self.dynkin_diagram().connected_components())
+        c = len(self.dynkin_diagram().connected_components(sort=False))
         if kern.dimension() < c:
             # the Cartan matrix is not symmetrizable
             return None
diff --git a/src/sage/geometry/polyhedral_complex.py b/src/sage/geometry/polyhedral_complex.py
index cff7b73c94e..bec18228aca 100644
--- a/src/sage/geometry/polyhedral_complex.py
+++ b/src/sage/geometry/polyhedral_complex.py
@@ -1234,7 +1234,7 @@ def connected_component(self, cell=None):
             if v not in g:
                 raise ValueError(
                     "the polyhedral complex does not contain the given cell")
-            vertices = g.connected_component_containing_vertex(v)
+            vertices = g.connected_component_containing_vertex(v, sort=False)
             facets = [f for f in self.maximal_cell_iterator()
                       if any(vf in f.vertices_matrix().columns()
                              for vf in vertices)]
@@ -1243,7 +1243,7 @@ def connected_component(self, cell=None):
             if cell not in g:
                 raise ValueError(
                     "the polyhedral complex does not contain the given cell")
-            faces = g.connected_component_containing_vertex(cell)
+            faces = g.connected_component_containing_vertex(cell, sort=False)
             facets = [f for f in self.maximal_cell_iterator()
                       if f in faces]
         return PolyhedralComplex(facets, maximality_check=False,
diff --git a/src/sage/geometry/polyhedron/base_QQ.py b/src/sage/geometry/polyhedron/base_QQ.py
index d279aa5ea7f..7aa1c568480 100644
--- a/src/sage/geometry/polyhedron/base_QQ.py
+++ b/src/sage/geometry/polyhedron/base_QQ.py
@@ -1180,7 +1180,7 @@ class functions of the acting group. A character `\rho` is effective if
             True
             sage: Hstar = p3.Hstar_function(S3)                         # optional - pynormaliz
             sage: Hlin  = p3.Hstar_function(S3,                         # optional - pynormaliz
-            ....:                           output='Hstar_as_lin_comb')]
+            ....:                           output='Hstar_as_lin_comb')
             sage: p3.is_effective(Hstar, Hlin)                          # optional - pynormaliz
             True
 
@@ -1242,7 +1242,7 @@ class functions of the acting group. A character `\rho` is effective if
             sage: p2 = Polyhedron(vertices=[[0], [1/2]],                     # optional - pynormaliz
             ....:                 backend='normaliz')
             sage: Hstar    = p2.Hstar_function()                             # optional - pynormaliz
-            sage: Hstarlin = p2.Hstar_function(output='Hstar_as_lin_comb')]  # optional - pynormaliz
+            sage: Hstarlin = p2.Hstar_function(output='Hstar_as_lin_comb')   # optional - pynormaliz
             sage: p1._is_effective_normaliz(Hstar, Hstarlin)                 # optional - pynormaliz
             Traceback (most recent call last):
             ...
diff --git a/src/sage/graphs/base/c_graph.pyx b/src/sage/graphs/base/c_graph.pyx
index d6fb119d3cf..a4f2efd6060 100644
--- a/src/sage/graphs/base/c_graph.pyx
+++ b/src/sage/graphs/base/c_graph.pyx
@@ -4819,7 +4819,7 @@ cdef class Search_iterator:
 
         Immutable graphs (see :trac:`16019`)::
 
-            sage: DiGraph([[1,2]], immutable=True).connected_components()
+            sage: DiGraph([(1, 2)], immutable=True).connected_components(sort=True)
             [[1, 2]]
 
         """
diff --git a/src/sage/graphs/base/static_sparse_graph.pyx b/src/sage/graphs/base/static_sparse_graph.pyx
index a1f246883ff..cfaff95ccde 100644
--- a/src/sage/graphs/base/static_sparse_graph.pyx
+++ b/src/sage/graphs/base/static_sparse_graph.pyx
@@ -702,7 +702,7 @@ def tarjan_strongly_connected_components(G):
         sage: tarjan_strongly_connected_components(digraphs.Path(3))
         [[2], [1], [0]]
         sage: D = DiGraph( { 0 : [1, 3], 1 : [2], 2 : [3], 4 : [5, 6], 5 : [6] } )
-        sage: D.connected_components()
+        sage: D.connected_components(sort=True)
         [[0, 1, 2, 3], [4, 5, 6]]
         sage: D = DiGraph( { 0 : [1, 3], 1 : [2], 2 : [3], 4 : [5, 6], 5 : [6] } )
         sage: D.strongly_connected_components()
diff --git a/src/sage/graphs/bipartite_graph.py b/src/sage/graphs/bipartite_graph.py
index 86b6f4bba15..635d5554329 100644
--- a/src/sage/graphs/bipartite_graph.py
+++ b/src/sage/graphs/bipartite_graph.py
@@ -1209,7 +1209,7 @@ def _check_bipartition_for_add_edges(self, edges):
             vertex_in_left[v] = False
 
         # Map each vertex to the connected component it belongs to
-        vertex_to_component = {v: comp for comp in self.connected_components()
+        vertex_to_component = {v: comp for comp in self.connected_components(sort=False)
                                    for v in comp}
 
         for e in edges:
diff --git a/src/sage/graphs/comparability.pyx b/src/sage/graphs/comparability.pyx
index 0271f43d9d8..b694b552ee8 100644
--- a/src/sage/graphs/comparability.pyx
+++ b/src/sage/graphs/comparability.pyx
@@ -257,7 +257,7 @@ def greedy_is_comparability(g, no_certificate=False, equivalence_class=False):
     # Each vertex can partition its neighbors into equivalence classes
     equivalence_classes = {}
     for v in g:
-        equivalence_classes[v] = g.subgraph(vertices=g.neighbors(v)).complement().connected_components()
+        equivalence_classes[v] = g.subgraph(vertices=g.neighbors(v)).complement().connected_components(sort=False)
 
     # We build a graph h with one vertex per (vertex of g + equivalence class)
     from sage.graphs.graph import Graph
@@ -288,7 +288,7 @@ def greedy_is_comparability(g, no_certificate=False, equivalence_class=False):
         if equivalence_class:
 
             # Returning the largest equivalence class
-            cc = sorted(h.connected_components(), key=len)[-1]
+            cc = sorted(h.connected_components(sort=False), key=len)[-1]
 
             edges = []
             for v, sid in cc:
diff --git a/src/sage/graphs/connectivity.pyx b/src/sage/graphs/connectivity.pyx
index 0cbcb87dd8d..c2100141f30 100644
--- a/src/sage/graphs/connectivity.pyx
+++ b/src/sage/graphs/connectivity.pyx
@@ -58,6 +58,18 @@ Methods
 -------
 """
 
+# ****************************************************************************
+#
+#       Copyright (C) 2023      David Coudert <david.coudert@inria.fr>
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 2 of the License, or
+# (at your option) any later version.
+#                  https://www.gnu.org/licenses/
+# ****************************************************************************
+
+from sage.misc.superseded import deprecation
 
 def is_connected(G):
     """
@@ -119,7 +131,7 @@ def is_connected(G):
         return len(conn_verts) == G.num_verts()
 
 
-def connected_components(G, sort=True):
+def connected_components(G, sort=None, key=None):
     """
     Return the list of connected components.
 
@@ -130,20 +142,30 @@ def connected_components(G, sort=True):
 
     - ``G`` -- the input graph
 
-    - ``sort`` -- boolean (default ``True``); whether to sort vertices inside
-      each component
+    - ``sort`` -- boolean (default: ``None``); if ``True``, vertices inside each
+      component are sorted according to the default ordering
+
+      As of :trac:`35889`, this argument must be explicitly specified (unless a
+      ``key`` is given); otherwise a warning is printed and ``sort=True`` is
+      used. The default will eventually be changed to ``False``.
+
+    - ``key`` -- a function (default: ``None``); a function that takes a
+      vertex as its one argument and returns a value that can be used for
+      comparisons in the sorting algorithm (we must have ``sort=True``)
 
     EXAMPLES::
 
         sage: from sage.graphs.connectivity import connected_components
         sage: G = Graph({0: [1, 3], 1: [2], 2: [3], 4: [5, 6], 5: [6]})
-        sage: connected_components(G)
+        sage: connected_components(G, sort=True)
         [[0, 1, 2, 3], [4, 5, 6]]
-        sage: G.connected_components()
+        sage: G.connected_components(sort=True)
         [[0, 1, 2, 3], [4, 5, 6]]
         sage: D = DiGraph({0: [1, 3], 1: [2], 2: [3], 4: [5, 6], 5: [6]})
-        sage: connected_components(D)
+        sage: connected_components(D, sort=True)
         [[0, 1, 2, 3], [4, 5, 6]]
+        sage: connected_components(D, sort=True, key=lambda x: -x)
+        [[3, 2, 1, 0], [6, 5, 4]]
 
     TESTS:
 
@@ -154,16 +176,40 @@ def connected_components(G, sort=True):
         Traceback (most recent call last):
         ...
         TypeError: the input must be a Sage graph
+
+    When parameter ``key`` is set, parameter ``sort`` must be ``True``::
+
+        sage: G = Graph(2)
+        sage: G.connected_components(sort=False, key=lambda x: x)
+        Traceback (most recent call last):
+        ...
+        ValueError: sort keyword is False, yet a key function is given
+
+    Deprecation warning for ``sort=None`` (:trac:`35889`)::
+
+        sage: G = graphs.HouseGraph()
+        sage: G.connected_components()
+        doctest:...: DeprecationWarning: parameter 'sort' will be set to False by default in the future
+        See https://github.com/sagemath/sage/issues/35889 for details.
+        [[0, 1, 2, 3, 4]]
     """
     from sage.graphs.generic_graph import GenericGraph
     if not isinstance(G, GenericGraph):
         raise TypeError("the input must be a Sage graph")
 
+    if sort is None:
+        if key is None:
+            deprecation(35889, "parameter 'sort' will be set to False by default in the future")
+        sort = True
+
+    if (not sort) and key:
+        raise ValueError('sort keyword is False, yet a key function is given')
+
     cdef set seen = set()
     cdef list components = []
     for v in G:
         if v not in seen:
-            c = connected_component_containing_vertex(G, v, sort=sort)
+            c = connected_component_containing_vertex(G, v, sort=sort, key=key)
             seen.update(c)
             components.append(c)
     components.sort(key=lambda comp: -len(comp))
@@ -236,7 +282,7 @@ def connected_components_subgraphs(G):
     return [G.subgraph(c, inplace=False) for c in connected_components(G, sort=False)]
 
 
-def connected_component_containing_vertex(G, vertex, sort=True):
+def connected_component_containing_vertex(G, vertex, sort=None, key=None):
     """
     Return a list of the vertices connected to vertex.
 
@@ -246,20 +292,30 @@ def connected_component_containing_vertex(G, vertex, sort=True):
 
     - ``v`` -- the vertex to search for
 
-    - ``sort`` -- boolean (default ``True``); whether to sort vertices inside
-      the component
+    - ``sort`` -- boolean (default: ``None``); if ``True``, vertices inside the
+      component are sorted according to the default ordering
+
+      As of :trac:`35889`, this argument must be explicitly specified (unless a
+      ``key`` is given); otherwise a warning is printed and ``sort=True`` is
+      used. The default will eventually be changed to ``False``.
+
+    - ``key`` -- a function (default: ``None``); a function that takes a
+      vertex as its one argument and returns a value that can be used for
+      comparisons in the sorting algorithm (we must have ``sort=True``)
 
     EXAMPLES::
 
         sage: from sage.graphs.connectivity import connected_component_containing_vertex
         sage: G = Graph({0: [1, 3], 1: [2], 2: [3], 4: [5, 6], 5: [6]})
-        sage: connected_component_containing_vertex(G, 0)
+        sage: connected_component_containing_vertex(G, 0, sort=True)
         [0, 1, 2, 3]
-        sage: G.connected_component_containing_vertex(0)
+        sage: G.connected_component_containing_vertex(0, sort=True)
         [0, 1, 2, 3]
         sage: D = DiGraph({0: [1, 3], 1: [2], 2: [3], 4: [5, 6], 5: [6]})
-        sage: connected_component_containing_vertex(D, 0)
+        sage: connected_component_containing_vertex(D, 0, sort=True)
         [0, 1, 2, 3]
+        sage: connected_component_containing_vertex(D, 0, sort=True, key=lambda x: -x)
+        [3, 2, 1, 0]
 
     TESTS:
 
@@ -270,18 +326,52 @@ def connected_component_containing_vertex(G, vertex, sort=True):
         Traceback (most recent call last):
         ...
         TypeError: the input must be a Sage graph
+
+    :trac:`35889` is fixed::
+
+        sage: G = Graph([('A', 1)])
+        sage: G.connected_component_containing_vertex(1, sort=False)
+        [1, 'A']
+        sage: G.connected_component_containing_vertex(1, sort=True)
+        Traceback (most recent call last):
+        ...
+        TypeError: '<' not supported between instances of 'str' and 'int'
+
+    When parameter ``key`` is set, parameter ``sort`` must be ``True``::
+
+        sage: G = Graph(2)
+        sage: G.connected_component_containing_vertex(1, sort=False, key=lambda x: x)
+        Traceback (most recent call last):
+        ...
+        ValueError: sort keyword is False, yet a key function is given
+
+    Deprecation warning for ``sort=None`` (:trac:`35889`)::
+
+        sage: G = graphs.HouseGraph()
+        sage: G.connected_component_containing_vertex(1)
+        doctest:...: DeprecationWarning: parameter 'sort' will be set to False by default in the future
+        See https://github.com/sagemath/sage/issues/35889 for details.
+        [0, 1, 2, 3, 4]
     """
     from sage.graphs.generic_graph import GenericGraph
     if not isinstance(G, GenericGraph):
         raise TypeError("the input must be a Sage graph")
 
+    if sort is None:
+        if key is None:
+            deprecation(35889, "parameter 'sort' will be set to False by default in the future")
+        sort = True
+
+    if (not sort) and key:
+        raise ValueError('sort keyword is False, yet a key function is given')
+
     try:
         c = list(G._backend.depth_first_search(vertex, ignore_direction=True))
     except AttributeError:
         c = list(G.depth_first_search(vertex, ignore_direction=True))
 
     if sort:
-        c.sort()
+        return sorted(c, key=key)
     return c
 
 
@@ -325,7 +415,7 @@ def connected_components_sizes(G):
     return [len(cc) for cc in connected_components(G, sort=False)]
 
 
-def blocks_and_cut_vertices(G, algorithm="Tarjan_Boost", sort=False):
+def blocks_and_cut_vertices(G, algorithm="Tarjan_Boost", sort=False, key=None):
     """
     Return the blocks and cut vertices of the graph.
 
@@ -351,6 +441,10 @@ def blocks_and_cut_vertices(G, algorithm="Tarjan_Boost", sort=False):
       the components and the list of cut vertices
       **currently only available for ``"Tarjan_Sage"``**
 
+    - ``key`` -- a function (default: ``None``); a function that takes a
+      vertex as its one argument and returns a value that can be used for
+      comparisons in the sorting algorithm (we must have ``sort=True``)
+
     OUTPUT: ``(B, C)``, where ``B`` is a list of blocks - each is a list of
     vertices and the blocks are the corresponding induced subgraphs - and
     ``C`` is a list of cut vertices.
@@ -439,6 +533,9 @@ def blocks_and_cut_vertices(G, algorithm="Tarjan_Boost", sort=False):
         raise NotImplementedError("blocks and cut vertices algorithm '%s' is not implemented" % algorithm)
 
     # If algorithm is "Tarjan_Sage"
+    if (not sort) and key:
+        raise ValueError('sort keyword is False, yet a key function is given')
+
     blocks = []
     cut_vertices = set()
 
@@ -531,7 +628,7 @@ def blocks_and_cut_vertices(G, algorithm="Tarjan_Boost", sort=False):
                         u1, u2 = edge_stack.pop()
                     new_block.add(u1)
                     if sort:
-                        this_block = sorted(new_block)
+                        this_block = sorted(new_block, key=key)
                     else:
                         this_block = list(new_block)
                     blocks.append(this_block)
@@ -546,9 +643,8 @@ def blocks_and_cut_vertices(G, algorithm="Tarjan_Boost", sort=False):
                         start_already_seen = True
 
     if sort:
-        return blocks, sorted(cut_vertices)
-    else:
-        return blocks, list(cut_vertices)
+        return blocks, sorted(cut_vertices, key=key)
+    return blocks, list(cut_vertices)
 
 
 def blocks_and_cuts_tree(G):
@@ -1042,7 +1138,9 @@ def edge_connectivity(G,
 
         sage: g = graphs.PetersenGraph()
         sage: edge_connectivity((2 * g), vertices=True)
-        [0, [], [[0, 1, 2, 3, 4, 5, 6, 7, 8, 9], [10, 11, 12, 13, 14, 15, 16, 17, 18, 19]]]
+        [0, [], [{0, 1, 2, 3, 4, 5, 6, 7, 8, 9}, {10, 11, 12, 13, 14, 15, 16, 17, 18, 19}]]
+        sage: edge_connectivity(Graph(), vertices=True)
+        [0, [], [{}, {}]]
 
     If ``G`` is not a Sage graph, an error is raised::
 
@@ -1078,7 +1176,7 @@ def edge_connectivity(G,
         if value_only:
             return 0
         elif vertices:
-            return [0, [], [[], []]]
+            return [0, [], [{}, {}]]
         else:
             return [0, []]
 
@@ -1101,9 +1199,9 @@ def edge_connectivity(G,
                 b = set(H).difference(a)
                 val.append([a, b])
             else:
-                val.append(connected_components(H))
+                val.append([set(c) for c in connected_components(H, sort=False)])
         elif vertices:
-            val.append(connected_components(G))
+            val.append([set(c) for c in connected_components(G, sort=False)])
 
         return val
 
@@ -1190,13 +1288,13 @@ def edge_connectivity(G,
         val.append(edges)
 
         if vertices:
-            a = []
-            b = []
+            a = {}
+            b = {}
             for v in g:
                 if in_set[0, v]:
-                    a.append(v)
+                    a.add(v)
                 else:
-                    b.append(v)
+                    b.add(v)
             val.append([a, b])
 
         return val
@@ -2125,7 +2223,7 @@ def cleave(G, cut_vertices=None, virtual_edges=True, solver=None, verbose=0,
 
     H = G.copy(immutable=False)
     H.delete_vertices(cut_vertices)
-    CC = H.connected_components()
+    CC = H.connected_components(sort=False)
     if len(CC) == 1:
         raise ValueError("the set cut_vertices is not a vertex cut of the graph")
 
diff --git a/src/sage/graphs/distances_all_pairs.pyx b/src/sage/graphs/distances_all_pairs.pyx
index 2451de506fd..46baaf8ca73 100644
--- a/src/sage/graphs/distances_all_pairs.pyx
+++ b/src/sage/graphs/distances_all_pairs.pyx
@@ -2572,7 +2572,7 @@ def antipodal_graph(G):
 
     if not G.is_connected():
         import itertools
-        CC = G.connected_components()
+        CC = G.connected_components(sort=False)
         for c1, c2 in itertools.combinations(CC, 2):
             A.add_edges(itertools.product(c1, c2))
         return A
diff --git a/src/sage/graphs/generators/world_map.py b/src/sage/graphs/generators/world_map.py
index 42d0a95372f..76c32335f60 100644
--- a/src/sage/graphs/generators/world_map.py
+++ b/src/sage/graphs/generators/world_map.py
@@ -95,7 +95,7 @@ def AfricaMap(continental=False, year=2018):
     G = Graph(common_border, format='dict_of_lists')
 
     if continental:
-        G = G.subgraph(G.connected_component_containing_vertex('Central Africa'))
+        G = G.subgraph(G.connected_component_containing_vertex('Central Africa', sort=False))
         G.name(new="Continental Africa Map")
     else:
         G.add_vertices(no_land_border)
@@ -172,7 +172,7 @@ def EuropeMap(continental=False, year=2018):
     G = Graph(common_border, format='dict_of_lists')
 
     if continental:
-        G = G.subgraph(G.connected_component_containing_vertex('Austria'))
+        G = G.subgraph(G.connected_component_containing_vertex('Austria', sort=False))
         G.name(new="Continental Europe Map")
     else:
         G.add_vertices(no_land_border)
@@ -308,7 +308,7 @@ def WorldMap():
         True
         sage: g.gps_coordinates["Bolivia"]
         [[17, 'S'], [65, 'W']]
-        sage: sorted(g.connected_component_containing_vertex('Ireland'))
+        sage: g.connected_component_containing_vertex('Ireland', sort=True)
         ['Ireland', 'United Kingdom']
 
     TESTS:
diff --git a/src/sage/graphs/generic_graph.py b/src/sage/graphs/generic_graph.py
index c3f1b940b5c..bf5f79bcc4c 100644
--- a/src/sage/graphs/generic_graph.py
+++ b/src/sage/graphs/generic_graph.py
@@ -6528,7 +6528,7 @@ def num_faces(self, embedding=None):
             else:
                 if self.is_planar():
                     # We use Euler's formula: V-E+F-C=1
-                    C = len(self.connected_components())
+                    C = self.connected_components_number()
                     return self.size() - self.order() + C + 1
                 else:
                     raise ValueError("no embedding is provided and the graph is not planar")
@@ -6735,7 +6735,7 @@ def steiner_tree(self, vertices, weighted=False, solver=None, verbose=0,
 
         # Can the problem be solved ? Are all the vertices in the same
         # connected component ?
-        cc = g.connected_component_containing_vertex(vertices[0])
+        cc = g.connected_component_containing_vertex(vertices[0], sort=False)
         if any(v not in cc for v in vertices):
             from sage.categories.sets_cat import EmptySetError
             raise EmptySetError("the given vertices do not all belong to the "
@@ -9141,6 +9141,12 @@ def feedback_vertex_set(self, value_only=False, solver=None, verbose=0,
             Traceback (most recent call last):
             ...
             ValueError: the only implementation available for undirected graphs is with constraint_generation set to True
+
+        :trac:`35889` is fixed::
+
+            sage: G = Graph([('A', 1)])
+            sage: G.feedback_vertex_set()
+            []
         """
         if not constraint_generation and not self.is_directed():
             raise ValueError("the only implementation available for "
@@ -14751,7 +14757,7 @@ def is_gallai_tree(self):
         a special vertex `-1` is a ''star-shaped'' Gallai tree::
 
             sage: g = 8 * graphs.CompleteGraph(6)
-            sage: g.add_edges([(-1, c[0]) for c in g.connected_components()])
+            sage: g.add_edges([(-1, c[0]) for c in g.connected_components(sort=False)])
             sage: g.is_gallai_tree()
             True
 
@@ -19368,7 +19374,7 @@ def transitive_reduction(self):
                 return Graph(self.min_spanning_tree(weight_function=lambda e: 1))
             else:
                 G = Graph(list(self))
-                for cc in self.connected_components():
+                for cc in self.connected_components(sort=False):
                     if len(cc) > 1:
                         edges = self.subgraph(cc).min_spanning_tree(weight_function=lambda e: 1)
                         G.add_edges(edges)
diff --git a/src/sage/graphs/graph.py b/src/sage/graphs/graph.py
index 3113fe0af38..e2c86f69d40 100644
--- a/src/sage/graphs/graph.py
+++ b/src/sage/graphs/graph.py
@@ -1634,9 +1634,9 @@ def is_forest(self, certificate=False, output='vertex'):
             sage: g.is_forest(certificate=True)
             (True, None)
             sage: (2*g + graphs.PetersenGraph() + g).is_forest(certificate=True)
-            (False, [68, 66, 69, 67, 65])
+            (False, [64, 69, 67, 65, 60])
         """
-        connected_components = self.connected_components()
+        connected_components = self.connected_components(sort=False)
         number_of_connected_components = len(connected_components)
         isit = (self.order() ==
                 self.size() + number_of_connected_components)
@@ -8910,7 +8910,7 @@ def perfect_matchings(self, labels=False):
         if not self:
             yield []
             return
-        if self.order() % 2 or any(len(cc) % 2 for cc in self.connected_components()):
+        if self.order() % 2 or any(len(cc) % 2 for cc in self.connected_components(sort=False)):
             return
 
         def rec(G):
@@ -9164,7 +9164,7 @@ def effective_resistance(self, i, j, *, base_ring=None):
 
         self._scream_if_not_simple()
         if not self.is_connected():
-            connected_i = self.connected_component_containing_vertex(i)
+            connected_i = self.connected_component_containing_vertex(i, sort=False)
             if j in connected_i:
                 component = self.subgraph(connected_i)
                 return component.effective_resistance(i, j)
diff --git a/src/sage/graphs/graph_coloring.pyx b/src/sage/graphs/graph_coloring.pyx
index 0b34088ef9c..60abd45d4d1 100644
--- a/src/sage/graphs/graph_coloring.pyx
+++ b/src/sage/graphs/graph_coloring.pyx
@@ -602,7 +602,7 @@ def vertex_coloring(g, k=None, value_only=False, hex_colors=False, solver=None,
         # by the test of degeneracy (as previously).
         if not g.is_connected():
             if value_only:
-                for component in g.connected_components():
+                for component in g.connected_components(sort=False):
                     tmp = vertex_coloring(g.subgraph(component), k=k,
                                           value_only=value_only,
                                           hex_colors=hex_colors,
@@ -612,7 +612,7 @@ def vertex_coloring(g, k=None, value_only=False, hex_colors=False, solver=None,
                         return False
                 return True
             colorings = []
-            for component in g.connected_components():
+            for component in g.connected_components(sort=False):
                 tmp = vertex_coloring(g.subgraph(component), k=k,
                                       value_only=value_only,
                                       hex_colors=False,
diff --git a/src/sage/graphs/graph_decompositions/clique_separators.pyx b/src/sage/graphs/graph_decompositions/clique_separators.pyx
index 73a26a6615c..b4f92db5172 100644
--- a/src/sage/graphs/graph_decompositions/clique_separators.pyx
+++ b/src/sage/graphs/graph_decompositions/clique_separators.pyx
@@ -423,7 +423,7 @@ def atoms_and_clique_separators(G, tree=False, rooted_tree=False, separators=Fal
     if not G.is_connected():
         from sage.graphs.graph import Graph
 
-        for cc in G.connected_components():
+        for cc in G.connected_components(sort=False):
             g = Graph([cc, G.edge_boundary(cc, cc, False, False)],
                       format='vertices_and_edges',
                       loops=True, multiedges=True)
diff --git a/src/sage/graphs/graph_decompositions/graph_products.pyx b/src/sage/graphs/graph_decompositions/graph_products.pyx
index fa7fdb2e0c8..7042eabbfaf 100644
--- a/src/sage/graphs/graph_decompositions/graph_products.pyx
+++ b/src/sage/graphs/graph_decompositions/graph_products.pyx
@@ -309,7 +309,7 @@ def is_cartesian_product(g, certificate=False, relabeling=False):
 
     # Gathering the connected components, relabeling the vertices on-the-fly
     edges = [[(int_to_vertex[u], int_to_vertex[v]) for u, v in cc]
-             for cc in h.connected_components()]
+             for cc in h.connected_components(sort=False)]
 
     # Only one connected component ?
     if len(edges) == 1:
@@ -320,7 +320,7 @@ def is_cartesian_product(g, certificate=False, relabeling=False):
     for cc in edges:
         tmp = Graph()
         tmp.add_edges(cc)
-        factors.append(tmp.subgraph(vertices=tmp.connected_components()[0]))
+        factors.append(tmp.subgraph(vertices=tmp.connected_components(sort=False)[0]))
 
     # Computing the product of these graphs
     answer = factors[0]
diff --git a/src/sage/graphs/graph_decompositions/modular_decomposition.py b/src/sage/graphs/graph_decompositions/modular_decomposition.py
index 794c1d95d7b..230525e00e1 100644
--- a/src/sage/graphs/graph_decompositions/modular_decomposition.py
+++ b/src/sage/graphs/graph_decompositions/modular_decomposition.py
@@ -451,7 +451,7 @@ def gamma_classes(graph):
     pieces = DisjointSet(frozenset(e) for e in graph.edge_iterator(labels=False))
     for v in graph:
         neighborhood = graph.subgraph(vertices=graph.neighbors(v))
-        for component in neighborhood.complement().connected_components():
+        for component in neighborhood.complement().connected_components(sort=False):
             v1 = component[0]
             e = frozenset([v1, v])
             for vi in component[1:]:
@@ -625,7 +625,7 @@ def habib_maurer_algorithm(graph, g_classes=None):
     elif not graph.is_connected():
         root = create_parallel_node()
         root.children = [habib_maurer_algorithm(graph.subgraph(vertices=sg), g_classes)
-                         for sg in graph.connected_components()]
+                         for sg in graph.connected_components(sort=False)]
         return root
 
     g_comp = graph.complement()
@@ -650,7 +650,7 @@ def habib_maurer_algorithm(graph, g_classes=None):
 
     root = create_series_node()
     root.children = [habib_maurer_algorithm(graph.subgraph(vertices=sg), g_classes)
-                     for sg in g_comp.connected_components()]
+                     for sg in g_comp.connected_components(sort=False)]
     return root
 
 
diff --git a/src/sage/graphs/graph_decompositions/tree_decomposition.pyx b/src/sage/graphs/graph_decompositions/tree_decomposition.pyx
index fb9b217a49d..1e364277c8c 100644
--- a/src/sage/graphs/graph_decompositions/tree_decomposition.pyx
+++ b/src/sage/graphs/graph_decompositions/tree_decomposition.pyx
@@ -1117,7 +1117,7 @@ cdef class TreelengthConnected:
 
                 # Removing v may have disconnected cc. We iterate on its
                 # connected components
-                for _cci in g.subgraph(ccv).connected_components():
+                for _cci in g.subgraph(ccv).connected_components(sort=False):
                     cci = frozenset(_cci)
 
                     # The recursive subcalls. We remove on-the-fly the vertices
diff --git a/src/sage/graphs/partial_cube.py b/src/sage/graphs/partial_cube.py
index f7bb3e25173..0ec81f6b200 100644
--- a/src/sage/graphs/partial_cube.py
+++ b/src/sage/graphs/partial_cube.py
@@ -342,7 +342,7 @@ def is_partial_cube(G, certificate=False):
 
         # Map vertices to components of labeled-edge graph
         component = {}
-        for i, SCC in enumerate(labeled.connected_components()):
+        for i, SCC in enumerate(labeled.connected_components(sort=False)):
             for v in SCC:
                 component[v] = i
 
diff --git a/src/sage/graphs/spanning_tree.pyx b/src/sage/graphs/spanning_tree.pyx
index 179625ab743..08f441d90a8 100644
--- a/src/sage/graphs/spanning_tree.pyx
+++ b/src/sage/graphs/spanning_tree.pyx
@@ -1178,8 +1178,8 @@ def spanning_trees(g, labels=False):
             # e=xy links the CC (connected component) of forest containing x
             # with the CC containing y. Any other edge which does that cannot be
             # added to forest anymore, and B is the list of them
-            c1 = forest.connected_component_containing_vertex(e[0])
-            c2 = forest.connected_component_containing_vertex(e[1])
+            c1 = forest.connected_component_containing_vertex(e[0], sort=False)
+            c2 = forest.connected_component_containing_vertex(e[1], sort=False)
             G.delete_edge(e)
             B = G.edge_boundary(c1, c2, sort=False)
             G.add_edge(e)
diff --git a/src/sage/graphs/traversals.pyx b/src/sage/graphs/traversals.pyx
index 9ebc253b6f9..a565d46fdcb 100644
--- a/src/sage/graphs/traversals.pyx
+++ b/src/sage/graphs/traversals.pyx
@@ -2009,7 +2009,7 @@ def maximum_cardinality_search_M(G, initial_vertex=None):
         ....:     if len(X) < k - 1:
         ....:         raise ValueError("something goes wrong")
         sage: G = graphs.RandomGNP(10, .2)
-        sage: cc = G.connected_components()
+        sage: cc = G.connected_components(sort=False)
         sage: _, _, X = G.maximum_cardinality_search_M()
         sage: len(X) >= len(cc) - 1
         True
diff --git a/src/sage/graphs/tutte_polynomial.py b/src/sage/graphs/tutte_polynomial.py
index 8f003084619..6794c960805 100644
--- a/src/sage/graphs/tutte_polynomial.py
+++ b/src/sage/graphs/tutte_polynomial.py
@@ -317,7 +317,7 @@ def find_ear(g):
                                in g.degree_iterator(labels=True)
                                if degree == 2]
         subgraph = g.subgraph(degree_two_vertices)
-        for component in subgraph.connected_components():
+        for component in subgraph.connected_components(sort=False):
             edges = g.edges_incident(vertices=component, labels=True)
             all_vertices = sorted(set(sum([e[:2] for e in edges], ())))
             if len(all_vertices) < 3:
diff --git a/src/sage/matroids/linear_matroid.pyx b/src/sage/matroids/linear_matroid.pyx
index 30bca7bea29..11f6b4a2656 100644
--- a/src/sage/matroids/linear_matroid.pyx
+++ b/src/sage/matroids/linear_matroid.pyx
@@ -578,7 +578,7 @@ cdef class LinearMatroid(BasisExchangeMatroid):
             sage: M.representation(lift_map=lift_map('sru'))                            # needs sage.rings.finite_rings
             [     1      0      0      1      0      1      1      1]
             [     0      1      0 -z + 1      1      0      0      1]
-            [     0      0      1      0     -z      z      1      0]
+            [     0      0      1      0      1      -1 z - 1      0]
         """
         cdef LeanMatrix A
         if order is None:
diff --git a/src/sage/matroids/matroid.pyx b/src/sage/matroids/matroid.pyx
index 3caac5bea62..dcee87ea49c 100644
--- a/src/sage/matroids/matroid.pyx
+++ b/src/sage/matroids/matroid.pyx
@@ -6114,8 +6114,8 @@ cdef class Matroid(SageObject):
         entries = [(e, f, (e, f)) for e in basis for f in self._fundamental_cocircuit(basis, e).difference([e])]
         G = Graph(entries)
         T = set()
-        for C in G.connected_components():
-            T.update(G.subgraph(C).min_spanning_tree())
+        for C in G.connected_components_subgraphs():
+            T.update(C.min_spanning_tree())
         for edge in T:
             e,f = edge[2]
             A.set(bdx[e],idx[f], 1)
diff --git a/src/sage/matroids/utilities.py b/src/sage/matroids/utilities.py
index 6351c536413..a8612dab29f 100644
--- a/src/sage/matroids/utilities.py
+++ b/src/sage/matroids/utilities.py
@@ -392,8 +392,8 @@ def spanning_forest(M):
         G.add_edge(x + m, y)
     T = []
     # find spanning tree in each component
-    for component in G.connected_components():
-        spanning_tree = kruskal(G.subgraph(component))
+    for component in G.connected_components_subgraphs():
+        spanning_tree = kruskal(component)
         for (x, y, z) in spanning_tree:
             if x < m:
                 t = x
@@ -546,9 +546,9 @@ def lift_cross_ratios(A, lift_map=None):
         [0 6 3 6 0]
         sage: Z = lift_cross_ratios(A, to_sixth_root_of_unity)                          # needs sage.rings.finite_rings sage.rings.number_field
         sage: Z                                                                         # needs sage.rings.finite_rings sage.rings.number_field
-        [ 1  0  1  1  1]
-        [ 1  1  0  0  z]
-        [ 0 -1  z  1  0]
+        [     1      0      1      1      1]
+        [-z + 1      1      0      0      1]
+        [     0     -1      1 -z + 1      0]
         sage: M = LinearMatroid(reduced_matrix=A)
         sage: sorted(M.cross_ratios())
         [3, 5]
@@ -573,8 +573,8 @@ def lift_cross_ratios(A, lift_map=None):
     G = Graph([((r, 0), (c, 1), (r, c)) for r, c in A.nonzero_positions()])
     # write the entries of (a scaled version of) A as products of cross ratios of A
     T = set()
-    for C in G.connected_components():
-        T.update(G.subgraph(C).min_spanning_tree())
+    for C in G.connected_components_subgraphs():
+        T.update(C.min_spanning_tree())
     # - fix a tree of the support graph G to units (= empty dict, product of 0 terms)
     F = {entry[2]: dict() for entry in T}
     W = set(G.edge_iterator()) - set(T)
diff --git a/src/sage/schemes/curves/zariski_vankampen.py b/src/sage/schemes/curves/zariski_vankampen.py
index 751ed6502e0..2edb29f365a 100644
--- a/src/sage/schemes/curves/zariski_vankampen.py
+++ b/src/sage/schemes/curves/zariski_vankampen.py
@@ -815,8 +815,9 @@ def geometric_basis(G, E, p):
             q = EC[-i]
             connecting_path = list(reversed(EC[-i:]))
             break
+    I_cc_q = set(I.connected_component_containing_vertex(q, sort=False))
     distancequotients = [(E.distance(q, v)**2 / I.distance(q, v), v) for v in E
-                         if v in I.connected_component_containing_vertex(q) and not v == q]
+                         if v in I_cc_q and not v == q]
     r = max(distancequotients)[1]
     cutpath = I.shortest_path(q, r)
     Gcut = copy(G)
diff --git a/src/sage/tests/books/computational-mathematics-with-sagemath/lp_doctest.py b/src/sage/tests/books/computational-mathematics-with-sagemath/lp_doctest.py
index f3aa2201ac8..da89f065214 100644
--- a/src/sage/tests/books/computational-mathematics-with-sagemath/lp_doctest.py
+++ b/src/sage/tests/books/computational-mathematics-with-sagemath/lp_doctest.py
@@ -221,7 +221,7 @@
 Sage example in ./lp.tex, line 906::
 
   sage: while not h.is_connected():
-  ....:     S = h.connected_components()[0]
+  ....:     S = h.connected_components(sort=False)[0]
   ....:     p.add_constraint(
   ....:            p.sum( B(u,v) for u,v
   ....:                 in g.edge_boundary(S, labels = False))
diff --git a/src/sage/topology/simplicial_complex.py b/src/sage/topology/simplicial_complex.py
index 20b89bb61c8..91f7deb1486 100644
--- a/src/sage/topology/simplicial_complex.py
+++ b/src/sage/topology/simplicial_complex.py
@@ -4078,7 +4078,7 @@ def connected_component(self, simplex=None):
             v = self.vertices()[0]
         else:
             v = simplex[0]
-        vertices = self.graph().connected_component_containing_vertex(v)
+        vertices = self.graph().connected_component_containing_vertex(v, sort=False)
         facets = [f for f in self.facets() if f.is_face(Simplex(vertices))]
         return SimplicialComplex(facets)
 
diff --git a/src/sage/topology/simplicial_set_constructions.py b/src/sage/topology/simplicial_set_constructions.py
index 316d30bb1a1..b2b6bbbae7c 100644
--- a/src/sage/topology/simplicial_set_constructions.py
+++ b/src/sage/topology/simplicial_set_constructions.py
@@ -1499,7 +1499,7 @@ def __init__(self, maps=None, vertex_name=None):
             for x in domain.n_cells(n):
                 edges.extend([[(x, -1), (f(x), i)] for i, f in enumerate(maps)])
             G = Graph([vertices, edges], format='vertices_and_edges')
-            data[n] = [set(_) for _ in G.connected_components()]
+            data[n] = [set(_) for _ in G.connected_components(sort=False)]
         # data is now a dictionary indexed by dimension, and data[n]
         # consists of sets of n-simplices of the domain and the
         # codomains, each set an equivalence class of n-simplices