Among the microscopic approaches to nuclear many-body problem, nuclear density functional theory is nowadays the only one that can be applied for nuclei throughout the nuclear chart. In the past years, we have established a method of Multi-Reference Density Functional Theory (MR-DFT) by introducing generator coordinate method (GCM), particle number and 3D-angular momentum projections into the covariant DFT for nuclei allowing for triaxial shapes. The MR-DFT provides us a state-of-the-art and full microscopic description of nuclear low-lying states with a universal energy density functional. In this talk, I am going to introduce this beyond-mean-field method briefly and focus more on its recent applications to several selected topics, including the collectivity evolution in N=40 isotones, triaxiality in nuclear low-lying states, dynamical correlation effects on nuclear density and transition densities, impurity effect of Lambda hyperon on nuclear collectivity, and a preliminary calculation of nuclear matrix element (NME) in double beta-decay.