Thermal and many-body localized phases are separated by a dynamical phase transition of a new kind. We analyze the distribution of off-diagonal matrix elements of local operators across the many-body localization transition (MBLT) in a disordered spin chain, and use it to characterize the breakdown of the eigenstate thermalization hypothesis and to extract the many-body Thouless energy. We find a wide critical region around the MBLT, where Thouless energy becomes smaller than the level spacing, while matrix elements show critical dependence on the energy difference. In the same region, matrix elements, viewed as amplitudes of a fictitious wave function, exhibit strong multifractality. Our findings show that the correlation length becomes larger than the accessible system sizes in a broad range of disorder strength values, and shed light on the critical behaviour of MBL systems.