Different behaviour of a mesh coming from generate_mesh() resp. UnitSquareMesh()?

As these questions are directly related to interface changes in pygmsh, the two last questions are best to ask in the issue tracker for GitHub - nschloe/pygmsh: Gmsh for Python.
In general, you can remove the z-coordinate from any 2D mesh when you are converting it with meshio
(I.e. if the read in points with meshio are of the shape (num_points, 3), remove the last column before saving it to xdmf).
I personally stopped using pygmsh after v 6.1.1, as it now is more or less the same as the GMSH python API, which I have referenced in previous posts on this topic.

My suggested pipeline for new users is the following:

  1. Generate an msh file using either GMSH graphical interface or GMSH python interface.
  2. Convert the mesh from msh to XDMF using meshio.

Optionally, replace gmsh with pygmsh or any other mesh generator that can produce a file that is readable by meshio.

There are various mesh formats available for representing unstructured meshes. 
meshio can read and write all of the following and smoothly converts between them:

    Abaqus, ANSYS msh, AVS-UCD, CGNS, DOLFIN XML, Exodus, FLAC3D, H5M, Kratos/MDPA, Medit, MED/Salome, Nastran (bulk data), Neuroglancer precomputed format, Gmsh (format versions 2.2, 4.0, and 4.1), OBJ, OFF, PERMAS, PLY, STL, Tecplot .dat, TetGen .node/.ele, SVG (2D only, output only), SU2, UGRID, VTK, VTU, WKT (TIN), XDMF.