The bond nature in beryllium dimer has been theoretically investigated using high-level ab initio methods. A series of ANO basis sets of increasing quality, going from sp to spdfghi contractions, has been employed, combined with HF, CAS-SCF, CISD and MRCI calculations with several different active spaces. The quality of these calculations has been checked by comparing the results with valence Full-CI calculations, performed with the same basis sets. It is shown that two quasi-degenerated partly- occupied orbitals play a crucial role in order to give a qualitatively correct description of the bond. Their nature is similar to the edge orbitals that give rise to the quasi- degenerated singlet-triplet states in longer beryllium chains.