We present full QCD correlator data and corresponding reconstructed spectral functions in the pseudoscalar channel. Correlators are obtained using clover-improved Wilson fermions on $N_f=2+1$ HISQ lattices. We use gradient flow to check whether it reduces cut-off and mixed action effects. Valence quark masses are tuned to their physical values by comparing the mass spectrum obtained from the lattice QCD with experimental values at each flow time. For the spectral reconstruction, we use models based on perturbative spectral functions from different frequency regions like resummed thermal contributions around the threshold from pNRQCD and vacuum contributions well above the threshold. We show preliminary results of the reconstructed spectral function obtained for the first time in our study for full QCD. In addition, we compare the results with the previous continuum extrapolated results in the quenched approximation.