So here the solution.
I had to rewrite the boundary conditions as follows:
facets_axis = mesh.locate_entities_boundary(domain, fdim, axis)
facets_bottom = mesh.locate_entities_boundary(domain, fdim, bottom)
# Locate dofs (IMPORTANT tuple syntax)
dofs_axis = fem.locate_dofs_topological(V.sub(0), fdim, facets_axis)
dofs_bottom = fem.locate_dofs_topological(V.sub(1), fdim, facets_bottom)
# Ur = 0 on axis (ID 3)
bc_axis = fem.dirichletbc(fem.Constant(domain, default_scalar_type(0.0)), dofs_axis, V.sub(0))
# Uz = 0 on bottom (ID 2)
bc_bottom = fem.dirichletbc(fem.Constant(domain,default_scalar_type(0.0)), dofs_bottom, V.sub(1))
meaning using the dofs and not directly the boundaries